/[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.37 by heimbach, Mon Apr 6 23:47:06 2009 UTC revision 1.48 by gforget, Sat Jun 28 22:36:04 2014 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "MOM_COMMON_OPTIONS.h"  #include "MOM_COMMON_OPTIONS.h"
5    #ifdef ALLOW_AUTODIFF
6    # include "AUTODIFF_OPTIONS.h"
7    #endif
8    
9    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
10    CBOP
11    C     !ROUTINE: MOM_CALC_VISC
12    
13    C     !INTERFACE:
14        SUBROUTINE MOM_CALC_VISC(        SUBROUTINE MOM_CALC_VISC(
15       I        bi,bj,k,       I        bi,bj,k,
16       O        viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,       O        viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
      O        harmonic,biharmonic,useVariableViscosity,  
17       I        hDiv,vort3,tension,strain,KE,hFacZ,       I        hDiv,vort3,tension,strain,KE,hFacZ,
18       I        myThid)       I        myThid)
19    
20        IMPLICIT NONE  C     !DESCRIPTION:
 C  
21  C     Calculate horizontal viscosities (L is typical grid width)  C     Calculate horizontal viscosities (L is typical grid width)
22  C     harmonic viscosity=  C     harmonic viscosity=
23  C       viscAh (or viscAhD on div pts and viscAhZ on zeta pts)  C       viscAh (or viscAhD on div pts and viscAhZ on zeta pts)
# Line 50  C     viscAhgridmin and viscA4gridmin ar Line 55  C     viscAhgridmin and viscA4gridmin ar
55  C       harmonic viscosity>0.25*viscAhgridmin*L**2/deltaT  C       harmonic viscosity>0.25*viscAhgridmin*L**2/deltaT
56  C       biharmonic viscosity>viscA4gridmin*L**4/32/deltaT (approx)  C       biharmonic viscosity>viscA4gridmin*L**4/32/deltaT (approx)
57    
   
 C  
58  C     RECOMMENDED VALUES  C     RECOMMENDED VALUES
59  C     viscC2Leith=1-3  C     viscC2Leith=1-3
60  C     viscC2LeithD=1-3  C     viscC2LeithD=1-3
# Line 69  C     viscA4grid<1 Line 72  C     viscA4grid<1
72  C     viscAhgridmin<<1  C     viscAhgridmin<<1
73  C     viscA4gridmin<<1  C     viscA4gridmin<<1
74    
75    C     !USES:
76          IMPLICIT NONE
77    
78  C     == Global variables ==  C     == Global variables ==
79  #include "SIZE.h"  #include "SIZE.h"
80  #include "GRID.h"  #include "GRID.h"
81  #include "EEPARAMS.h"  #include "EEPARAMS.h"
82  #include "PARAMS.h"  #include "PARAMS.h"
83  #ifdef ALLOW_NONHYDROSTATIC  #include "MOM_VISC.h"
84  #include "NH_VARS.h"  #ifdef ALLOW_AUTODIFF
 #endif  
 #ifdef ALLOW_AUTODIFF_TAMC  
85  #include "tamc.h"  #include "tamc.h"
86  #include "tamc_keys.h"  #include "tamc_keys.h"
87  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF */
88    
89  C     == Routine arguments ==  C     !INPUT/OUTPUT PARAMETERS:
90    C     myThid               :: my thread Id number
91        INTEGER bi,bj,k        INTEGER bi,bj,k
92        _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93        _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 95  C     == Routine arguments == Line 100  C     == Routine arguments ==
100        _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        INTEGER myThid        INTEGER myThid
103        LOGICAL harmonic,biharmonic,useVariableViscosity  CEOP
104    
105  C     == Local variables ==  C     !LOCAL VARIABLES:
106        INTEGER I,J        INTEGER i,j
107  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
108        INTEGER kp1        _RL shiftAh, shiftA4
109  #endif  #endif
110    #ifdef ALLOW_AUTODIFF_TAMC
111        INTEGER lockey_1, lockey_2        INTEGER lockey_1, lockey_2
112    #endif
113        _RL smag2fac, smag4fac        _RL smag2fac, smag4fac
114        _RL leith2fac, leith4fac        _RL leith2fac, leith4fac
115        _RL leithD2fac, leithD4fac        _RL leithD2fac, leithD4fac
116        _RL viscAhRe_max, viscA4Re_max        _RL viscAhRe_max, viscA4Re_max
117        _RL Alin,grdVrt,grdDiv, keZpt        _RL Alin,grdVrt,grdDiv, keZpt
118        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt        _RL L2, L3, L5, L2rdt, L4rdt, recip_dt
119        _RL Uscl,U4scl        _RL Uscl,U4scl
120        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 151  C     == Local variables == Line 158  C     == Local variables ==
158  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
159    
160  C--   Set flags which are used in this S/R and elsewhere :  C--   Set flags which are used in this S/R and elsewhere :
161        useVariableViscosity=  C     useVariableVisc, useHarmonicVisc and useBiharmonicVisc
162       &      (viscAhGrid.NE.0.)  C     are now set early on (in S/R SET_PARAMS)
      &  .OR.(viscA4Grid.NE.0.)  
      &  .OR.(viscC2leith.NE.0.)  
      &  .OR.(viscC2leithD.NE.0.)  
      &  .OR.(viscC4leith.NE.0.)  
      &  .OR.(viscC4leithD.NE.0.)  
      &  .OR.(viscC2smag.NE.0.)  
      &  .OR.(viscC4smag.NE.0.)  
   
       harmonic=  
      &      (viscAh.NE.0.)  
      &  .OR.(viscAhD.NE.0.)  
      &  .OR.(viscAhZ.NE.0.)  
      &  .OR.(viscAhGrid.NE.0.)  
      &  .OR.(viscC2leith.NE.0.)  
      &  .OR.(viscC2leithD.NE.0.)  
      &  .OR.(viscC2smag.NE.0.)  
   
       biharmonic=  
      &      (viscA4.NE.0.)  
      &  .OR.(viscA4D.NE.0.)  
      &  .OR.(viscA4Z.NE.0.)  
      &  .OR.(viscA4Grid.NE.0.)  
      &  .OR.(viscC4leith.NE.0.)  
      &  .OR.(viscC4leithD.NE.0.)  
      &  .OR.(viscC4smag.NE.0.)  
163    
164        IF (useVariableViscosity) THEN  c     IF ( useVariableVisc ) THEN
165  C---- variable viscosity :  C---- variable viscosity :
166    
167         IF ((harmonic).AND.(viscAhReMax.NE.0.)) THEN         recip_dt = 1. _d 0
168           IF ( deltaTMom.NE.0. ) recip_dt = 1. _d 0/deltaTMom
169    
170           IF ( useHarmonicVisc .AND. viscAhReMax.NE.0. ) THEN
171          viscAhRe_max=SQRT(2. _d 0)/viscAhReMax          viscAhRe_max=SQRT(2. _d 0)/viscAhReMax
172         ELSE         ELSE
173          viscAhRe_max=0. _d 0          viscAhRe_max=0. _d 0
174         ENDIF         ENDIF
175    
176         IF ((biharmonic).AND.(viscA4ReMax.NE.0.)) THEN         IF ( useBiharmonicVisc .AND. viscA4ReMax.NE.0. ) THEN
177          viscA4Re_max=0.125 _d 0*SQRT(2. _d 0)/viscA4ReMax          viscA4Re_max=0.125 _d 0*SQRT(2. _d 0)/viscA4ReMax
178         ELSE         ELSE
179          viscA4Re_max=0. _d 0          viscA4Re_max=0. _d 0
# Line 204  C---- variable viscosity : Line 189  C---- variable viscosity :
189       &      (viscC2smag.NE.0.)       &      (viscC2smag.NE.0.)
190       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
191    
        IF (deltaTmom.NE.0.) THEN  
         recip_dt=1. _d 0/deltaTmom  
        ELSE  
         recip_dt=0. _d 0  
        ENDIF  
   
192         IF (calcSmag) THEN         IF (calcSmag) THEN
193          smag2fac=(viscC2smag/pi)**2          smag2fac=(viscC2smag/pi)**2
194          smag4fac=0.125 _d 0*(viscC4smag/pi)**2          smag4fac=0.125 _d 0*(viscC4smag/pi)**2
# Line 237  C---- variable viscosity : Line 216  C---- variable viscosity :
216          leithD4fac=0. _d 0          leithD4fac=0. _d 0
217         ENDIF         ENDIF
218    
219  #ifdef ALLOW_AUTODIFF_TAMC         DO j=1-OLy,sNy+OLy
220  cphtest       IF ( calcLeith .OR. calcSmag ) THEN          DO i=1-OLx,sNx+OLx
221  cphtest        STOP 'calcLeith or calcSmag not implemented for ADJOINT'  C-    viscosity arrays have been initialised everywhere before calling this S/R
222  cphtest       ENDIF  c         viscAh_D(i,j) = viscAhD
223  #endif  c         viscAh_Z(i,j) = viscAhZ
224         DO j=1-Oly,sNy+Oly  c         viscA4_D(i,j) = viscA4D
225          DO i=1-Olx,sNx+Olx  c         viscA4_Z(i,j) = viscA4Z
226            viscAh_D(i,j)=viscAhD  
           viscAh_Z(i,j)=viscAhZ  
           viscA4_D(i,j)=viscA4D  
           viscA4_Z(i,j)=viscA4Z  
 c  
227            visca4_zsmg(i,j) = 0. _d 0            visca4_zsmg(i,j) = 0. _d 0
228            viscah_zsmg(i,j) = 0. _d 0            viscah_zsmg(i,j) = 0. _d 0
229  c  
230            viscAh_Dlth(i,j) = 0. _d 0            viscAh_Dlth(i,j) = 0. _d 0
231            viscA4_Dlth(i,j) = 0. _d 0            viscA4_Dlth(i,j) = 0. _d 0
232            viscAh_DlthD(i,j)= 0. _d 0            viscAh_DlthD(i,j)= 0. _d 0
233            viscA4_DlthD(i,j)= 0. _d 0            viscA4_DlthD(i,j)= 0. _d 0
234  c  
235            viscAh_DSmg(i,j) = 0. _d 0            viscAh_DSmg(i,j) = 0. _d 0
236            viscA4_DSmg(i,j) = 0. _d 0            viscA4_DSmg(i,j) = 0. _d 0
237  c  
238            viscAh_ZLth(i,j) = 0. _d 0            viscAh_ZLth(i,j) = 0. _d 0
239            viscA4_ZLth(i,j) = 0. _d 0            viscA4_ZLth(i,j) = 0. _d 0
240            viscAh_ZLthD(i,j)= 0. _d 0            viscAh_ZLthD(i,j)= 0. _d 0
# Line 267  c Line 242  c
242          ENDDO          ENDDO
243         ENDDO         ENDDO
244    
245  C-     Initialise to zero gradient of vorticity & divergence:  C-    Initialise to zero gradient of vorticity & divergence:
246         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
247          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
248            divDx(i,j) = 0.            divDx(i,j) = 0.
249            divDy(i,j) = 0.            divDy(i,j) = 0.
250            vrtDx(i,j) = 0.            vrtDx(i,j) = 0.
# Line 277  C-     Initialise to zero gradient of vo Line 252  C-     Initialise to zero gradient of vo
252          ENDDO          ENDDO
253         ENDDO         ENDDO
254    
255         IF (calcLeith) THEN         IF ( calcLeith ) THEN
256  C      horizontal gradient of horizontal divergence:  C--   horizontal gradient of horizontal divergence:
   
257  C-       gradient in x direction:  C-       gradient in x direction:
 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC  
258           IF (useCubedSphereExchange) THEN           IF (useCubedSphereExchange) THEN
259  C        to compute d/dx(hDiv), fill corners with appropriate values:  C        to compute d/dx(hDiv), fill corners with appropriate values:
260             CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,             CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
261       &                                hDiv, bi,bj, myThid )       &                                hDiv, bi,bj, myThid )
262           ENDIF           ENDIF
263  cph-exch2#endif           DO j=2-OLy,sNy+OLy-1
264           DO j=2-Oly,sNy+Oly-1            DO i=2-OLx,sNx+OLx-1
265            DO i=2-Olx,sNx+Olx-1              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)  
266            ENDDO            ENDDO
267           ENDDO           ENDDO
268    
269  C-       gradient in y direction:  C-       gradient in y direction:
 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC  
