/[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.2 by baylor, Mon Sep 19 20:46:23 2005 UTC revision 1.33 by heimbach, Fri Apr 18 18:11:43 2008 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "MOM_COMMON_OPTIONS.h"  #include "MOM_COMMON_OPTIONS.h"
5    
6    
7        SUBROUTINE MOM_CALC_VISC(        SUBROUTINE MOM_CALC_VISC(
8       I        bi,bj,k,       I        bi,bj,k,
9       O        viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,       O        viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
10       O        harmonic,biharmonic,useVariableViscosity,       O        harmonic,biharmonic,useVariableViscosity,
11       I        hDiv,vort3,tension,strain,KE,       I        hDiv,vort3,tension,strain,KE,hFacZ,
12       I        myThid)       I        myThid)
13    
14        IMPLICIT NONE        IMPLICIT NONE
15    C
16    C     Calculate horizontal viscosities (L is typical grid width)
17    C     harmonic viscosity=
18    C       viscAh (or viscAhD on div pts and viscAhZ on zeta pts)
19    C       +0.25*L**2*viscAhGrid/deltaT
20    C       +sqrt((viscC2leith/pi)**6*grad(Vort3)**2
21    C             +(viscC2leithD/pi)**6*grad(hDiv)**2)*L**3
22    C       +(viscC2smag/pi)**2*L**2*sqrt(Tension**2+Strain**2)
23    C
24    C     biharmonic viscosity=
25    C       viscA4 (or viscA4D on div pts and viscA4Z on zeta pts)
26    C       +0.25*0.125*L**4*viscA4Grid/deltaT (approx)
27    C       +0.125*L**5*sqrt((viscC4leith/pi)**6*grad(Vort3)**2
28    C                        +(viscC4leithD/pi)**6*grad(hDiv)**2)
29    C       +0.125*L**4*(viscC4smag/pi)**2*sqrt(Tension**2+Strain**2)
30    C
31    C     Note that often 0.125*L**2 is the scale between harmonic and
32    C     biharmonic (see Griffies and Hallberg (2000))
33    C     This allows the same value of the coefficient to be used
34    C     for roughly similar results with biharmonic and harmonic
35    C
36    C     LIMITERS -- limit min and max values of viscosities
37    C     viscAhReMax is min value for grid point harmonic Reynolds num
38    C      harmonic viscosity>sqrt(2*KE)*L/viscAhReMax
39    C
40    C     viscA4ReMax is min value for grid point biharmonic Reynolds num
41    C      biharmonic viscosity>sqrt(2*KE)*L**3/8/viscA4ReMax
42    C
43    C     viscAhgridmax is CFL stability limiter for harmonic viscosity
44    C      harmonic viscosity<0.25*viscAhgridmax*L**2/deltaT
45    C
46    C     viscA4gridmax is CFL stability limiter for biharmonic viscosity
47    C      biharmonic viscosity<viscA4gridmax*L**4/32/deltaT (approx)
48    C
49    C     viscAhgridmin and viscA4gridmin are lower limits for viscosity:
50    C       harmonic viscosity>0.25*viscAhgridmin*L**2/deltaT
51    C       biharmonic viscosity>viscA4gridmin*L**4/32/deltaT (approx)
52    
53    
54    C
55    C     RECOMMENDED VALUES
56    C     viscC2Leith=1-3
57    C     viscC2LeithD=1-3
58    C     viscC4Leith=1-3
59    C     viscC4LeithD=1.5-3
60    C     viscC2smag=2.2-4 (Griffies and Hallberg,2000)
61    C               0.2-0.9 (Smagorinsky,1993)
62    C     viscC4smag=2.2-4 (Griffies and Hallberg,2000)
63    C     viscAhReMax>=1, (<2 suppresses a computational mode)
64    C     viscA4ReMax>=1, (<2 suppresses a computational mode)
65    C     viscAhgridmax=1
66    C     viscA4gridmax=1
67    C     viscAhgrid<1
68    C     viscA4grid<1
69    C     viscAhgridmin<<1
70    C     viscA4gridmin<<1
71    
72  C     == Global variables ==  C     == Global variables ==
73  #include "SIZE.h"  #include "SIZE.h"
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 35  C     == Routine arguments == Line 99  C     == Routine arguments ==
99    
100  C     == Local variables ==  C     == Local variables ==
101        INTEGER I,J        INTEGER I,J
102        _RL ASmag2, ASmag4, smag2fac, smag4fac  #ifdef ALLOW_NONHYDROSTATIC
103        _RL vg2Min, vg2Max, AlinMax, AlinMin        INTEGER kp1
104        _RL lenA2, lenAz2  #endif
105        _RL Alin,Alth2,Alth4,grdVrt,grdDiv        INTEGER lockey
106        _RL vg2,vg4,vg4Min,vg4Max        _RL smag2fac, smag4fac
107          _RL leith2fac, leith4fac
108          _RL leithD2fac, leithD4fac
109          _RL viscAhRe_max, viscA4Re_max
110          _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
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)
118          _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119          _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120          _RL viscA4_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121          _RL viscAh_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122          _RL viscAh_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123          _RL viscA4_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124          _RL viscA4_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125          _RL viscAh_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126          _RL viscAh_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127          _RL viscA4_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128          _RL viscA4_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129          _RL viscAh_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130          _RL viscAh_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131          _RL viscA4_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132          _RL viscA4_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133          _RL viscAh_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134          _RL viscAh_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135          _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136          _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137          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 71  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 (deltaTmom.NE.0.) THEN        IF (useVariableViscosity) THEN
183         recip_dt=1./deltaTmom  C---- variable viscosity :
       ELSE  
        recip_dt=0.  
       ENDIF  
