/[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.1 by baylor, Fri Sep 16 19:33:59 2005 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,
17       O        harmonic,biharmonic,useVariableViscosity,       I        hDiv,vort3,tension,strain,KE,hFacZ,
      I        hDiv,vort3,tension,strain,KE,  
18       I        myThid)       I        myThid)
19    
20    C     !DESCRIPTION:
21    C     Calculate horizontal viscosities (L is typical grid width)
22    C     harmonic viscosity=
23    C       viscAh (or viscAhD on div pts and viscAhZ on zeta pts)
24    C       +0.25*L**2*viscAhGrid/deltaT
25    C       +sqrt((viscC2leith/pi)**6*grad(Vort3)**2
26    C             +(viscC2leithD/pi)**6*grad(hDiv)**2)*L**3
27    C       +(viscC2smag/pi)**2*L**2*sqrt(Tension**2+Strain**2)
28    C
29    C     biharmonic viscosity=
30    C       viscA4 (or viscA4D on div pts and viscA4Z on zeta pts)
31    C       +0.25*0.125*L**4*viscA4Grid/deltaT (approx)
32    C       +0.125*L**5*sqrt((viscC4leith/pi)**6*grad(Vort3)**2
33    C                        +(viscC4leithD/pi)**6*grad(hDiv)**2)
34    C       +0.125*L**4*(viscC4smag/pi)**2*sqrt(Tension**2+Strain**2)
35    C
36    C     Note that often 0.125*L**2 is the scale between harmonic and
37    C     biharmonic (see Griffies and Hallberg (2000))
38    C     This allows the same value of the coefficient to be used
39    C     for roughly similar results with biharmonic and harmonic
40    C
41    C     LIMITERS -- limit min and max values of viscosities
42    C     viscAhReMax is min value for grid point harmonic Reynolds num
43    C      harmonic viscosity>sqrt(2*KE)*L/viscAhReMax
44    C
45    C     viscA4ReMax is min value for grid point biharmonic Reynolds num
46    C      biharmonic viscosity>sqrt(2*KE)*L**3/8/viscA4ReMax
47    C
48    C     viscAhgridmax is CFL stability limiter for harmonic viscosity
49    C      harmonic viscosity<0.25*viscAhgridmax*L**2/deltaT
50    C
51    C     viscA4gridmax is CFL stability limiter for biharmonic viscosity
52    C      biharmonic viscosity<viscA4gridmax*L**4/32/deltaT (approx)
53    C
54    C     viscAhgridmin and viscA4gridmin are lower limits for viscosity:
55    C       harmonic viscosity>0.25*viscAhgridmin*L**2/deltaT
56    C       biharmonic viscosity>viscA4gridmin*L**4/32/deltaT (approx)
57    
58    C     RECOMMENDED VALUES
59    C     viscC2Leith=1-3
60    C     viscC2LeithD=1-3
61    C     viscC4Leith=1-3
62    C     viscC4LeithD=1.5-3
63    C     viscC2smag=2.2-4 (Griffies and Hallberg,2000)
64    C               0.2-0.9 (Smagorinsky,1993)
65    C     viscC4smag=2.2-4 (Griffies and Hallberg,2000)
66    C     viscAhReMax>=1, (<2 suppresses a computational mode)
67    C     viscA4ReMax>=1, (<2 suppresses a computational mode)
68    C     viscAhgridmax=1
69    C     viscA4gridmax=1
70    C     viscAhgrid<1
71    C     viscA4grid<1
72    C     viscAhgridmin<<1
73    C     viscA4gridmin<<1
74    
75    C     !USES:
76        IMPLICIT NONE        IMPLICIT NONE
77    
78  C     == Global variables ==  C     == Global variables ==
# Line 17  C     == Global variables == Line 80  C     == Global variables ==
80  #include "GRID.h"  #include "GRID.h"
81  #include "EEPARAMS.h"  #include "EEPARAMS.h"
82  #include "PARAMS.h"  #include "PARAMS.h"
83    #include "MOM_VISC.h"
84    #ifdef ALLOW_AUTODIFF
85    #include "tamc.h"
86    #include "tamc_keys.h"
87    #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 31  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        _RL ASmag2, ASmag4, smag2fac, smag4fac  #ifdef ALLOW_NONHYDROSTATIC
108        _RL vg2Min, vg2Max, AlinMax, AlinMin        _RL shiftAh, shiftA4
109        _RL lenA2, lenAz2  #endif
110        _RL Alin,Alth2,Alth4,grdVrt,grdDiv  #ifdef ALLOW_AUTODIFF_TAMC
111        _RL vg2,vg4,vg4Min,vg4Max        INTEGER lockey_1, lockey_2
112        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt  #endif
113          _RL smag2fac, smag4fac
114          _RL leith2fac, leith4fac
115        useVariableViscosity=        _RL leithD2fac, leithD4fac
116       &      (viscAhGrid.NE.0.)        _RL viscAhRe_max, viscA4Re_max
117       &  .OR.(viscA4Grid.NE.0.)        _RL Alin,grdVrt,grdDiv, keZpt
118       &  .OR.(viscC2leith.NE.0.)        _RL L2, L3, L5, L2rdt, L4rdt, recip_dt
119          _RL Uscl,U4scl
120          _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121          _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122          _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123          _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124          _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125          _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126          _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127          _RL viscA4_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128          _RL viscAh_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129          _RL viscAh_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130          _RL viscA4_ZMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131          _RL viscA4_DMin(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132          _RL viscAh_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133          _RL viscAh_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134          _RL viscA4_ZLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135          _RL viscA4_DLth(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136          _RL viscAh_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137          _RL viscAh_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138          _RL viscA4_ZLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139          _RL viscA4_DLthD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140          _RL viscAh_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141          _RL viscAh_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
142          _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
143          _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
144          LOGICAL calcLeith, calcSmag
145    
146    #ifdef ALLOW_AUTODIFF_TAMC
147              act1 = bi - myBxLo(myThid)
148              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
149              act2 = bj - myByLo(myThid)
150              max2 = myByHi(myThid) - myByLo(myThid) + 1
151              act3 = myThid - 1
152              max3 = nTx*nTy
153              act4 = ikey_dynamics - 1
154              ikey = (act1 + 1) + act2*max1
155         &                      + act3*max1*max2
156         &                      + act4*max1*max2*max3
157              lockey_1 = (ikey-1)*Nr + k
158    #endif /* ALLOW_AUTODIFF_TAMC */
159    
160    C--   Set flags which are used in this S/R and elsewhere :
161    C     useVariableVisc, useHarmonicVisc and useBiharmonicVisc
162    C     are now set early on (in S/R SET_PARAMS)
163    
164    c     IF ( useVariableVisc ) THEN
165    C---- variable viscosity :
166    
167           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
172           ELSE
173            viscAhRe_max=0. _d 0
174           ENDIF
175    
176           IF ( useBiharmonicVisc .AND. viscA4ReMax.NE.0. ) THEN
177            viscA4Re_max=0.125 _d 0*SQRT(2. _d 0)/viscA4ReMax
178           ELSE
179            viscA4Re_max=0. _d 0
180           ENDIF
181    
182           calcLeith=
183         &      (viscC2leith.NE.0.)
184       &  .OR.(viscC2leithD.NE.0.)       &  .OR.(viscC2leithD.NE.0.)
185       &  .OR.(viscC4leith.NE.0.)       &  .OR.(viscC4leith.NE.0.)
186       &  .OR.(viscC4leithD.NE.0.)       &  .OR.(viscC4leithD.NE.0.)
187       &  .OR.(viscC2smag.NE.0.)  
188           calcSmag=
189         &      (viscC2smag.NE.0.)
190       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
191    
192        harmonic=         IF (calcSmag) THEN
193       &      (viscAh.NE.0.)          smag2fac=(viscC2smag/pi)**2
194       &  .OR.(viscAhD.NE.0.)          smag4fac=0.125 _d 0*(viscC4smag/pi)**2
195       &  .OR.(viscAhZ.NE.0.)         ELSE
196       &  .OR.(viscAhGrid.NE.0.)          smag2fac=0. _d 0
197       &  .OR.(viscC2leith.NE.0.)          smag4fac=0. _d 0
198       &  .OR.(viscC2leithD.NE.0.)         ENDIF
199       &  .OR.(viscC2smag.NE.0.)  
200           IF (calcLeith) THEN
201            IF (useFullLeith) THEN
202             leith2fac =(viscC2leith /pi)**6
203             leithD2fac=(viscC2leithD/pi)**6
204             leith4fac =0.015625 _d 0*(viscC4leith /pi)**6
205             leithD4fac=0.015625 _d 0*(viscC4leithD/pi)**6
206            ELSE
207             leith2fac =(viscC2leith /pi)**3
208             leithD2fac=(viscC2leithD/pi)**3
209             leith4fac =0.125 _d 0*(viscC4leith /pi)**3
210             leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3
211            ENDIF
212           ELSE
213            leith2fac=0. _d 0
214            leith4fac=0. _d 0
215            leithD2fac=0. _d 0
216            leithD4fac=0. _d 0
217           ENDIF
218    
219           DO j=1-OLy,sNy+OLy
220            DO i=1-OLx,sNx+OLx
221    C-    viscosity arrays have been initialised everywhere before calling this S/R
222    c         viscAh_D(i,j) = viscAhD
223    c         viscAh_Z(i,j) = viscAhZ
224    c         viscA4_D(i,j) = viscA4D
225    c         viscA4_Z(i,j) = viscA4Z
226    
227              visca4_zsmg(i,j) = 0. _d 0
228              viscah_zsmg(i,j) = 0. _d 0
229    
230              viscAh_Dlth(i,j) = 0. _d 0
231              viscA4_Dlth(i,j) = 0. _d 0
232              viscAh_DlthD(i,j)= 0. _d 0
233              viscA4_DlthD(i,j)= 0. _d 0
234    
235              viscAh_DSmg(i,j) = 0. _d 0
236              viscA4_DSmg(i,j) = 0. _d 0
237    
238              viscAh_ZLth(i,j) = 0. _d 0
239              viscA4_ZLth(i,j) = 0. _d 0
240              viscAh_ZLthD(i,j)= 0. _d 0
241              viscA4_ZLthD(i,j)= 0. _d 0
242            ENDDO
243           ENDDO
244    
245        biharmonic=  C-    Initialise to zero gradient of vorticity & divergence:
246       &      (viscA4.NE.0.)         DO j=1-OLy,sNy+OLy
247       &  .OR.(viscA4D.NE.0.)          DO i=1-OLx,sNx+OLx
248       &  .OR.(viscA4Z.NE.0.)            divDx(i,j) = 0.
249       &  .OR.(viscA4Grid.NE.0.)            divDy(i,j) = 0.
250       &  .OR.(viscC4leith.NE.0.)            vrtDx(i,j) = 0.
251       &  .OR.(viscC4leithD.NE.0.)            vrtDy(i,j) = 0.
252       &  .OR.(viscC4smag.NE.0.)          ENDDO
253           ENDDO
254    
255        IF (deltaTmom.NE.0.) THEN         IF ( calcLeith ) THEN
256         recip_dt=1./deltaTmom  C--   horizontal gradient of horizontal divergence:
257        ELSE  C-       gradient in x direction:
258         recip_dt=0.           IF (useCubedSphereExchange) THEN
259        ENDIF  C        to compute d/dx(hDiv), fill corners with appropriate values:
260               CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
261         &                                hDiv, bi,bj, myThid )
262             ENDIF
263             DO j=2-OLy,sNy+OLy-1
264              DO i=2-OLx,sNx+OLx-1
265                divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_dxC(i,j,bi,bj)
266              ENDDO
267             ENDDO
268    
269    C-       gradient in y direction:
270             IF (useCubedSphereExchange) THEN
271    C        to compute d/dy(hDiv), fill corners with appropriate values:
272               CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
273         &                                hDiv, bi,bj, myThid )
274             ENDIF
275             DO j=2-OLy,sNy+OLy-1
276              DO i=2-OLx,sNx+OLx-1
277                divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_dyC(i,j,bi,bj)
278              ENDDO
279             ENDDO
280    
281    C--   horizontal gradient of vertical vorticity:
282    C-       gradient in x direction:
283             DO j=2-OLy,sNy+OLy
284              DO i=2-OLx,sNx+OLx-1
285                vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
286         &                  *recip_dxG(i,j,bi,bj)
287         &                  *maskS(i,j,k,bi,bj)
288    #ifdef ALLOW_OBCS
289         &                  *maskInS(i,j,bi,bj)
290    #endif
291              ENDDO
292             ENDDO
293    C-       gradient in y direction:
294             DO j=2-OLy,sNy+OLy-1
295              DO i=2-OLx,sNx+OLx
296                vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
297         &                  *recip_dyG(i,j,bi,bj)
298         &                  *maskW(i,j,k,bi,bj)
299    #ifdef ALLOW_OBCS
300         &                  *maskInW(i,j,bi,bj)
301    #endif
302              ENDDO
303             ENDDO
304    
305    C--   end if calcLeith
306           ENDIF
307    
308        vg2=viscAhGrid*recip_dt         DO j=2-OLy,sNy+OLy-1
309        vg2Min=viscAhGridMin*recip_dt          DO i=2-OLx,sNx+OLx-1
       vg2Max=viscAhGridMax*recip_dt  
       vg4=viscA4Grid*recip_dt  
       vg4Min=viscA4GridMin*recip_dt  
       vg4Max=viscA4GridMax*recip_dt  
   
       smag2fac=(viscC2smag/pi)**2  
       smag4fac=0.125*(viscC4smag/pi)**2  
   
 C     - Viscosity  
       IF (useVariableViscosity) THEN  
        DO j=2-Oly,sNy+Oly-1  
         DO i=2-Olx,sNx+Olx-1  