270           IF (useCubedSphereExchange) THEN           IF (useCubedSphereExchange) THEN
271  C        to compute d/dy(hDiv), fill corners with appropriate values:  C        to compute d/dy(hDiv), fill corners with appropriate values:
272             CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,             CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
273       &                                hDiv, bi,bj, myThid )       &                                hDiv, bi,bj, myThid )
274           ENDIF           ENDIF
275  cph-exch2#endif           DO j=2-OLy,sNy+OLy-1
276           DO j=2-Oly,sNy+Oly-1            DO i=2-OLx,sNx+OLx-1
277            DO i=2-Olx,sNx+Olx-1              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)  
278            ENDDO            ENDDO
279           ENDDO           ENDDO
280    
281  C      horizontal gradient of vertical vorticity:  C--   horizontal gradient of vertical vorticity:
282  C-       gradient in x direction:  C-       gradient in x direction:
283           DO j=2-Oly,sNy+Oly           DO j=2-OLy,sNy+OLy
284            DO i=2-Olx,sNx+Olx-1            DO i=2-OLx,sNx+OLx-1
285              vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))              vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
286       &                  *recip_DXG(i,j,bi,bj)       &                  *recip_dxG(i,j,bi,bj)
287       &                  *maskS(i,j,k,bi,bj)       &                  *maskS(i,j,k,bi,bj)
288    #ifdef ALLOW_OBCS
289         &                  *maskInS(i,j,bi,bj)
290    #endif
291            ENDDO            ENDDO
292           ENDDO           ENDDO
293  C-       gradient in y direction:  C-       gradient in y direction:
294           DO j=2-Oly,sNy+Oly-1           DO j=2-OLy,sNy+OLy-1
295            DO i=2-Olx,sNx+Olx            DO i=2-OLx,sNx+OLx
296              vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))              vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
297       &                  *recip_DYG(i,j,bi,bj)       &                  *recip_dyG(i,j,bi,bj)
298       &                  *maskW(i,j,k,bi,bj)       &                  *maskW(i,j,k,bi,bj)
299    #ifdef ALLOW_OBCS
300         &                  *maskInW(i,j,bi,bj)
301    #endif
302            ENDDO            ENDDO
303           ENDDO           ENDDO
304    
305    C--   end if calcLeith
306         ENDIF         ENDIF
307    
308         DO j=2-Oly,sNy+Oly-1         DO j=2-OLy,sNy+OLy-1
309          DO i=2-Olx,sNx+Olx-1          DO i=2-OLx,sNx+OLx-1
310  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
311    
312  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
313  # ifndef AUTODIFF_DISABLE_LEITH  # ifndef AUTODIFF_DISABLE_LEITH
314             lockey_2 = i+olx + (sNx+2*olx)*(j+oly-1)             lockey_2 = i+olx + (sNx+2*olx)*(j+oly-1)
315       &                      + (sNx+2*olx)*(sNy+2*oly)*(lockey_1-1)       &                      + (sNx+2*olx)*(sNy+2*oly)*(lockey_1-1)
316  CADJ STORE viscA4_ZSmg(i,j)  CADJ STORE viscA4_ZSmg(i,j)
317  CADJ &     = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte  CADJ &     = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte
318  CADJ STORE viscAh_ZSmg(i,j)  CADJ STORE viscAh_ZSmg(i,j)
319  CADJ &     = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte  CADJ &     = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte
320  # endif  # endif
321  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
322    
323  C These are (powers of) length scales  C These are (powers of) length scales
324           IF (useAreaViscLength) THEN           L2 = L2_D(i,j,bi,bj)
325            L2=rA(i,j,bi,bj)           L2rdt = 0.25 _d 0*recip_dt*L2
326            L4rdt=0.03125 _d 0*recip_dt*L2**2           L3 = L3_D(i,j,bi,bj)
327           ELSE           L4rdt = L4rdt_D(i,j,bi,bj)
328            L2=2. _d 0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))           L5 = (L2*L3)
           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  
          L3=(L2**1.5)  
          L4=(L2**2)  
          L5=(L2*L3)  
   
          L2rdt=0.25 _d 0*recip_dt*L2  