184    
185        vg2=viscAhGrid*recip_dt         IF ((harmonic).AND.(viscAhReMax.NE.0.)) THEN
186        vg2Min=viscAhGridMin*recip_dt          viscAhRe_max=SQRT(2. _d 0)/viscAhReMax
187        vg2Max=viscAhGridMax*recip_dt         ELSE
188        vg4=viscA4Grid*recip_dt          viscAhRe_max=0. _d 0
189        vg4Min=viscA4GridMin*recip_dt         ENDIF
190        vg4Max=viscA4GridMax*recip_dt  
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        smag2fac=(viscC2smag/pi)**2         calcLeith=
198        smag4fac=0.125*(viscC4smag/pi)**2       &      (viscC2leith.NE.0.)
199         &  .OR.(viscC2leithD.NE.0.)
200         &  .OR.(viscC4leith.NE.0.)
201         &  .OR.(viscC4leithD.NE.0.)
202    
203           calcSmag=
204         &      (viscC2smag.NE.0.)
205         &  .OR.(viscC4smag.NE.0.)
206    
207           IF (deltaTmom.NE.0.) THEN
208            recip_dt=1. _d 0/deltaTmom
209           ELSE
210            recip_dt=0. _d 0
211           ENDIF
212    
213           IF (calcSmag) THEN
214            smag2fac=(viscC2smag/pi)**2
215            smag4fac=0.125 _d 0*(viscC4smag/pi)**2
216           ELSE
217            smag2fac=0. _d 0
218            smag4fac=0. _d 0
219           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  C These are (powers of) length scales  
342           L2=rA(i,j,bi,bj)  C These are (powers of) length scales
343             IF (useAreaViscLength) THEN
344              L2=rA(i,j,bi,bj)
345              L4rdt=0.03125 _d 0*recip_dt*L2**2
346             ELSE
347              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
353           L3=(L2**1.5)           L3=(L2**1.5)
354           L4=(L2**2)           L4=(L2**2)
355           L5=0.125*(L2**2.5)           L5=(L2*L3)
356           IF (useAnisotropicViscAGridMax) THEN  
357            L2rdt=recip_dt/( 2.*(recip_DXF(I,J,bi,bj)**2           L2rdt=0.25 _d 0*recip_dt*L2
358       &                       +recip_DYF(I,J,bi,bj)**2) )  
359            L4rdt=recip_dt/( 6.*(recip_DXF(I,J,bi,bj)**4  C Velocity Reynolds Scale
360       &                       +recip_DYF(I,J,bi,bj)**4)           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
361       &                   +8.*((recip_DXF(I,J,bi,bj)             Uscl=SQRT(KE(i,j)*L2)*viscAhRe_max
362       &                        *recip_DYF(I,J,bi,bj))**2) )           ELSE
363               Uscl=0.
364             ENDIF
365             IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
366               U4scl=SQRT(KE(i,j))*L3*viscA4Re_max
367             ELSE
368               U4scl=0.
369           ENDIF           ENDIF
370    
371           IF (useFullLeith) 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*(            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))*recip_DXG(i,j+1,bi,bj))**2       &                        + vrtDy(i,j)*vrtDy(i,j) )  )
      &     +((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*(            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) )  )
