/[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.42 by jmc, Tue Dec 8 21:43:43 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    #include "MOM_VISC.h"
85    
86  C     == Routine arguments ==  C     == Routine arguments ==
87        INTEGER bi,bj,k        INTEGER bi,bj,k
# Line 90  C     == Routine arguments == Line 100  C     == Routine arguments ==
100    
101  C     == Local variables ==  C     == Local variables ==
102        INTEGER I,J        INTEGER I,J
103    #ifdef ALLOW_NONHYDROSTATIC
104          INTEGER kp1
105    #endif
106    #ifdef ALLOW_AUTODIFF_TAMC
107          INTEGER lockey_1, lockey_2
108    #endif
109        _RL smag2fac, smag4fac        _RL smag2fac, smag4fac
110        _RL leith2fac, leith4fac        _RL leith2fac, leith4fac
111        _RL leithD2fac, leithD4fac        _RL leithD2fac, leithD4fac
112        _RL viscAhRe_max, viscA4Re_max        _RL viscAhRe_max, viscA4Re_max
113        _RL Alin,grdVrt,grdDiv, keZpt        _RL Alin,grdVrt,grdDiv, keZpt
114        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt        _RL L2,L3,L5,L2rdt,L4rdt
115        _RL Uscl,U4scl        _RL Uscl,U4scl
116        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118          _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119          _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _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 137  C     == Local variables ==
137        _RL viscAh_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138        _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139        _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140        LOGICAL calcLeith,calcSmag        LOGICAL calcLeith, calcSmag
141    
142    #ifdef ALLOW_AUTODIFF_TAMC
143              act1 = bi - myBxLo(myThid)
144              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
145              act2 = bj - myByLo(myThid)
146              max2 = myByHi(myThid) - myByLo(myThid) + 1
147              act3 = myThid - 1
148              max3 = nTx*nTy
149              act4 = ikey_dynamics - 1
150              ikey = (act1 + 1) + act2*max1
151         &                      + act3*max1*max2
152         &                      + act4*max1*max2*max3
153              lockey_1 = (ikey-1)*Nr + k
154    #endif /* ALLOW_AUTODIFF_TAMC */
155    
156    C--   Set flags which are used in this S/R and elsewhere :
157        useVariableViscosity=        useVariableViscosity=
158       &      (viscAhGrid.NE.0.)       &      (viscAhGrid.NE.0.)
159       &  .OR.(viscA4Grid.NE.0.)       &  .OR.(viscA4Grid.NE.0.)
# Line 140  C     == Local variables == Line 173  C     == Local variables ==
173       &  .OR.(viscC2leithD.NE.0.)       &  .OR.(viscC2leithD.NE.0.)
174       &  .OR.(viscC2smag.NE.0.)       &  .OR.(viscC2smag.NE.0.)
175    
       IF ((harmonic).and.(viscAhremax.ne.0.)) THEN  
         viscAhre_max=sqrt(2. _d 0)/viscAhRemax  
       ELSE  
         viscAhre_max=0. _d 0  
       ENDIF  
   
176        biharmonic=        biharmonic=
177       &      (viscA4.NE.0.)       &      (viscA4.NE.0.)
178       &  .OR.(viscA4D.NE.0.)       &  .OR.(viscA4D.NE.0.)
# Line 155  C     == Local variables == Line 182  C     == Local variables ==
182       &  .OR.(viscC4leithD.NE.0.)       &  .OR.(viscC4leithD.NE.0.)
183       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
184    
185        IF ((biharmonic).and.(viscA4remax.ne.0.)) THEN        IF (useVariableViscosity) THEN
186          viscA4re_max=0.125 _d 0*sqrt(2. _d 0)/viscA4Remax  C---- variable viscosity :
187        ELSE  
188          viscA4re_max=0. _d 0         IF ((harmonic).AND.(viscAhReMax.NE.0.)) THEN
189        ENDIF          viscAhRe_max=SQRT(2. _d 0)/viscAhReMax
190           ELSE
191            viscAhRe_max=0. _d 0
192           ENDIF
193    
194           IF ((biharmonic).AND.(viscA4ReMax.NE.0.)) THEN
195            viscA4Re_max=0.125 _d 0*SQRT(2. _d 0)/viscA4ReMax
196           ELSE
197            viscA4Re_max=0. _d 0
198           ENDIF
199    
200        calcleith=         calcLeith=
201       &      (viscC2leith.NE.0.)       &      (viscC2leith.NE.0.)
202       &  .OR.(viscC2leithD.NE.0.)       &  .OR.(viscC2leithD.NE.0.)
203       &  .OR.(viscC4leith.NE.0.)       &  .OR.(viscC4leith.NE.0.)
204       &  .OR.(viscC4leithD.NE.0.)       &  .OR.(viscC4leithD.NE.0.)
205    
206        calcsmag=         calcSmag=
207       &      (viscC2smag.NE.0.)       &      (viscC2smag.NE.0.)
208       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
209    
210        IF (deltaTmom.NE.0.) THEN         IF (calcSmag) THEN
        recip_dt=1. _d 0/deltaTmom  
       ELSE  
        recip_dt=0. _d 0  
       ENDIF  
   
       IF (calcsmag) THEN  
211          smag2fac=(viscC2smag/pi)**2          smag2fac=(viscC2smag/pi)**2
212          smag4fac=0.125 _d 0*(viscC4smag/pi)**2          smag4fac=0.125 _d 0*(viscC4smag/pi)**2
213        ELSE         ELSE
214          smag2fac=0. _d 0          smag2fac=0. _d 0
215          smag4fac=0. _d 0          smag4fac=0. _d 0
216        ENDIF         ENDIF
217    
218        IF (calcleith) THEN         IF (calcLeith) THEN
219          IF (useFullLeith) THEN          IF (useFullLeith) THEN
220           leith2fac =(viscC2leith /pi)**6           leith2fac =(viscC2leith /pi)**6
221           leithD2fac=(viscC2leithD/pi)**6           leithD2fac=(viscC2leithD/pi)**6
# Line 197  C     == Local variables == Line 227  C     == Local variables ==
227           leith4fac =0.125 _d 0*(viscC4leith /pi)**3           leith4fac =0.125 _d 0*(viscC4leith /pi)**3
228           leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3           leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3
229          ENDIF          ENDIF
230        ELSE         ELSE
231          leith2fac=0. _d 0          leith2fac=0. _d 0
232          leith4fac=0. _d 0          leith4fac=0. _d 0
233          leithD2fac=0. _d 0          leithD2fac=0. _d 0
234          leithD4fac=0. _d 0          leithD4fac=0. _d 0
235        ENDIF         ENDIF
236    
237  C     - Viscosity  #ifdef ALLOW_AUTODIFF_TAMC
238        IF (useVariableViscosity) THEN  cphtest       IF ( calcLeith .OR. calcSmag ) THEN
239    cphtest        STOP 'calcLeith or calcSmag not implemented for ADJOINT'
240    cphtest       ENDIF
241    #endif
242           DO j=1-Oly,sNy+Oly
243            DO i=1-Olx,sNx+Olx
244              viscAh_D(i,j)=viscAhD
245              viscAh_Z(i,j)=viscAhZ
246              viscA4_D(i,j)=viscA4D
247              viscA4_Z(i,j)=viscA4Z
248    c
249              visca4_zsmg(i,j) = 0. _d 0
250              viscah_zsmg(i,j) = 0. _d 0
251    c
252              viscAh_Dlth(i,j) = 0. _d 0
253              viscA4_Dlth(i,j) = 0. _d 0
254              viscAh_DlthD(i,j)= 0. _d 0
255              viscA4_DlthD(i,j)= 0. _d 0
256    c
257              viscAh_DSmg(i,j) = 0. _d 0
258              viscA4_DSmg(i,j) = 0. _d 0
259    c
260              viscAh_ZLth(i,j) = 0. _d 0
261              viscA4_ZLth(i,j) = 0. _d 0
262              viscAh_ZLthD(i,j)= 0. _d 0
263              viscA4_ZLthD(i,j)= 0. _d 0
264            ENDDO
265           ENDDO
266    
267  C      horizontal gradient of horizontal divergence:  C-     Initialise to zero gradient of vorticity & divergence:
268         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
269          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
270            divDx(i,j) = 0.            divDx(i,j) = 0.
271            divDy(i,j) = 0.            divDy(i,j) = 0.
272              vrtDx(i,j) = 0.
273              vrtDy(i,j) = 0.
274          ENDDO          ENDDO
275         ENDDO         ENDDO
276         IF (calcleith) THEN  
277           IF (calcLeith) THEN
278    C      horizontal gradient of horizontal divergence:
279    
280  C-       gradient in x direction:  C-       gradient in x direction:
281  #ifndef ALLOW_AUTODIFF_TAMC  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
282           IF (useCubedSphereExchange) THEN           IF (useCubedSphereExchange) THEN
283  C        to compute d/dx(hDiv), fill corners with appropriate values:  C        to compute d/dx(hDiv), fill corners with appropriate values:
284             CALL FILL_CS_CORNER_TR_RL( .TRUE., hDiv, bi,bj, myThid )             CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
285         &                                hDiv, bi,bj, myThid )
286           ENDIF           ENDIF
287  #endif  cph-exch2#endif
288           DO j=2-Oly,sNy+Oly-1           DO j=2-Oly,sNy+Oly-1
289            DO i=2-Olx,sNx+Olx-1            DO i=2-Olx,sNx+Olx-1
290              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 292  C        to compute d/dx(hDiv), fill cor
292           ENDDO           ENDDO
293    
294  C-       gradient in y direction:  C-       gradient in y direction:
295  #ifndef ALLOW_AUTODIFF_TAMC  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
296           IF (useCubedSphereExchange) THEN           IF (useCubedSphereExchange) THEN
297  C        to compute d/dy(hDiv), fill corners with appropriate values:  C        to compute d/dy(hDiv), fill corners with appropriate values:
298             CALL FILL_CS_CORNER_TR_RL(.FALSE., hDiv, bi,bj, myThid )             CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
299         &                                hDiv, bi,bj, myThid )
300           ENDIF           ENDIF
301  #endif  cph-exch2#endif
302           DO j=2-Oly,sNy+Oly-1           DO j=2-Oly,sNy+Oly-1
303            DO i=2-Olx,sNx+Olx-1            DO i=2-Olx,sNx+Olx-1
304              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)
305            ENDDO            ENDDO
306           ENDDO           ENDDO
307    
308    C      horizontal gradient of vertical vorticity:
309    C-       gradient in x direction:
310             DO j=2-Oly,sNy+Oly
311              DO i=2-Olx,sNx+Olx-1
312                vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
313         &                  *recip_DXG(i,j,bi,bj)
314         &                  *maskS(i,j,k,bi,bj)
315              ENDDO
316             ENDDO
317    C-       gradient in y direction:
318             DO j=2-Oly,sNy+Oly-1
319              DO i=2-Olx,sNx+Olx
320                vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
321         &                  *recip_DYG(i,j,bi,bj)
322         &                  *maskW(i,j,k,bi,bj)
323              ENDDO
324             ENDDO
325    
326         ENDIF         ENDIF
327    
328         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
329          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
330  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
331    
332  C These are (powers of) length scales  #ifdef    ALLOW_AUTODIFF_TAMC
333           IF (useAreaViscLength) THEN  # ifndef AUTODIFF_DISABLE_LEITH
334            L2=rA(i,j,bi,bj)             lockey_2 = i+olx + (sNx+2*olx)*(j+oly-1)
335           ELSE       &                      + (sNx+2*olx)*(sNy+2*oly)*(lockey_1-1)
336            L2=2. _d 0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))  CADJ STORE viscA4_ZSmg(i,j)
337           ENDIF  CADJ &     = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte
338           L3=(L2**1.5)  CADJ STORE viscAh_ZSmg(i,j)
339           L4=(L2**2)  CADJ &     = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte
340           L5=(L2**2.5)  # endif
341    #endif    /* ALLOW_AUTODIFF_TAMC */
342           L2rdt=0.25 _d 0*recip_dt*L2  
343    C These are (powers of) length scales
344           IF (useAreaViscLength) THEN           L2 = L2_D(i,j,bi,bj)
345            L4rdt=0.125 _d 0*recip_dt*rA(i,j,bi,bj)**2           L2rdt = 0.25 _d 0*recip_dt*L2
346           ELSE           L3 = L3_D(i,j,bi,bj)
347            L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4           L4rdt = L4rdt_D(i,j,bi,bj)
348       &                            +recip_DYF(I,J,bi,bj)**4)           L5 = (L2*L3)
      &                   +8. _d 0*((recip_DXF(I,J,bi,bj)  
      &                             *recip_DYF(I,J,bi,bj))**2) )  
          ENDIF  