329    
330    #ifndef AUTODIFF_DISABLE_REYNOLDS_SCALE
331  C Velocity Reynolds Scale  C Velocity Reynolds Scale
332           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
333             Uscl=SQRT(KE(i,j)*L2)*viscAhRe_max             Uscl=SQRT(KE(i,j)*L2)*viscAhRe_max
# Line 371  C Velocity Reynolds Scale Line 339  C Velocity Reynolds Scale
339           ELSE           ELSE
340             U4scl=0.             U4scl=0.
341           ENDIF           ENDIF
342    #endif /* ndef AUTODIFF_DISABLE_REYNOLDS_SCALE */
343    
 cph-leith#ifndef ALLOW_AUTODIFF_TAMC  
344  #ifndef AUTODIFF_DISABLE_LEITH  #ifndef AUTODIFF_DISABLE_LEITH
345           IF (useFullLeith.AND.calcLeith) THEN           IF (useFullLeith.AND.calcLeith) THEN
346  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
# Line 397  C Using it in Leith serves to damp insta Line 365  C Using it in Leith serves to damp insta
365            viscA4_DLthd(i,j)=            viscA4_DLthd(i,j)=
366       &     SQRT(leithD4fac*grdDiv)*L5       &     SQRT(leithD4fac*grdDiv)*L5
367           ELSEIF (calcLeith) THEN           ELSEIF (calcLeith) THEN
368  C but this approximation will work on cube  C but this approximation will work on cube (and differs by as much as 4X)
 c (and differs by as much as 4X)  