310  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
311  C These are (powers of) length scales  
312           L2=rA(i,j,bi,bj)  #ifdef ALLOW_AUTODIFF_TAMC
313           L3=(L2**1.5)  # ifndef AUTODIFF_DISABLE_LEITH
314           L4=(L2**2)             lockey_2 = i+olx + (sNx+2*olx)*(j+oly-1)
315           L5=0.125*(L2**2.5)       &                      + (sNx+2*olx)*(sNy+2*oly)*(lockey_1-1)
316           IF (useAnisotropicViscAGridMax) THEN  CADJ STORE viscA4_ZSmg(i,j)
317            L2rdt=recip_dt/( 2.*(recip_DXF(I,J,bi,bj)**2  CADJ &     = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte
318       &                       +recip_DYF(I,J,bi,bj)**2) )  CADJ STORE viscAh_ZSmg(i,j)
319            L4rdt=recip_dt/( 6.*(recip_DXF(I,J,bi,bj)**4  CADJ &     = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte
320       &                       +recip_DYF(I,J,bi,bj)**4)  # endif
321       &                   +8.*((recip_DXF(I,J,bi,bj)  #endif /* ALLOW_AUTODIFF_TAMC */
322       &                        *recip_DYF(I,J,bi,bj))**2) )  
323    C These are (powers of) length scales
324             L2 = L2_D(i,j,bi,bj)
325             L2rdt = 0.25 _d 0*recip_dt*L2
326             L3 = L3_D(i,j,bi,bj)
327             L4rdt = L4rdt_D(i,j,bi,bj)
328             L5 = (L2*L3)
329    
330    #ifndef AUTODIFF_DISABLE_REYNOLDS_SCALE
331    C Velocity Reynolds Scale
332             IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
333               Uscl=SQRT(KE(i,j)*L2)*viscAhRe_max
334             ELSE
335               Uscl=0.
336           ENDIF           ENDIF
337             IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
338               U4scl=SQRT(KE(i,j))*L3*viscA4Re_max
339             ELSE
340               U4scl=0.
341             ENDIF
342    #endif /* ndef AUTODIFF_DISABLE_REYNOLDS_SCALE */
343    
344           IF (useFullLeith) THEN  #ifndef AUTODIFF_DISABLE_LEITH
345             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
347            grdVrt=0.25*(            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
348       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
349       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
350       &     +((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)  
351    
352  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
353  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
354             grdDiv=0.25*(            grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j)
355       &      ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
356       &      +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (divDy(i,j+1)*divDy(i,j+1)
357       &      +((hDiv(i-1,j)-hDiv(i,j))*recip_DXG(i-1,j,bi,bj))**2       &                        + divDy(i,j)*divDy(i,j) )  )
358       &      +((hDiv(i,j-1)-hDiv(i,j))*recip_DYG(i,j-1,bi,bj))**2)  
359              viscAh_DLth(i,j)=
360             IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)       &     SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
361       &          .NE. 0. ) THEN            viscA4_DLth(i,j)=
362              Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
363             ELSE            viscAh_DLthd(i,j)=
364              Alth2=0. _d 0       &     SQRT(leithD2fac*grdDiv)*L3
365             ENDIF            viscA4_DLthd(i,j)=
366             IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)       &     SQRT(leithD4fac*grdDiv)*L5
367       &          .NE. 0. ) THEN           ELSEIF (calcLeith) THEN
368              Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5  C but this approximation will work on cube (and differs by as much as 4X)
369             ELSE            grdVrt=MAX( ABS(vrtDx(i,j+1)), ABS(vrtDx(i,j)) )
370              Alth4=0. _d 0            grdVrt=MAX( grdVrt, ABS(vrtDy(i+1,j)) )
371             ENDIF            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j))   )
372           ELSE  
373  C but this approximation will work on cube  C This approximation is good to the same order as above...
374  c (and differs by as much as 4X)            grdDiv=MAX( ABS(divDx(i+1,j)), ABS(divDx(i,j)) )
375             grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))            grdDiv=MAX( grdDiv, ABS(divDy(i,j+1)) )
376             grdVrt=max(grdVrt,            grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
377       &      abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))  
378             grdVrt=max(grdVrt,            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
379       &      abs((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj)))            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
380             grdVrt=max(grdVrt,            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
381       &      abs((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj)))            viscA4_DlthD(i,j)=((leithD4fac*grdDiv))*L5
   
            grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj)))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i-1,j)-hDiv(i,j))*recip_DXG(i-1,j,bi,bj)))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i,j-1)-hDiv(i,j))*recip_DYG(i,j-1,bi,bj)))  
   
 c This approximation is good to the same order as above...  
            Alth2=(viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3  
            Alth4=(viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5  
          ENDIF  
   
          IF (smag2fac.NE.0.) THEN  
           Asmag2=smag2fac*L2  
      &           *sqrt(tension(i,j)**2  
      &           +0.25*(strain(i+1, j )**2+strain( i ,j+1)**2  
      &                 +strain(i-1, j )**2+strain( i ,j-1)**2))  
382           ELSE           ELSE
383            Asmag2=0d0            viscAh_Dlth(i,j)=0. _d 0
384              viscA4_Dlth(i,j)=0. _d 0
385              viscAh_DlthD(i,j)=0. _d 0
386              viscA4_DlthD(i,j)=0. _d 0
387           ENDIF           ENDIF
388    
389           IF (smag4fac.NE.0.) THEN           IF (calcSmag) THEN
390            Asmag4=smag4fac*L4            viscAh_DSmg(i,j)=L2
391       &             *sqrt(tension(i,j)**2       &       *SQRT(tension(i,j)**2
392       &             +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
393       &                   +strain(i-1, j )**2+strain( i ,j-1)**2))       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))
394              viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
395              viscAh_DSmg(i,j)=smag2fac*viscAh_DSmg(i,j)
396           ELSE           ELSE
397            Asmag4=0d0            viscAh_DSmg(i,j)=0. _d 0
398              viscA4_DSmg(i,j)=0. _d 0
399           ENDIF           ENDIF
400    #endif /* AUTODIFF_DISABLE_LEITH */
401    
402  C  Harmonic on Div.u points  C  Harmonic on Div.u points
403           Alin=viscAhD+vg2*L2+Alth2+Asmag2           Alin=viscAhD+viscAhGrid*L2rdt
404           viscAh_D(i,j)=min(viscAhMax,Alin)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
405           IF (useAnisotropicViscAGridMax) THEN  #ifdef ALLOW_3D_VISCAH
406            AlinMax=viscAhGridMax*L2rdt       &          +viscAhDfld(i,j,k,bi,bj)
407            viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))  #ifdef ALLOW_AUTODIFF
408           ELSE       &          *viscFacAdj
409            IF (vg2Max.GT.0.) THEN  #endif
410             AlinMax=vg2Max*L2  #endif
411             viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))           viscAh_DMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
412            ENDIF           viscAh_D(i,j)=MAX(viscAh_DMin(i,j),Alin)
413           ENDIF           viscAh_DMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
414           AlinMin=vg2Min*L2           viscAh_D(i,j)=MIN(viscAh_DMax(i,j),viscAh_D(i,j))
          viscAh_D(i,j)=max(AlinMin,viscAh_D(i,j))  