386       &      +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2)  
387              viscAh_DLth(i,j)=
388             IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)       &     SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
389       &          .NE. 0. ) THEN            viscA4_DLth(i,j)=
390              Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
391             ELSE            viscAh_DLthd(i,j)=
392              Alth2=0. _d 0       &     SQRT(leithD2fac*grdDiv)*L3
393             ENDIF            viscA4_DLthd(i,j)=
394             IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)       &     SQRT(leithD4fac*grdDiv)*L5
395       &          .NE. 0. ) THEN           ELSEIF (calcLeith) THEN
             Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5  
            ELSE  
             Alth4=0. _d 0  
            ENDIF  
          ELSE  
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             Alth2=(viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3            grdDiv=MAX( ABS(divDx(i+1,j)), ABS(divDx(i,j)) )
404             Alth4=(viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5            grdDiv=MAX( grdDiv, ABS(divDy(i,j+1)) )
405           ENDIF            grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
406    
407           IF (smag2fac.NE.0.) THEN            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
408            Asmag2=smag2fac*L2            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
409       &           *sqrt(tension(i,j)**2            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
410       &           +0.25*(strain(i+1, j )**2+strain( i ,j+1)**2            viscA4_DlthD(i,j)=((leithD4fac*grdDiv))*L5
      &                 +strain(i-1, j )**2+strain( i ,j-1)**2))  
411           ELSE           ELSE
412            Asmag2=0d0            viscAh_Dlth(i,j)=0. _d 0
413              viscA4_Dlth(i,j)=0. _d 0
414              viscAh_DlthD(i,j)=0. _d 0
415              viscA4_DlthD(i,j)=0. _d 0
416           ENDIF           ENDIF
417    
418           IF (smag4fac.NE.0.) THEN           IF (calcSmag) THEN
419            Asmag4=smag4fac*L4            viscAh_DSmg(i,j)=L2
420       &             *sqrt(tension(i,j)**2       &       *SQRT(tension(i,j)**2
421       &             +0.25*(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-1, j )**2+strain( i ,j-1)**2))       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))
423              viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
424              viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)
425           ELSE           ELSE
426            Asmag4=0d0            viscAh_DSmg(i,j)=0. _d 0
427              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+vg2*L2+Alth2+Asmag2           Alin=viscAhD+viscAhGrid*L2rdt
433           viscAh_D(i,j)=min(viscAhMax,Alin)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
434           IF (useAnisotropicViscAGridMax) THEN           viscAh_DMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
435            AlinMax=viscAhGridMax*L2rdt           viscAh_D(i,j)=MAX(viscAh_DMin(i,j),Alin)
436            viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))           viscAh_DMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
437           ELSE           viscAh_D(i,j)=MIN(viscAh_DMax(i,j),viscAh_D(i,j))
           IF (vg2Max.GT.0.) THEN  
            AlinMax=vg2Max*L2  
            viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))  
           ENDIF  
          ENDIF  
          AlinMin=vg2Min*L2  
          viscAh_D(i,j)=max(AlinMin,viscAh_D(i,j))  
438    
439  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
440           Alin=viscA4D+vg4*L4+Alth4+Asmag4           Alin=viscA4D+viscA4Grid*L4rdt
441           viscA4_D(i,j)=min(viscA4Max,Alin)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
442           IF (useAnisotropicViscAGridMax) THEN           viscA4_DMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
443              AlinMax=viscA4GridMax*L4rdt           viscA4_D(i,j)=MAX(viscA4_DMin(i,j),Alin)
444              viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))           viscA4_DMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
445             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           ELSE
466           IF (vg4Max.GT.0.) THEN  C     Prepare for next level (next call)
467              AlinMax=vg4Max*L4            viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
468              viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))            viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
469           ENDIF  
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           ENDIF
475           AlinMin=vg4Min*L4  #endif /* ALLOW_NONHYDROSTATIC */
          viscA4_D(i,j)=max(AlinMin,viscA4_D(i,j))  
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           L2=rAz(i,j,bi,bj)           IF (useAreaViscLength) THEN
480              L2=rAz(i,j,bi,bj)
481              L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2
482             ELSE
483              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
488    
489           L3=(L2**1.5)           L3=(L2**1.5)
490           L4=(L2**2)           L4=(L2**2)
491           L5=0.125*(L2**2.5)           L5=(L2*L3)
          IF (useAnisotropicViscAGridMax) THEN  
          L2rdt=recip_dt/( 2.*(recip_DXV(I,J,bi,bj)**2  
      &                       +recip_DYU(I,J,bi,bj)**2) )  
          L4rdt=recip_dt/( 6.*(recip_DXV(I,J,bi,bj)**4  
      &                       +recip_DYU(I,J,bi,bj)**4)  
      &                   +8.*((recip_DXV(I,J,bi,bj)  
      &                        *recip_DYU(I,J,bi,bj))**2) )  
          ENDIF  
