/[MITgcm]/MITgcm/pkg/mom_common/mom_calc_visc.F
ViewVC logotype

Diff of /MITgcm/pkg/mom_common/mom_calc_visc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.19 by baylor, Mon Oct 10 19:49:48 2005 UTC revision 1.30 by mlosch, Wed Nov 14 15:55:26 2007 UTC
# 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
# Line 72  C     == Global variables == Line 74  C     == Global variables ==
74  #include "GRID.h"  #include "GRID.h"
75  #include "EEPARAMS.h"  #include "EEPARAMS.h"
76  #include "PARAMS.h"  #include "PARAMS.h"
77    #ifdef ALLOW_NONHYDROSTATIC
78    #include "NH_VARS.h"
79    #endif
80    
81  C     == Routine arguments ==  C     == Routine arguments ==
82        INTEGER bi,bj,k        INTEGER bi,bj,k
# Line 90  C     == Routine arguments == Line 95  C     == Routine arguments ==
95    
96  C     == Local variables ==  C     == Local variables ==
97        INTEGER I,J        INTEGER I,J
98    #ifdef ALLOW_NONHYDROSTATIC
99          INTEGER kp1
100    #endif
101        _RL smag2fac, smag4fac        _RL smag2fac, smag4fac
102        _RL leith2fac, leith4fac        _RL leith2fac, leith4fac
103        _RL leithD2fac, leithD4fac        _RL leithD2fac, leithD4fac
# Line 99  C     == Local variables == Line 107  C     == Local variables ==
107        _RL Uscl,U4scl        _RL Uscl,U4scl
108        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110          _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111          _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114        _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 204  C     == Local variables == Line 214  C     == Local variables ==
214          leithD4fac=0. _d 0          leithD4fac=0. _d 0
215        ENDIF        ENDIF
216    
217    #ifdef ALLOW_AUTODIFF_TAMC
218          IF ( calcLeith .OR. calcSmag ) THEN
219           STOP 'calcLeith or calcSmag not implemented for ADJOINT'
220          ENDIF
221    #endif
222          DO j=1-Oly,sNy+Oly
223            DO i=1-Olx,sNx+Olx
224              viscAh_D(i,j)=viscAhD
225              viscAh_Z(i,j)=viscAhZ
226              viscA4_D(i,j)=viscA4D
227              viscA4_Z(i,j)=viscA4Z
228    c
229              visca4_zsmg(i,j) = 0. _d 0
230              viscah_zsmg(i,j) = 0. _d 0
231    c
232              viscAh_Dlth(i,j) = 0. _d 0
233              viscA4_Dlth(i,j) = 0. _d 0
234              viscAh_DlthD(i,j)= 0. _d 0
235              viscA4_DlthD(i,j)= 0. _d 0
236    c
237              viscAh_DSmg(i,j) = 0. _d 0
238              viscA4_DSmg(i,j) = 0. _d 0
239    c
240              viscAh_ZLth(i,j) = 0. _d 0
241              viscA4_ZLth(i,j) = 0. _d 0
242              viscAh_ZLthD(i,j)= 0. _d 0
243              viscA4_ZLthD(i,j)= 0. _d 0
244            ENDDO
245          ENDDO
246    
247  C     - Viscosity  C     - Viscosity
248        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
249    
250  C      horizontal gradient of horizontal divergence:  C-     Initialise to zero gradient of vorticity & divergence:
251         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
252          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
253            divDx(i,j) = 0.            divDx(i,j) = 0.
254            divDy(i,j) = 0.            divDy(i,j) = 0.
255              vrtDx(i,j) = 0.
256              vrtDy(i,j) = 0.
257          ENDDO          ENDDO
258         ENDDO         ENDDO
259    
260         IF (calcleith) THEN         IF (calcleith) THEN
261    C      horizontal gradient of horizontal divergence:
262    
263  C-       gradient in x direction:  C-       gradient in x direction:
264  #ifndef ALLOW_AUTODIFF_TAMC  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
265           IF (useCubedSphereExchange) THEN           IF (useCubedSphereExchange) THEN
266  C        to compute d/dx(hDiv), fill corners with appropriate values:  C        to compute d/dx(hDiv), fill corners with appropriate values:
267             CALL FILL_CS_CORNER_TR_RL( .TRUE., hDiv, bi,bj, myThid )             CALL FILL_CS_CORNER_TR_RL( .TRUE., .FALSE.,
268         &                                hDiv, bi,bj, myThid )
269           ENDIF           ENDIF
270  #endif  cph-exch2#endif
271           DO j=2-Oly,sNy+Oly-1           DO j=2-Oly,sNy+Oly-1
272            DO i=2-Olx,sNx+Olx-1            DO i=2-Olx,sNx+Olx-1
273              divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)              divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)
# Line 229  C        to compute d/dx(hDiv), fill cor Line 275  C        to compute d/dx(hDiv), fill cor
275           ENDDO           ENDDO
276    
277  C-       gradient in y direction:  C-       gradient in y direction:
278  #ifndef ALLOW_AUTODIFF_TAMC  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
279           IF (useCubedSphereExchange) THEN           IF (useCubedSphereExchange) THEN
280  C        to compute d/dy(hDiv), fill corners with appropriate values:  C        to compute d/dy(hDiv), fill corners with appropriate values:
281             CALL FILL_CS_CORNER_TR_RL(.FALSE., hDiv, bi,bj, myThid )             CALL FILL_CS_CORNER_TR_RL(.FALSE., .FALSE.,
282         &                                hDiv, bi,bj, myThid )
283           ENDIF           ENDIF
284  #endif  cph-exch2#endif
285           DO j=2-Oly,sNy+Oly-1           DO j=2-Oly,sNy+Oly-1
286            DO i=2-Olx,sNx+Olx-1            DO i=2-Olx,sNx+Olx-1
287              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)
288            ENDDO            ENDDO
289           ENDDO           ENDDO
290    
291    C      horizontal gradient of vertical vorticity:
292    C-       gradient in x direction:
293             DO j=2-Oly,sNy+Oly
294              DO i=2-Olx,sNx+Olx-1
295                vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
296         &                  *recip_DXG(i,j,bi,bj)
297         &                  *maskS(i,j,k,bi,bj)
298              ENDDO
299             ENDDO
300    C-       gradient in y direction:
301             DO j=2-Oly,sNy+Oly-1
302              DO i=2-Olx,sNx+Olx
303                vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
304         &                  *recip_DYG(i,j,bi,bj)
305         &                  *maskW(i,j,k,bi,bj)
306              ENDDO
307             ENDDO
308    
309         ENDIF         ENDIF
310    
311         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
# Line 259  C These are (powers of) length scales Line 325  C These are (powers of) length scales
325           L2rdt=0.25 _d 0*recip_dt*L2           L2rdt=0.25 _d 0*recip_dt*L2
326    
327           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
328            L4rdt=0.125 _d 0*recip_dt*rA(i,j,bi,bj)**2            L4rdt=0.03125 _d 0*recip_dt*L2**2
329           ELSE           ELSE
330            L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4            L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4
331       &                            +recip_DYF(I,J,bi,bj)**4)       &                            +recip_DYF(I,J,bi,bj)**4)
# Line 279  C Velocity Reynolds Scale Line 345  C Velocity Reynolds Scale
345             U4scl=0.             U4scl=0.
346           ENDIF           ENDIF
347    
348    #ifndef ALLOW_AUTODIFF_TAMC
349           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
350  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
351            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
352       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
353       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
354       &     +((vort3(i+1,j+1)-vort3(i,j+1))       &                        + vrtDy(i,j)*vrtDy(i,j) )  )
      &               *recip_DXG(i,j+1,bi,bj))**2  
      &     +((vort3(i+1,j+1)-vort3(i+1,j))  
      &               *recip_DYG(i+1,j,bi,bj))**2)  
