/[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.23 by baylor, Fri Jul 7 18:52:10 2006 UTC revision 1.36 by jmc, Wed Oct 22 00:28:47 2008 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 75  C     == Global variables == Line 77  C     == Global variables ==
77  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
78  #include "NH_VARS.h"  #include "NH_VARS.h"
79  #endif  #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 93  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        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 125  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 146  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 161  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 203  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  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
241        IF ( calcLeith .OR. calcSmag ) THEN  cphtest       IF ( calcLeith .OR. calcSmag ) THEN
242         STOP 'calcLeith or calcSmag not implemented for ADJOINT'  cphtest        STOP 'calcLeith or calcSmag not implemented for ADJOINT'
243        ENDIF  cphtest       ENDIF
244        DO j=1-Oly,sNy+Oly  #endif
245           DO j=1-Oly,sNy+Oly
246          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
247            viscAh_D(i,j)=viscAhD            viscAh_D(i,j)=viscAhD
248            viscAh_Z(i,j)=viscAhZ            viscAh_Z(i,j)=viscAhZ
# Line 237  c Line 265  c
265            viscAh_ZLthD(i,j)= 0. _d 0            viscAh_ZLthD(i,j)= 0. _d 0
266            viscA4_ZLthD(i,j)= 0. _d 0            viscA4_ZLthD(i,j)= 0. _d 0
267          ENDDO          ENDDO
268        ENDDO         ENDDO
 #endif  
   
   
   
 C     - Viscosity  
       IF (useVariableViscosity) THEN  
269    
270  C-     Initialise to zero gradient of vorticity & divergence:  C-     Initialise to zero gradient of vorticity & divergence:
271         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
# Line 255  C-     Initialise to zero gradient of vo Line 277  C-     Initialise to zero gradient of vo
277          ENDDO          ENDDO
278         ENDDO         ENDDO
279    
280         IF (calcleith) THEN         IF (calcLeith) THEN
281  C      horizontal gradient of horizontal divergence:  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 272  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)
# Line 308  C-       gradient in y direction: Line 332  C-       gradient in y direction:
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 + sNx*(j-1) + sNx*sNy*(lockey_1-1)
338    CADJ STORE viscA4_ZSmg(i,j)
339    CADJ &     = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte
340    CADJ STORE viscAh_ZSmg(i,j)
341    CADJ &     = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte
342    # endif
343    #endif    /* ALLOW_AUTODIFF_TAMC */
344    
345    C These are (powers of) length scales
346           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
347            L2=rA(i,j,bi,bj)            L2=rA(i,j,bi,bj)
348              L4rdt=0.03125 _d 0*recip_dt*L2**2
349           ELSE           ELSE
350            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))
351              L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4
352         &                             +recip_DYF(I,J,bi,bj)**4)
353         &                   +8. _d 0*((recip_DXF(I,J,bi,bj)
354         &                             *recip_DYF(I,J,bi,bj))**2) )
355           ENDIF           ENDIF
356           L3=(L2**1.5)           L3=(L2**1.5)
357           L4=(L2**2)           L4=(L2**2)
358           L5=(L2**2.5)           L5=(L2*L3)
359    
360           L2rdt=0.25 _d 0*recip_dt*L2           L2rdt=0.25 _d 0*recip_dt*L2
361    
          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  
   