369            grdVrt=MAX( ABS(vrtDx(i,j+1)), ABS(vrtDx(i,j)) )            grdVrt=MAX( ABS(vrtDx(i,j+1)), ABS(vrtDx(i,j)) )
370            grdVrt=MAX( grdVrt, ABS(vrtDy(i+1,j)) )            grdVrt=MAX( grdVrt, ABS(vrtDy(i+1,j)) )
371            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j))   )            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j))   )
372    
373  c This approximation is good to the same order as above...  C This approximation is good to the same order as above...
374            grdDiv=MAX( ABS(divDx(i+1,j)), ABS(divDx(i,j)) )            grdDiv=MAX( ABS(divDx(i+1,j)), ABS(divDx(i,j)) )
375            grdDiv=MAX( grdDiv, ABS(divDy(i,j+1)) )            grdDiv=MAX( grdDiv, ABS(divDy(i,j+1)) )
376            grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )            grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
# Line 435  c This approximation is good to the same Line 402  c This approximation is good to the same
402  C  Harmonic on Div.u points  C  Harmonic on Div.u points
403           Alin=viscAhD+viscAhGrid*L2rdt           Alin=viscAhD+viscAhGrid*L2rdt
404       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
405    #ifdef ALLOW_3D_VISCAH
406         &          +viscAhDfld(i,j,k,bi,bj)
407    #ifdef ALLOW_AUTODIFF
408         &          *viscFacAdj
409    #endif
410    #endif
411           viscAh_DMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)           viscAh_DMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
412           viscAh_D(i,j)=MAX(viscAh_DMin(i,j),Alin)           viscAh_D(i,j)=MAX(viscAh_DMin(i,j),Alin)
413           viscAh_DMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)           viscAh_DMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
# Line 443  C  Harmonic on Div.u points Line 416  C  Harmonic on Div.u points
416  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
417           Alin=viscA4D+viscA4Grid*L4rdt           Alin=viscA4D+viscA4Grid*L4rdt
418       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
419    #ifdef ALLOW_3D_VISCA4
420         &          +viscA4Dfld(i,j,k,bi,bj)
421    #ifdef ALLOW_AUTODIFF
422         &          *viscFacAdj
423    #endif
424    #endif
425           viscA4_DMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)           viscA4_DMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
426           viscA4_D(i,j)=MAX(viscA4_DMin(i,j),Alin)           viscA4_D(i,j)=MAX(viscA4_DMin(i,j),Alin)
427           viscA4_DMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)           viscA4_DMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
428           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))
429    
 #ifdef ALLOW_NONHYDROSTATIC  
 C--   Pass Viscosities to calc_gw, if constant, not necessary  
   
          kp1  = MIN(k+1,Nr)  
   
          IF ( k.EQ.1 ) THEN  
 C     Prepare for next level (next call)  
           viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)  
           viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)  
   
 C     These values dont get used  
           viscAh_W(i,j,k,bi,bj)=viscAh_D(i,j)  
           viscA4_W(i,j,k,bi,bj)=viscA4_D(i,j)  
   
          ELSEIF ( k.EQ.Nr ) THEN  
           viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)  
           viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)  
   
          ELSE  
 C     Prepare for next level (next call)  
           viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)  
           viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)  
   
 C     Note that previous call of this function has already added half.  
           viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)  
           viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)  
   
          ENDIF  
 #endif /* ALLOW_NONHYDROSTATIC */  
   