355    
356  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
357  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
# Line 307  C Using it in Leith serves to damp insta Line 371  C Using it in Leith serves to damp insta
371           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
372  C but this approximation will work on cube  C but this approximation will work on cube
373  c (and differs by as much as 4X)  c (and differs by as much as 4X)
374            grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))            grdVrt=max( abs(vrtDx(i,j+1)), abs(vrtDx(i,j)) )
375            grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i+1,j)) )
376       &     abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))            grdVrt=max( grdVrt, abs(vrtDy(i,j))   )
           grdVrt=max(grdVrt,  
      &     abs((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj)))  
           grdVrt=max(grdVrt,  
      &     abs((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj)))  
377    
378    c This approximation is good to the same order as above...
379            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )
380            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )
381            grdDiv=max( grdDiv, abs(divDy(i,j))   )            grdDiv=max( grdDiv, abs(divDy(i,j))   )
382    
 c This approximation is good to the same order as above...  
383            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
384            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
385            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
# Line 342  c This approximation is good to the same Line 402  c This approximation is good to the same
402            viscAh_DSmg(i,j)=0. _d 0            viscAh_DSmg(i,j)=0. _d 0
403            viscA4_DSmg(i,j)=0. _d 0            viscA4_DSmg(i,j)=0. _d 0
404           ENDIF           ENDIF
405    #endif /* ALLOW_AUTODIFF_TAMC */
406    
407  C  Harmonic on Div.u points  C  Harmonic on Div.u points
408           Alin=viscAhD+viscAhGrid*L2rdt           Alin=viscAhD+viscAhGrid*L2rdt
# Line 359  C  BiHarmonic on Div.u points Line 420  C  BiHarmonic on Div.u points
420           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
421           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))
422    
423    #ifdef ALLOW_NONHYDROSTATIC
424    C     /* Pass Viscosities to calc_gw, if constant, not necessary */
425    
426             kp1  = MIN(k+1,Nr)
427    
428             if (k .eq. 1) then
429              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
430              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
431    
432    C     /* These values dont get used */
433              viscAh_W(i,j,k,bi,bj)=viscAh_D(i,j)
434              viscA4_W(i,j,k,bi,bj)=viscA4_D(i,j)
435             else
436    C Note that previous call of this function has already added half.
437              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
438              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
439    
440              viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)
441              viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)
442             endif
443    #endif /* ALLOW_NONHYDROSTATIC */
444    
445  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
446  C These are (powers of) length scales  C These are (powers of) length scales
447           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
# Line 396  C Velocity Reynolds Scale (Pb here at CS Line 479  C Velocity Reynolds Scale (Pb here at CS
479             U4scl=0.             U4scl=0.
480           ENDIF           ENDIF
481    
482    #ifndef ALLOW_AUTODIFF_TAMC
483  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
484           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
485            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
486       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
487       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
488       &     +((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)  
489    
490  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
491            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
# Line 421  C This is the vector magnitude of grad(d Line 504  C This is the vector magnitude of grad(d
504    
505           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
506  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
507            grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))            grdVrt=max( abs(vrtDx(i-1,j)), abs(vrtDx(i,j)) )
508            grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i,j-1)) )
509       &     abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))            grdVrt=max( grdVrt, abs(vrtDy(i,j))   )
           grdVrt=max(grdVrt,  
      &     abs((vort3(i-1,j)-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)))  
510    
511            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )
512            grdDiv=max( grdDiv, abs(divDy(i,j))   )            grdDiv=max( grdDiv, abs(divDy(i,j))   )
# Line 452  C but this approximation will work on cu Line 531  C but this approximation will work on cu
531            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
532            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
533           ENDIF           ENDIF
534    #endif /* ALLOW_AUTODIFF_TAMC */
535    
536  C  Harmonic on Zeta points  C  Harmonic on Zeta points
537           Alin=viscAhZ+viscAhGrid*L2rdt           Alin=viscAhZ+viscAhGrid*L2rdt
# Line 487  C  BiHarmonic on Zeta points Line 567  C  BiHarmonic on Zeta points
567         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
568         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
569         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
570    #ifdef ALLOW_NONHYDROSTATIC
571           CALL DIAGNOSTICS_FILL(viscAh_W,'VISCAHW ',k,1,2,bi,bj,myThid)
572           CALL DIAGNOSTICS_FILL(viscA4_W,'VISCA4W ',k,1,2,bi,bj,myThid)
573    #endif
574    
575         CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
576         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.19  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22