362  C Velocity Reynolds Scale  C Velocity Reynolds Scale
363           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
364             Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max             Uscl=SQRT(KE(i,j)*L2)*viscAhRe_max
365           ELSE           ELSE
366             Uscl=0.             Uscl=0.
367           ENDIF           ENDIF
368           IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN           IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
369             U4scl=sqrt(KE(i,j))*L3*viscA4Re_max             U4scl=SQRT(KE(i,j))*L3*viscA4Re_max
370           ELSE           ELSE
371             U4scl=0.             U4scl=0.
372           ENDIF           ENDIF
373    
374  #ifndef ALLOW_AUTODIFF_TAMC  cph-leith#ifndef ALLOW_AUTODIFF_TAMC
375           IF (useFullLeith.and.calcleith) THEN  #ifndef AUTODIFF_DISABLE_LEITH
376             IF (useFullLeith.AND.calcLeith) THEN
377  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
378            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
379       &                        + vrtDx(i,j)*vrtDx(i,j) )       &                        + vrtDx(i,j)*vrtDx(i,j) )
# Line 357  C Using it in Leith serves to damp insta Line 388  C Using it in Leith serves to damp insta
388       &                        + divDy(i,j)*divDy(i,j) )  )       &                        + divDy(i,j)*divDy(i,j) )  )
389    
390            viscAh_DLth(i,j)=            viscAh_DLth(i,j)=
391       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3       &     SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
392            viscA4_DLth(i,j)=            viscA4_DLth(i,j)=
393       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5       &     SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
394            viscAh_DLthd(i,j)=            viscAh_DLthd(i,j)=
395       &     sqrt(leithD2fac*grdDiv)*L3       &     SQRT(leithD2fac*grdDiv)*L3
396            viscA4_DLthd(i,j)=            viscA4_DLthd(i,j)=
397       &     sqrt(leithD4fac*grdDiv)*L5       &     SQRT(leithD4fac*grdDiv)*L5
398           ELSEIF (calcleith) THEN           ELSEIF (calcLeith) THEN
399  C but this approximation will work on cube  C but this approximation will work on cube
400  c (and differs by as much as 4X)  c (and differs by as much as 4X)
401            grdVrt=max( abs(vrtDx(i,j+1)), abs(vrtDx(i,j)) )            grdVrt=MAX( ABS(vrtDx(i,j+1)), ABS(vrtDx(i,j)) )
402            grdVrt=max( grdVrt, abs(vrtDy(i+1,j)) )            grdVrt=MAX( grdVrt, ABS(vrtDy(i+1,j)) )
403            grdVrt=max( grdVrt, abs(vrtDy(i,j))   )            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j))   )
404    
405  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
406            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )            grdDiv=MAX( ABS(divDx(i+1,j)), ABS(divDx(i,j)) )
407            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )            grdDiv=MAX( grdDiv, ABS(divDy(i,j+1)) )
408            grdDiv=max( grdDiv, abs(divDy(i,j))   )            grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
409    
410            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
411            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
# Line 387  c This approximation is good to the same Line 418  c This approximation is good to the same
418            viscA4_DlthD(i,j)=0. _d 0            viscA4_DlthD(i,j)=0. _d 0
419           ENDIF           ENDIF
420    
421           IF (calcsmag) THEN           IF (calcSmag) THEN
422            viscAh_DSmg(i,j)=L2            viscAh_DSmg(i,j)=L2
423       &       *sqrt(tension(i,j)**2       &       *SQRT(tension(i,j)**2
424       &       +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
425       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))
426            viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)            viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
# Line 398  c This approximation is good to the same Line 429  c This approximation is good to the same
429            viscAh_DSmg(i,j)=0. _d 0            viscAh_DSmg(i,j)=0. _d 0
430            viscA4_DSmg(i,j)=0. _d 0            viscA4_DSmg(i,j)=0. _d 0
431           ENDIF           ENDIF
432  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* AUTODIFF_DISABLE_LEITH */
433    
434  C  Harmonic on Div.u points  C  Harmonic on Div.u points
435           Alin=viscAhD+viscAhGrid*L2rdt           Alin=viscAhD+viscAhGrid*L2rdt
436       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
437           viscAh_DMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)           viscAh_DMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
438           viscAh_D(i,j)=max(viscAh_DMin(i,j),Alin)           viscAh_D(i,j)=MAX(viscAh_DMin(i,j),Alin)
439           viscAh_DMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)           viscAh_DMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
440           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))
441    
442  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
443           Alin=viscA4D+viscA4Grid*L4rdt           Alin=viscA4D+viscA4Grid*L4rdt
444       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
445           viscA4_DMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)           viscA4_DMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
446           viscA4_D(i,j)=max(viscA4_DMin(i,j),Alin)           viscA4_D(i,j)=MAX(viscA4_DMin(i,j),Alin)
447           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)           viscA4_DMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
448           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))
449    
450  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
451  C     /* Pass Viscosities to calc_gw, if constant, not necessary */  C--   Pass Viscosities to calc_gw, if constant, not necessary
452    
453           kp1  = MIN(k+1,Nr)           kp1  = MIN(k+1,Nr)
454    
455           if (k .eq. 1) then           IF ( k.EQ.1 ) THEN
456    C     Prepare for next level (next call)
457            viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)            viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
458            viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)            viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
459    
460            viscAh_W(i,j,k,bi,bj)=viscAh_D(i,j) /* These values dont get used */  C     These values dont get used
461              viscAh_W(i,j,k,bi,bj)=viscAh_D(i,j)
462            viscA4_W(i,j,k,bi,bj)=viscA4_D(i,j)            viscA4_W(i,j,k,bi,bj)=viscA4_D(i,j)
463           else  
464  C Note that previous call of this function has already added half.           ELSEIF ( k.EQ.Nr ) THEN
465              viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)
466              viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)
467    
468             ELSE
469    C     Prepare for next level (next call)
470            viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)            viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
471            viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)            viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
472    
473    C     Note that previous call of this function has already added half.
474            viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)            viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)
475            viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)            viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)
476           endif  
477             ENDIF
478  #endif /* ALLOW_NONHYDROSTATIC */  #endif /* ALLOW_NONHYDROSTATIC */
479    
480  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
481  C These are (powers of) length scales  C These are (powers of) length scales
482           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
483            L2=rAz(i,j,bi,bj)            L2=rAz(i,j,bi,bj)
484              L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2
485           ELSE           ELSE
486            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))
487              L4rdt=recip_dt/
488         &     ( 6. _d 0*(recip_DXV(I,J,bi,bj)**4+recip_DYU(I,J,bi,bj)**4)
489         &      +8. _d 0*((recip_DXV(I,J,bi,bj)*recip_DYU(I,J,bi,bj))**2))
490           ENDIF           ENDIF
491    
492           L3=(L2**1.5)           L3=(L2**1.5)
493           L4=(L2**2)           L4=(L2**2)
494           L5=(L2**2.5)           L5=(L2*L3)
495    
496           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  
497    
498  C Velocity Reynolds Scale (Pb here at CS-grid corners !)  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
499           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
500             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))
501       &                      +(KE(i-1,j)+KE(i,j-1)) )       &                      +(KE(i-1,j)+KE(i,j-1)) )
502             IF ( keZpt.GT.0. ) THEN             IF ( keZpt.GT.0. ) THEN
503               Uscl = sqrt(keZpt*L2)*viscAhRe_max               Uscl = SQRT(keZpt*L2)*viscAhRe_max
504               U4scl= sqrt(keZpt)*L3*viscA4Re_max               U4scl= SQRT(keZpt)*L3*viscA4Re_max
505             ELSE             ELSE
506               Uscl =0.               Uscl =0.
507               U4scl=0.               U4scl=0.
# Line 474  C Velocity Reynolds Scale (Pb here at CS Line 511  C Velocity Reynolds Scale (Pb here at CS
511             U4scl=0.             U4scl=0.
512           ENDIF           ENDIF
513    
514  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef AUTODIFF_DISABLE_LEITH
515  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
516           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.AND.calcLeith) THEN
517            grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)            grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
518       &                        + vrtDx(i,j)*vrtDx(i,j) )       &                        + vrtDx(i,j)*vrtDx(i,j) )
519       &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)       &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
# Line 489  C This is the vector magnitude of grad(d Line 526  C This is the vector magnitude of grad(d
526       &                        + divDy(i,j)*divDy(i,j) )  )       &                        + divDy(i,j)*divDy(i,j) )  )
527    
528            viscAh_ZLth(i,j)=            viscAh_ZLth(i,j)=
529       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3       &     SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
530            viscA4_ZLth(i,j)=            viscA4_ZLth(i,j)=
531       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5       &     SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
532            viscAh_ZLthD(i,j)=            viscAh_ZLthD(i,j)=
533       &     sqrt(leithD2fac*grdDiv)*L3       &     SQRT(leithD2fac*grdDiv)*L3
534            viscA4_ZLthD(i,j)=            viscA4_ZLthD(i,j)=
535       &     sqrt(leithD4fac*grdDiv)*L5       &     SQRT(leithD4fac*grdDiv)*L5
536    
537           ELSEIF (calcleith) THEN           ELSEIF (calcLeith) THEN
538  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
539            grdVrt=max( abs(vrtDx(i-1,j)), abs(vrtDx(i,j)) )            grdVrt=MAX( ABS(vrtDx(i-1,j)), ABS(vrtDx(i,j)) )
540            grdVrt=max( grdVrt, abs(vrtDy(i,j-1)) )            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j-1)) )
541            grdVrt=max( grdVrt, abs(vrtDy(i,j))   )            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j))   )
542    
543            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )            grdDiv=MAX( ABS(divDx(i,j)), ABS(divDx(i,j-1)) )
544            grdDiv=max( grdDiv, abs(divDy(i,j))   )            grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
545            grdDiv=max( grdDiv, abs(divDy(i-1,j)) )            grdDiv=MAX( grdDiv, ABS(divDy(i-1,j)) )
546    
547            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
548            viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5            viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
# Line 518  C but this approximation will work on cu Line 555  C but this approximation will work on cu
555            viscA4_ZLthD(i,j)=0. _d 0            viscA4_ZLthD(i,j)=0. _d 0
556           ENDIF           ENDIF
557    
558           IF (calcsmag) THEN           IF (calcSmag) THEN
559            viscAh_ZSmg(i,j)=L2            viscAh_ZSmg(i,j)=L2
560       &      *sqrt(strain(i,j)**2       &      *SQRT(strain(i,j)**2
561       &        +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
562       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))
563            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
564            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
565           ENDIF           ENDIF
566  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* AUTODIFF_DISABLE_LEITH */
567    
568  C  Harmonic on Zeta points  C  Harmonic on Zeta points
569           Alin=viscAhZ+viscAhGrid*L2rdt           Alin=viscAhZ+viscAhGrid*L2rdt
570       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
571           viscAh_ZMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)           viscAh_ZMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
572           viscAh_Z(i,j)=max(viscAh_ZMin(i,j),Alin)           viscAh_Z(i,j)=MAX(viscAh_ZMin(i,j),Alin)
573           viscAh_ZMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)           viscAh_ZMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
574           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))
575    
576  C  BiHarmonic on Zeta points  C  BiHarmonic on Zeta points
577           Alin=viscA4Z+viscA4Grid*L4rdt           Alin=viscA4Z+viscA4Grid*L4rdt
578       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
579           viscA4_ZMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)           viscA4_ZMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
580           viscA4_Z(i,j)=max(viscA4_ZMin(i,j),Alin)           viscA4_Z(i,j)=MAX(viscA4_ZMin(i,j),Alin)
581           viscA4_ZMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)           viscA4_ZMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
582           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))
583          ENDDO          ENDDO
584         ENDDO         ENDDO
585    
586        ELSE        ELSE
587    C---- use constant viscosity (useVariableViscosity=F):
588    
589         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
590          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
591           viscAh_D(i,j)=viscAhD           viscAh_D(i,j)=viscAhD
# Line 554  C  BiHarmonic on Zeta points Line 594  C  BiHarmonic on Zeta points
594           viscA4_Z(i,j)=viscA4Z           viscA4_Z(i,j)=viscA4Z
595          ENDDO          ENDDO
596         ENDDO         ENDDO
597    
598    C---- variable/constant viscosity : end if/else block
599        ENDIF        ENDIF
600    
601  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
# Line 562  C  BiHarmonic on Zeta points Line 604  C  BiHarmonic on Zeta points
604         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
605         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
606         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
 #ifdef ALLOW_NONHYDROSTATIC  
        CALL DIAGNOSTICS_FILL(viscAh_W,'VISCAHW ',k,1,2,bi,bj,myThid)  
        CALL DIAGNOSTICS_FILL(viscA4_W,'VISCA4W ',k,1,2,bi,bj,myThid)  
 #endif  
607    
608         CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
609         CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.36
  ViewVC Help
Powered by ViewVC 1.1.22