415    
416  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
417           Alin=viscA4D+vg4*L4+Alth4+Asmag4           Alin=viscA4D+viscA4Grid*L4rdt
418           viscA4_D(i,j)=min(viscA4Max,Alin)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
419           IF (useAnisotropicViscAGridMax) THEN  #ifdef ALLOW_3D_VISCA4
420              AlinMax=viscA4GridMax*L4rdt       &          +viscA4Dfld(i,j,k,bi,bj)
421              viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))  #ifdef ALLOW_AUTODIFF
422           ELSE       &          *viscFacAdj
423           IF (vg4Max.GT.0.) THEN  #endif
424              AlinMax=vg4Max*L4  #endif
425              viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))           viscA4_DMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
426           ENDIF           viscA4_D(i,j)=MAX(viscA4_DMin(i,j),Alin)
427           ENDIF           viscA4_DMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
428           AlinMin=vg4Min*L4           viscA4_D(i,j)=MIN(viscA4_DMax(i,j),viscA4_D(i,j))
          viscA4_D(i,j)=max(AlinMin,viscA4_D(i,j))  
429    
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           L2=rAz(i,j,bi,bj)           L2 = L2_Z(i,j,bi,bj)
433           L3=(L2**1.5)           L2rdt = 0.25 _d 0*recip_dt*L2
434           L4=(L2**2)           L3 = L3_Z(i,j,bi,bj)
435           L5=0.125*(L2**2.5)           L4rdt = L4rdt_Z(i,j,bi,bj)
436           IF (useAnisotropicViscAGridMax) THEN           L5 = (L2*L3)
437           L2rdt=recip_dt/( 2.*(recip_DXV(I,J,bi,bj)**2  
438       &                       +recip_DYU(I,J,bi,bj)**2) )  #ifndef AUTODIFF_DISABLE_REYNOLDS_SCALE
439           L4rdt=recip_dt/( 6.*(recip_DXV(I,J,bi,bj)**4  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
440       &                       +recip_DYU(I,J,bi,bj)**4)           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
441       &                   +8.*((recip_DXV(I,J,bi,bj)             keZpt=0.25 _d 0*( (KE(i,j)+KE(i-1,j-1))
442       &                        *recip_DYU(I,J,bi,bj))**2) )       &                      +(KE(i-1,j)+KE(i,j-1)) )
443               IF ( keZpt.GT.0. ) THEN
444                 Uscl = SQRT(keZpt*L2)*viscAhRe_max
445                 U4scl= SQRT(keZpt)*L3*viscA4Re_max
446               ELSE
447                 Uscl =0.
448                 U4scl=0.
449               ENDIF
450             ELSE
451               Uscl =0.
452               U4scl=0.
453           ENDIF           ENDIF
454    #endif /* ndef AUTODIFF_DISABLE_REYNOLDS_SCALE */
455    
456    #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
458           IF (useFullLeith) THEN           IF (useFullLeith.AND.calcLeith) THEN
459             grdVrt=0.25*(            grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
460       &      ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
461       &      +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
462       &      +((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj))**2       &                        + vrtDy(i,j)*vrtDy(i,j) )  )
      &      +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)  