430  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
431  C These are (powers of) length scales  C These are (powers of) length scales
432           IF (useAreaViscLength) THEN           L2 = L2_Z(i,j,bi,bj)
433            L2=rAz(i,j,bi,bj)           L2rdt = 0.25 _d 0*recip_dt*L2
434            L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2           L3 = L3_Z(i,j,bi,bj)
435           ELSE           L4rdt = L4rdt_Z(i,j,bi,bj)
436            L2=2. _d 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))           L5 = (L2*L3)
           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  
   
          L3=(L2**1.5)  
          L4=(L2**2)  
          L5=(L2*L3)  
   
          L2rdt=0.25 _d 0*recip_dt*L2  
437    
438    #ifndef AUTODIFF_DISABLE_REYNOLDS_SCALE
439  C Velocity Reynolds Scale (Pb here at CS-grid corners !)  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
440           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
441             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))
# Line 511  C Velocity Reynolds Scale (Pb here at CS Line 451  C Velocity Reynolds Scale (Pb here at CS
451             Uscl =0.             Uscl =0.
452             U4scl=0.             U4scl=0.
453           ENDIF           ENDIF
454    #endif /* ndef AUTODIFF_DISABLE_REYNOLDS_SCALE */
455    
456  #ifndef AUTODIFF_DISABLE_LEITH  #ifndef AUTODIFF_DISABLE_LEITH
457  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
# Line 569  C but this approximation will work on cu Line 510  C but this approximation will work on cu
510  C  Harmonic on Zeta points  C  Harmonic on Zeta points
511           Alin=viscAhZ+viscAhGrid*L2rdt           Alin=viscAhZ+viscAhGrid*L2rdt
512       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
513    #ifdef ALLOW_3D_VISCAH
514         &          +viscAhZfld(i,j,k,bi,bj)
515    #endif
516           viscAh_ZMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)           viscAh_ZMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
517           viscAh_Z(i,j)=MAX(viscAh_ZMin(i,j),Alin)           viscAh_Z(i,j)=MAX(viscAh_ZMin(i,j),Alin)
518           viscAh_ZMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)           viscAh_ZMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
# Line 577  C  Harmonic on Zeta points Line 521  C  Harmonic on Zeta points
521  C  BiHarmonic on Zeta points  C  BiHarmonic on Zeta points
522           Alin=viscA4Z+viscA4Grid*L4rdt           Alin=viscA4Z+viscA4Grid*L4rdt
523       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
524    #ifdef ALLOW_3D_VISCA4
525         &          +viscA4Zfld(i,j,k,bi,bj)
526    #endif
527           viscA4_ZMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)           viscA4_ZMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
528           viscA4_Z(i,j)=MAX(viscA4_ZMin(i,j),Alin)           viscA4_Z(i,j)=MAX(viscA4_ZMin(i,j),Alin)
529           viscA4_ZMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)           viscA4_ZMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
# Line 584  C  BiHarmonic on Zeta points Line 531  C  BiHarmonic on Zeta points
531          ENDDO          ENDDO
532         ENDDO         ENDDO
533    
534        ELSE  #ifdef ALLOW_NONHYDROSTATIC
535  C---- use constant viscosity (useVariableViscosity=F):         IF ( nonHydrostatic ) THEN
536    C--   Pass Viscosities to calc_gw (if constant, not necessary)
537    
538         DO j=1-Oly,sNy+Oly          IF ( k.LT.Nr ) THEN
539          DO i=1-Olx,sNx+Olx  C     Prepare for next level (next call)
540           viscAh_D(i,j)=viscAhD           DO j=1-OLy,sNy+OLy
541           viscAh_Z(i,j)=viscAhZ            DO i=1-OLx,sNx+OLx
542           viscA4_D(i,j)=viscA4D              viscAh_W(i,j,k+1,bi,bj) = halfRL*viscAh_D(i,j)
543           viscA4_Z(i,j)=viscA4Z              viscA4_W(i,j,k+1,bi,bj) = halfRL*viscA4_D(i,j)
544          ENDDO            ENDDO
545         ENDDO           ENDDO
546            ENDIF
547    
548            shiftAh = viscAhW - viscAhD
549            shiftA4 = viscA4W - viscA4D
550            IF ( k.EQ.1 ) THEN
551    C     These values dont get used
552             DO j=1-OLy,sNy+OLy
553              DO i=1-OLx,sNx+OLx
554                viscAh_W(i,j,k,bi,bj) = shiftAh + viscAh_D(i,j)
555                viscA4_W(i,j,k,bi,bj) = shiftA4 + viscA4_D(i,j)
556              ENDDO
557             ENDDO
558            ELSE
559    C     Note that previous call of this function has already added half.
560             DO j=1-OLy,sNy+OLy
561              DO i=1-OLx,sNx+OLx
562                viscAh_W(i,j,k,bi,bj) = shiftAh + viscAh_W(i,j,k,bi,bj)
563         &                                      + halfRL*viscAh_D(i,j)
564                viscA4_W(i,j,k,bi,bj) = shiftA4 + viscA4_W(i,j,k,bi,bj)
565         &                                      + halfRL*viscA4_D(i,j)
566              ENDDO
567             ENDDO
568            ENDIF
569    
570           ENDIF
571    #endif /* ALLOW_NONHYDROSTATIC */
572    
573    c     ELSE
574    C---- use constant viscosity (useVariableVisc=F):
575    c      DO j=1-OLy,sNy+OLy
576    c       DO i=1-OLx,sNx+OLx
577    c        viscAh_D(i,j) = viscAhD
578    c        viscAh_Z(i,j) = viscAhZ
579    c        viscA4_D(i,j) = viscA4D
580    c        viscA4_Z(i,j) = viscA4Z
581    c       ENDDO
582    c      ENDDO
583  C---- variable/constant viscosity : end if/else block  C---- variable/constant viscosity : end if/else block
584        ENDIF  c     ENDIF
585    
586  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
587        IF (useDiagnostics) THEN        IF (useDiagnostics) THEN
# Line 621  C---- variable/constant viscosity : end Line 605  C---- variable/constant viscosity : end
605         CALL DIAGNOSTICS_FILL(viscAh_ZLth,'VAHZLTH ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_ZLth,'VAHZLTH ',k,1,2,bi,bj,myThid)
606         CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_ZLth,'VA4ZLTH ',k,1,2,bi,bj,myThid)
607    
608         CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD'         CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD',
609       &   ,k,1,2,bi,bj,myThid)       &                       k,1,2,bi,bj,myThid)
610         CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD'         CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD',
611       &   ,k,1,2,bi,bj,myThid)       &                       k,1,2,bi,bj,myThid)
612         CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD'         CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD',
613       &   ,k,1,2,bi,bj,myThid)       &                       k,1,2,bi,bj,myThid)
614         CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD'         CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD',
615       &   ,k,1,2,bi,bj,myThid)       &                       k,1,2,bi,bj,myThid)
616    
617         CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_DSmg,'VAHDSMAG',k,1,2,bi,bj,myThid)
618         CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_DSmg,'VA4DSMAG',k,1,2,bi,bj,myThid)
# Line 639  C---- variable/constant viscosity : end Line 623  C---- variable/constant viscosity : end
623    
624        RETURN        RETURN
625        END        END
   
Legend:
Removed from v.1.37  
changed lines
  Added in v.1.48
  ViewVC Help
Powered by ViewVC 1.1.22