349    
350    #ifndef AUTODIFF_DISABLE_REYNOLDS_SCALE
351  C Velocity Reynolds Scale  C Velocity Reynolds Scale
352           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
353             Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max             Uscl=SQRT(KE(i,j)*L2)*viscAhRe_max
354           ELSE           ELSE
355             Uscl=0.             Uscl=0.
356           ENDIF           ENDIF
357           IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN           IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
358             U4scl=sqrt(KE(i,j))*L3*viscA4Re_max             U4scl=SQRT(KE(i,j))*L3*viscA4Re_max
359           ELSE           ELSE
360             U4scl=0.             U4scl=0.
361           ENDIF           ENDIF
362    #endif /* ndef AUTODIFF_DISABLE_REYNOLDS_SCALE */
363    
364           IF (useFullLeith.and.calcleith) THEN  cph-leith#ifndef ALLOW_AUTODIFF_TAMC
365    #ifndef AUTODIFF_DISABLE_LEITH
366             IF (useFullLeith.AND.calcLeith) THEN
367  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
368            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
369       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
370       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
371       &     +((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)  
372    
373  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
374  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 378  C Using it in Leith serves to damp insta
378       &                        + divDy(i,j)*divDy(i,j) )  )       &                        + divDy(i,j)*divDy(i,j) )  )
379    
380            viscAh_DLth(i,j)=            viscAh_DLth(i,j)=
381       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3       &     SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
382            viscA4_DLth(i,j)=            viscA4_DLth(i,j)=
383       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5       &     SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
384            viscAh_DLthd(i,j)=            viscAh_DLthd(i,j)=
385       &     sqrt(leithD2fac*grdDiv)*L3       &     SQRT(leithD2fac*grdDiv)*L3
386            viscA4_DLthd(i,j)=            viscA4_DLthd(i,j)=
387       &     sqrt(leithD4fac*grdDiv)*L5       &     SQRT(leithD4fac*grdDiv)*L5
388           ELSEIF (calcleith) THEN           ELSEIF (calcLeith) THEN
389  C but this approximation will work on cube  C but this approximation will work on cube
390  c (and differs by as much as 4X)  c (and differs by as much as 4X)
391            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)) )
392            grdVrt=max(grdVrt,            grdVrt=MAX( grdVrt, ABS(vrtDy(i+1,j)) )
393       &     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))   )  
394    
395  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
396              grdDiv=MAX( ABS(divDx(i+1,j)), ABS(divDx(i,j)) )
397              grdDiv=MAX( grdDiv, ABS(divDy(i,j+1)) )
398              grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
399    
400            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
401            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
402            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 408  c This approximation is good to the same
408            viscA4_DlthD(i,j)=0. _d 0            viscA4_DlthD(i,j)=0. _d 0
409           ENDIF           ENDIF
410    
411           IF (calcsmag) THEN           IF (calcSmag) THEN
412            viscAh_DSmg(i,j)=L2            viscAh_DSmg(i,j)=L2
413       &       *sqrt(tension(i,j)**2       &       *SQRT(tension(i,j)**2
414       &       +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
415       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))
416            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 419  c This approximation is good to the same
419            viscAh_DSmg(i,j)=0. _d 0            viscAh_DSmg(i,j)=0. _d 0
420            viscA4_DSmg(i,j)=0. _d 0            viscA4_DSmg(i,j)=0. _d 0
421           ENDIF           ENDIF
422    #endif /* AUTODIFF_DISABLE_LEITH */
423    
424  C  Harmonic on Div.u points  C  Harmonic on Div.u points
425           Alin=viscAhD+viscAhGrid*L2rdt           Alin=viscAhD+viscAhGrid*L2rdt
426       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
427           viscAh_DMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)           viscAh_DMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
428           viscAh_D(i,j)=max(viscAh_DMin(i,j),Alin)           viscAh_D(i,j)=MAX(viscAh_DMin(i,j),Alin)
429           viscAh_DMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)           viscAh_DMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
430           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))
431    
432  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
433           Alin=viscA4D+viscA4Grid*L4rdt           Alin=viscA4D+viscA4Grid*L4rdt
434       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
435           viscA4_DMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)           viscA4_DMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
436           viscA4_D(i,j)=max(viscA4_DMin(i,j),Alin)           viscA4_D(i,j)=MAX(viscA4_DMin(i,j),Alin)
437           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)           viscA4_DMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
438           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))
439    
440    #ifdef ALLOW_NONHYDROSTATIC
441    C--   Pass Viscosities to calc_gw, if constant, not necessary
442    
443             kp1  = MIN(k+1,Nr)
444    
445             IF ( k.EQ.1 ) THEN
446    C     Prepare for next level (next call)
447              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
448              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
449    
450    C     These values dont get used
451              viscAh_W(i,j,k,bi,bj)=viscAh_D(i,j)
452              viscA4_W(i,j,k,bi,bj)=viscA4_D(i,j)
453    
454             ELSEIF ( k.EQ.Nr ) THEN
455              viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)
456              viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)
457    
 CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  
 C These are (powers of) length scales  
          IF (useAreaViscLength) THEN  
           L2=rAz(i,j,bi,bj)  