492    
493  C This is the vector magnitude of the vorticity gradient squared           L2rdt=0.25 _d 0*recip_dt*L2
          IF (useFullLeith) THEN  
            grdVrt=0.25*(  
      &      ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2  
      &      +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2  
      &      +((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj))**2  
      &      +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)  
494    
495  C This is the vector magnitude of grad(div.v) squared  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
496             grdDiv=0.25*(           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
497       &       ((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2             keZpt=0.25 _d 0*( (KE(i,j)+KE(i-1,j-1))
498       &      +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2       &                      +(KE(i-1,j)+KE(i,j-1)) )
499       &      +((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj))**2             IF ( keZpt.GT.0. ) THEN
500       &      +((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYC(i-1,j,bi,bj))**2)               Uscl = SQRT(keZpt*L2)*viscAhRe_max
501                 U4scl= SQRT(keZpt)*L3*viscA4Re_max
            IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)  
      &          .NE. 0. ) THEN  
             Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3  
            ELSE  
             Alth2=0. _d 0  
            ENDIF  
            IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)  
      &          .NE. 0. ) THEN  
             Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5  
502             ELSE             ELSE
503              Alth4=0. _d 0               Uscl =0.
504                 U4scl=0.
505             ENDIF             ENDIF
506           ELSE           ELSE
507               Uscl =0.
508               U4scl=0.
509             ENDIF
510    
511    #ifndef AUTODIFF_DISABLE_LEITH
512    C This is the vector magnitude of the vorticity gradient squared
513             IF (useFullLeith.AND.calcLeith) THEN
514              grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
515         &                        + vrtDx(i,j)*vrtDx(i,j) )
516         &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
517         &                        + vrtDy(i,j)*vrtDy(i,j) )  )
518    
519    C This is the vector magnitude of grad(div.v) squared
520              grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
521         &                        + divDx(i,j)*divDx(i,j) )
522         &                     + (divDy(i-1,j)*divDy(i-1,j)
523         &                        + divDy(i,j)*divDy(i,j) )  )
524    
525              viscAh_ZLth(i,j)=
526         &     SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
527              viscA4_ZLth(i,j)=
528         &     SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
529              viscAh_ZLthD(i,j)=
530         &     SQRT(leithD2fac*grdDiv)*L3
531              viscA4_ZLthD(i,j)=
532         &     SQRT(leithD4fac*grdDiv)*L5
533    
534             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_DXG(i,j-1,bi,bj)))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYG(i-1,j,bi,bj)))  
   
   
 C This if statement is just to prevent bitwise changes when leithd=0  
            Alth2=(viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3  
            Alth4=(viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5  
          ENDIF  
   
          IF (smag2fac.NE.0.) THEN  
           Asmag2=smag2fac*L2  
      &           *sqrt(strain(i,j)**2  
      &              +0.25*(tension( i , j )**2+tension( i ,j+1)**2  
      &                +tension(i+1, j )**2+tension(i+1,j+1)**2))  
548           ELSE           ELSE
549            Asmag2=0d0            viscAh_ZLth(i,j)=0. _d 0
550              viscA4_ZLth(i,j)=0. _d 0
551              viscAh_ZLthD(i,j)=0. _d 0
552              viscA4_ZLthD(i,j)=0. _d 0
553           ENDIF           ENDIF
554    
555           IF (smag4fac.NE.0.) THEN           IF (calcSmag) THEN
556            Asmag4=smag4fac*L4            viscAh_ZSmg(i,j)=L2
557       &             *sqrt(strain(i,j)**2       &      *SQRT(strain(i,j)**2
558       &             +0.25*(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           ELSE            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
561            Asmag4=0d0            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+vg2*L2+Alth2+Asmag2           Alin=viscAhZ+viscAhGrid*L2rdt
567           viscAh_Z(i,j)=min(viscAhMax,Alin)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
568           IF (useAnisotropicViscAGridMax) THEN           viscAh_ZMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
569              AlinMax=viscAhGridMax*L2rdt           viscAh_Z(i,j)=MAX(viscAh_ZMin(i,j),Alin)
570              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           viscAh_ZMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
571           ELSE           viscAh_Z(i,j)=MIN(viscAh_ZMax(i,j),viscAh_Z(i,j))
572           IF (vg2Max.GT.0.) THEN  
573              AlinMax=vg2Max*L2  C  BiHarmonic on Zeta points
574              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           Alin=viscA4Z+viscA4Grid*L4rdt
575           ENDIF       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
576           ENDIF           viscA4_ZMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
577           AlinMin=vg2Min*L2           viscA4_Z(i,j)=MAX(viscA4_ZMin(i,j),Alin)
578           viscAh_Z(i,j)=max(AlinMin,viscAh_Z(i,j))           viscA4_ZMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
579             viscA4_Z(i,j)=MIN(viscA4_ZMax(i,j),viscA4_Z(i,j))
 C  BiHarmonic on Zeta Points  
          Alin=viscA4Z+vg4*L4+Alth4+Asmag4  
          viscA4_Z(i,j)=min(viscA4Max,Alin)  
          IF (useAnisotropicViscAGridMax) THEN  
             AlinMax=viscA4GridMax*L4rdt  
             viscA4_Z(i,j)=min(AlinMax,viscA4_Z(i,j))  
          ELSE  
          IF (vg4Max.GT.0.) THEN  
             AlinMax=vg4Max*L4  
             viscA4_Z(i,j)=min(AlinMax,viscA4_Z(i,j))  
          ENDIF  
          ENDIF  
          AlinMin=vg4Min*L4  
          viscA4_Z(i,j)=max(AlinMin,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 330  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
# Line 338  C  BiHarmonic on Zeta Points Line 601  C  BiHarmonic on Zeta Points
601         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
602         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
603         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
604    
605           CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
606           CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)
607           CALL DIAGNOSTICS_FILL(viscAh_ZMax,'VAHZMAX ',k,1,2,bi,bj,myThid)
608           CALL DIAGNOSTICS_FILL(viscA4_ZMax,'VA4ZMAX ',k,1,2,bi,bj,myThid)
609    
610           CALL DIAGNOSTICS_FILL(viscAh_DMin,'VAHDMIN ',k,1,2,bi,bj,myThid)
611           CALL DIAGNOSTICS_FILL(viscA4_DMin,'VA4DMIN ',k,1,2,bi,bj,myThid)
612           CALL DIAGNOSTICS_FILL(viscAh_ZMin,'VAHZMIN ',k,1,2,bi,bj,myThid)
613           CALL DIAGNOSTICS_FILL(viscA4_ZMin,'VA4ZMIN ',k,1,2,bi,bj,myThid)
614    
615           CALL DIAGNOSTICS_FILL(viscAh_DLth,'VAHDLTH ',k,1,2,bi,bj,myThid)
616           CALL DIAGNOSTICS_FILL(viscA4_DLth,'VA4DLTH ',k,1,2,bi,bj,myThid)
617           CALL DIAGNOSTICS_FILL(viscAh_ZLth,'VAHZLTH ',k,1,2,bi,bj,myThid)
618           CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)
619    
620           CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD'
621         &   ,k,1,2,bi,bj,myThid)
622           CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD'
623         &   ,k,1,2,bi,bj,myThid)
624           CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD'
625         &   ,k,1,2,bi,bj,myThid)
626           CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD'
627         &   ,k,1,2,bi,bj,myThid)
628    
629           CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)
630           CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)
631           CALL DIAGNOSTICS_FILL(viscAh_ZSmg,'VAHZSMAG',k,1,2,bi,bj,myThid)
632           CALL DIAGNOSTICS_FILL(viscA4_ZSmg,'VA4ZSMAG',k,1,2,bi,bj,myThid)
633        ENDIF        ENDIF
634  #endif  #endif
635    
636        RETURN        RETURN
637        END        END
638    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22