463    
464  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
465             grdDiv=0.25*(            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
466       &       ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
467       &      +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (divDy(i-1,j)*divDy(i-1,j)
468       &      +((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i,j+1,bi,bj))**2       &                        + divDy(i,j)*divDy(i,j) )  )
469       &      +((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)  
470              viscAh_ZLth(i,j)=
471             IF ( (viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)       &     SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
472       &          .NE. 0. ) THEN            viscA4_ZLth(i,j)=
473              Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
474             ELSE            viscAh_ZLthD(i,j)=
475              Alth2=0. _d 0       &     SQRT(leithD2fac*grdDiv)*L3
476             ENDIF            viscA4_ZLthD(i,j)=
477             IF ( (viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)       &     SQRT(leithD4fac*grdDiv)*L5
      &          .NE. 0. ) THEN  
             Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5  
            ELSE  
             Alth4=0. _d 0  
            ENDIF  
          ELSE  
 C but this approximation will work on cube (and differs by 4X)  
            grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))  
            grdVrt=max(grdVrt,  
      &       abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))  
            grdVrt=max(grdVrt,  
      &       abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))  
            grdVrt=max(grdVrt,  
      &       abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))  
   
            grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj)))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i-1,j,bi,bj)))  
            grdDiv=max(grdDiv,  
      &      abs((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i,j-1,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  
478    
479           IF (smag2fac.NE.0.) THEN           ELSEIF (calcLeith) THEN
480            Asmag2=smag2fac*L2  C but this approximation will work on cube (and differs by 4X)
481       &           *sqrt(strain(i,j)**2            grdVrt=MAX( ABS(vrtDx(i-1,j)), ABS(vrtDx(i,j)) )
482       &              +0.25*(tension( i , j )**2+tension( i ,j+1)**2            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j-1)) )
483       &                +tension(i+1, j )**2+tension(i+1,j+1)**2))            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j))   )
484    
485              grdDiv=MAX( ABS(divDx(i,j)), ABS(divDx(i,j-1)) )
486              grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
487              grdDiv=MAX( grdDiv, ABS(divDy(i-1,j)) )
488    
489              viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
490              viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
491              viscAh_ZLthD(i,j)=(leithD2fac*grdDiv)*L3
492              viscA4_ZLthD(i,j)=(leithD4fac*grdDiv)*L5
493           ELSE           ELSE
494            Asmag2=0d0            viscAh_ZLth(i,j)=0. _d 0
495              viscA4_ZLth(i,j)=0. _d 0
496              viscAh_ZLthD(i,j)=0. _d 0
497              viscA4_ZLthD(i,j)=0. _d 0
498           ENDIF           ENDIF
499    
500           IF (smag4fac.NE.0.) THEN           IF (calcSmag) THEN
501            Asmag4=smag4fac*L4            viscAh_ZSmg(i,j)=L2
502       &             *sqrt(strain(i,j)**2       &      *SQRT(strain(i,j)**2
503       &             +0.25*(tension( i , j )**2+tension( i ,j+1)**2       &        +0.25 _d 0*(tension( i , j )**2+tension( i ,j-1)**2
504       &                +tension(i+1, j )**2+tension(i+1,j+1)**2))       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))
505           ELSE            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
506            Asmag4=0d0            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
507           ENDIF           ENDIF
508    #endif /* AUTODIFF_DISABLE_LEITH */
509    
510  C  Harmonic on Zeta points  C  Harmonic on Zeta points
511           Alin=viscAhZ+vg2*L2+Alth2+Asmag2           Alin=viscAhZ+viscAhGrid*L2rdt
512           viscAh_Z(i,j)=min(viscAhMax,Alin)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
513           IF (useAnisotropicViscAGridMax) THEN  #ifdef ALLOW_3D_VISCAH
514              AlinMax=viscAhGridMax*L2rdt       &          +viscAhZfld(i,j,k,bi,bj)
515              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))  #endif
516           ELSE           viscAh_ZMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
517           IF (vg2Max.GT.0.) THEN           viscAh_Z(i,j)=MAX(viscAh_ZMin(i,j),Alin)
518              AlinMax=vg2Max*L2           viscAh_ZMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
519              viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))           viscAh_Z(i,j)=MIN(viscAh_ZMax(i,j),viscAh_Z(i,j))
520           ENDIF  
521           ENDIF  C  BiHarmonic on Zeta points
522           AlinMin=vg2Min*L2           Alin=viscA4Z+viscA4Grid*L4rdt
523           viscAh_Z(i,j)=max(AlinMin,viscAh_Z(i,j))       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
524    #ifdef ALLOW_3D_VISCA4
525  C  BiHarmonic on Zeta Points       &          +viscA4Zfld(i,j,k,bi,bj)
526           Alin=viscA4Z+vg4*L4+Alth4+Asmag4  #endif
527           viscA4_Z(i,j)=min(viscA4Max,Alin)           viscA4_ZMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
528           IF (useAnisotropicViscAGridMax) THEN           viscA4_Z(i,j)=MAX(viscA4_ZMin(i,j),Alin)
529              AlinMax=viscA4GridMax*L4rdt           viscA4_ZMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
530              viscA4_Z(i,j)=min(AlinMax,viscA4_Z(i,j))           viscA4_Z(i,j)=MIN(viscA4_ZMax(i,j),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))  
         ENDDO  
        ENDDO  
       ELSE  
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          viscAh_D(i,j)=viscAhD  
          viscAh_Z(i,j)=viscAhZ  
          viscA4_D(i,j)=viscA4D  
          viscA4_Z(i,j)=viscA4Z  