458           ELSE           ELSE
459            L2=2. _d 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))  C     Prepare for next level (next call)
460           ENDIF            viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
461              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
462    
463    C     Note that previous call of this function has already added half.
464              viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)
465              viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)
466    
          L3=(L2**1.5)  
          L4=(L2**2)  
          L5=(L2**2.5)  
   
          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))  
467           ENDIF           ENDIF
468    #endif /* ALLOW_NONHYDROSTATIC */
469    
470    CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
471    C These are (powers of) length scales
472             L2 = L2_Z(i,j,bi,bj)
473             L2rdt = 0.25 _d 0*recip_dt*L2
474             L3 = L3_Z(i,j,bi,bj)
475             L4rdt = L4rdt_Z(i,j,bi,bj)
476             L5 = (L2*L3)
477    
478    #ifndef AUTODIFF_DISABLE_REYNOLDS_SCALE
479  C Velocity Reynolds Scale (Pb here at CS-grid corners !)  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
480           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
481             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))
482       &                      +(KE(i-1,j)+KE(i,j-1)) )       &                      +(KE(i-1,j)+KE(i,j-1)) )
483             IF ( keZpt.GT.0. ) THEN             IF ( keZpt.GT.0. ) THEN
484               Uscl = sqrt(keZpt*L2)*viscAhRe_max               Uscl = SQRT(keZpt*L2)*viscAhRe_max
485               U4scl= sqrt(keZpt)*L3*viscA4Re_max               U4scl= SQRT(keZpt)*L3*viscA4Re_max
486             ELSE             ELSE
487               Uscl =0.               Uscl =0.
488               U4scl=0.               U4scl=0.
# Line 395  C Velocity Reynolds Scale (Pb here at CS Line 491  C Velocity Reynolds Scale (Pb here at CS
491             Uscl =0.             Uscl =0.
492             U4scl=0.             U4scl=0.
493           ENDIF           ENDIF
494    #endif /* ndef AUTODIFF_DISABLE_REYNOLDS_SCALE */
495    
496    #ifndef AUTODIFF_DISABLE_LEITH
497  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
498           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.AND.calcLeith) THEN
499            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
500       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
501       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
502       &     +((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)  
503    
504  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
505            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 508  C This is the vector magnitude of grad(d
508       &                        + divDy(i,j)*divDy(i,j) )  )       &                        + divDy(i,j)*divDy(i,j) )  )
509    
510            viscAh_ZLth(i,j)=            viscAh_ZLth(i,j)=
511       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3       &     SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
512            viscA4_ZLth(i,j)=            viscA4_ZLth(i,j)=
513       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5       &     SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
514            viscAh_ZLthD(i,j)=            viscAh_ZLthD(i,j)=
515       &     sqrt(leithD2fac*grdDiv)*L3       &     SQRT(leithD2fac*grdDiv)*L3
516            viscA4_ZLthD(i,j)=            viscA4_ZLthD(i,j)=
517       &     sqrt(leithD4fac*grdDiv)*L5       &     SQRT(leithD4fac*grdDiv)*L5
518    
519           ELSEIF (calcleith) THEN           ELSEIF (calcLeith) THEN
520  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
521            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)) )
522            grdVrt=max(grdVrt,            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j-1)) )
523       &     abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j))   )
524            grdVrt=max(grdVrt,  
525       &     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)) )
526            grdVrt=max(grdVrt,            grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
527       &     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)) )  
528    
529            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
530            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 537  C but this approximation will work on cu
537            viscA4_ZLthD(i,j)=0. _d 0            viscA4_ZLthD(i,j)=0. _d 0
538           ENDIF           ENDIF
539    
540           IF (calcsmag) THEN           IF (calcSmag) THEN
541            viscAh_ZSmg(i,j)=L2            viscAh_ZSmg(i,j)=L2
542       &      *sqrt(strain(i,j)**2       &      *SQRT(strain(i,j)**2
543       &        +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
544       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))
545            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
546            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
547           ENDIF           ENDIF
548    #endif /* AUTODIFF_DISABLE_LEITH */
549    
550  C  Harmonic on Zeta points  C  Harmonic on Zeta points
551           Alin=viscAhZ+viscAhGrid*L2rdt           Alin=viscAhZ+viscAhGrid*L2rdt
552       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
553           viscAh_ZMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)           viscAh_ZMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
554           viscAh_Z(i,j)=max(viscAh_ZMin(i,j),Alin)           viscAh_Z(i,j)=MAX(viscAh_ZMin(i,j),Alin)
555           viscAh_ZMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)           viscAh_ZMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
556           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))
557    
558  C  BiHarmonic on Zeta points  C  BiHarmonic on Zeta points
559           Alin=viscA4Z+viscA4Grid*L4rdt           Alin=viscA4Z+viscA4Grid*L4rdt
560       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
561           viscA4_ZMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)           viscA4_ZMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
562           viscA4_Z(i,j)=max(viscA4_ZMin(i,j),Alin)           viscA4_Z(i,j)=MAX(viscA4_ZMin(i,j),Alin)
563           viscA4_ZMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)           viscA4_ZMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
564           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))
565          ENDDO          ENDDO
566         ENDDO         ENDDO
567    
568        ELSE        ELSE
569    C---- use constant viscosity (useVariableViscosity=F):
570    
571         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
572          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
573           viscAh_D(i,j)=viscAhD           viscAh_D(i,j)=viscAhD
# Line 479  C  BiHarmonic on Zeta points Line 576  C  BiHarmonic on Zeta points
576           viscA4_Z(i,j)=viscA4Z           viscA4_Z(i,j)=viscA4Z
577          ENDDO          ENDDO
578         ENDDO         ENDDO
579    
580    C---- variable/constant viscosity : end if/else block
581        ENDIF        ENDIF
582    
583  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.42
  ViewVC Help
Powered by ViewVC 1.1.22