531          ENDDO          ENDDO
532         ENDDO         ENDDO
533        ENDIF  
534    #ifdef ALLOW_NONHYDROSTATIC
535           IF ( nonHydrostatic ) THEN
536    C--   Pass Viscosities to calc_gw (if constant, not necessary)
537    
538            IF ( k.LT.Nr ) THEN
539    C     Prepare for next level (next call)
540             DO j=1-OLy,sNy+OLy
541              DO i=1-OLx,sNx+OLx
542                viscAh_W(i,j,k+1,bi,bj) = halfRL*viscAh_D(i,j)
543                viscA4_W(i,j,k+1,bi,bj) = halfRL*viscA4_D(i,j)
544              ENDDO
545             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
584    c     ENDIF
585    
586  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
587        IF (useDiagnostics) THEN        IF (useDiagnostics) THEN
# Line 337  C  BiHarmonic on Zeta Points Line 589  C  BiHarmonic on Zeta Points
589         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
590         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
591         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
592    
593           CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
594           CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)
595           CALL DIAGNOSTICS_FILL(viscAh_ZMax,'VAHZMAX ',k,1,2,bi,bj,myThid)
596           CALL DIAGNOSTICS_FILL(viscA4_ZMax,'VA4ZMAX ',k,1,2,bi,bj,myThid)
597    
598           CALL DIAGNOSTICS_FILL(viscAh_DMin,'VAHDMIN ',k,1,2,bi,bj,myThid)
599           CALL DIAGNOSTICS_FILL(viscA4_DMin,'VA4DMIN ',k,1,2,bi,bj,myThid)
600           CALL DIAGNOSTICS_FILL(viscAh_ZMin,'VAHZMIN ',k,1,2,bi,bj,myThid)
601           CALL DIAGNOSTICS_FILL(viscA4_ZMin,'VA4ZMIN ',k,1,2,bi,bj,myThid)
602    
603           CALL DIAGNOSTICS_FILL(viscAh_DLth,'VAHDLTH ',k,1,2,bi,bj,myThid)
604           CALL DIAGNOSTICS_FILL(viscA4_DLth,'VA4DLTH ',k,1,2,bi,bj,myThid)
605           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)
607    
608           CALL DIAGNOSTICS_FILL(viscAh_DLthD,'VAHDLTHD',
609         &                       k,1,2,bi,bj,myThid)
610           CALL DIAGNOSTICS_FILL(viscA4_DLthD,'VA4DLTHD',
611         &                       k,1,2,bi,bj,myThid)
612           CALL DIAGNOSTICS_FILL(viscAh_ZLthD,'VAHZLTHD',
613         &                       k,1,2,bi,bj,myThid)
614           CALL DIAGNOSTICS_FILL(viscA4_ZLthD,'VA4ZLTHD',
615         &                       k,1,2,bi,bj,myThid)
616    
617           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)
619           CALL DIAGNOSTICS_FILL(viscAh_ZSmg,'VAHZSMAG',k,1,2,bi,bj,myThid)
620           CALL DIAGNOSTICS_FILL(viscA4_ZSmg,'VA4ZSMAG',k,1,2,bi,bj,myThid)
621        ENDIF        ENDIF
622  #endif  #endif
623    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.48
  ViewVC Help
Powered by ViewVC 1.1.22