/[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.24 by mlosch, Thu Nov 23 09:49:36 2006 UTC
# Line 72  C     == Global variables == Line 72  C     == Global variables ==
72  #include "GRID.h"  #include "GRID.h"
73  #include "EEPARAMS.h"  #include "EEPARAMS.h"
74  #include "PARAMS.h"  #include "PARAMS.h"
75    #ifdef ALLOW_NONHYDROSTATIC
76    #include "NH_VARS.h"
77    #endif
78    
79  C     == Routine arguments ==  C     == Routine arguments ==
80        INTEGER bi,bj,k        INTEGER bi,bj,k
# Line 90  C     == Routine arguments == Line 93  C     == Routine arguments ==
93    
94  C     == Local variables ==  C     == Local variables ==
95        INTEGER I,J        INTEGER I,J
96          INTEGER kp1
97        _RL smag2fac, smag4fac        _RL smag2fac, smag4fac
98        _RL leith2fac, leith4fac        _RL leith2fac, leith4fac
99        _RL leithD2fac, leithD4fac        _RL leithD2fac, leithD4fac
# Line 99  C     == Local variables == Line 103  C     == Local variables ==
103        _RL Uscl,U4scl        _RL Uscl,U4scl
104        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106          _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107          _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110        _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 210  C     == Local variables ==
210          leithD4fac=0. _d 0          leithD4fac=0. _d 0
211        ENDIF        ENDIF
212    
213    #ifdef ALLOW_AUTODIFF_TAMC
214          IF ( calcLeith .OR. calcSmag ) THEN
215           STOP 'calcLeith or calcSmag not implemented for ADJOINT'
216          ENDIF
217    #endif
218          DO j=1-Oly,sNy+Oly
219            DO i=1-Olx,sNx+Olx
220              viscAh_D(i,j)=viscAhD
221              viscAh_Z(i,j)=viscAhZ
222              viscA4_D(i,j)=viscA4D
223              viscA4_Z(i,j)=viscA4Z
224    c
225              visca4_zsmg(i,j) = 0. _d 0
226              viscah_zsmg(i,j) = 0. _d 0
227    c
228              viscAh_Dlth(i,j) = 0. _d 0
229              viscA4_Dlth(i,j) = 0. _d 0
230              viscAh_DlthD(i,j)= 0. _d 0
231              viscA4_DlthD(i,j)= 0. _d 0
232    c
233              viscAh_DSmg(i,j) = 0. _d 0
234              viscA4_DSmg(i,j) = 0. _d 0
235    c
236              viscAh_ZLth(i,j) = 0. _d 0
237              viscA4_ZLth(i,j) = 0. _d 0
238              viscAh_ZLthD(i,j)= 0. _d 0
239              viscA4_ZLthD(i,j)= 0. _d 0
240            ENDDO
241          ENDDO
242    
243  C     - Viscosity  C     - Viscosity
244        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
245    
246  C      horizontal gradient of horizontal divergence:  C-     Initialise to zero gradient of vorticity & divergence:
247         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
248          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
249            divDx(i,j) = 0.            divDx(i,j) = 0.
250            divDy(i,j) = 0.            divDy(i,j) = 0.
251              vrtDx(i,j) = 0.
252              vrtDy(i,j) = 0.
253          ENDDO          ENDDO
254         ENDDO         ENDDO
255    
256         IF (calcleith) THEN         IF (calcleith) THEN
257    C      horizontal gradient of horizontal divergence:
258    
259  C-       gradient in x direction:  C-       gradient in x direction:
260  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
261           IF (useCubedSphereExchange) THEN           IF (useCubedSphereExchange) THEN
# Line 240  C        to compute d/dy(hDiv), fill cor Line 281  C        to compute d/dy(hDiv), fill cor
281              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)
282            ENDDO            ENDDO
283           ENDDO           ENDDO
284    
285    C      horizontal gradient of vertical vorticity:
286    C-       gradient in x direction:
287             DO j=2-Oly,sNy+Oly
288              DO i=2-Olx,sNx+Olx-1
289                vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
290         &                  *recip_DXG(i,j,bi,bj)
291         &                  *maskS(i,j,k,bi,bj)
292              ENDDO
293             ENDDO
294    C-       gradient in y direction:
295             DO j=2-Oly,sNy+Oly-1
296              DO i=2-Olx,sNx+Olx
297                vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
298         &                  *recip_DYG(i,j,bi,bj)
299         &                  *maskW(i,j,k,bi,bj)
300              ENDDO
301             ENDDO
302    
303         ENDIF         ENDIF
304    
305         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
# Line 279  C Velocity Reynolds Scale Line 339  C Velocity Reynolds Scale
339             U4scl=0.             U4scl=0.
340           ENDIF           ENDIF
341    
342    #ifndef ALLOW_AUTODIFF_TAMC
343           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
344  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
345            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
346       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
347       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
348       &     +((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)  
349    
350  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
351  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 365  C Using it in Leith serves to damp insta
365           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
366  C but this approximation will work on cube  C but this approximation will work on cube
367  c (and differs by as much as 4X)  c (and differs by as much as 4X)
368            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)) )
369            grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i+1,j)) )
370       &     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)))  
371    
372    c This approximation is good to the same order as above...
373            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )
374            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )
375            grdDiv=max( grdDiv, abs(divDy(i,j))   )            grdDiv=max( grdDiv, abs(divDy(i,j))   )
376    
 c This approximation is good to the same order as above...  
377            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
378            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
379            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 396  c This approximation is good to the same
396            viscAh_DSmg(i,j)=0. _d 0            viscAh_DSmg(i,j)=0. _d 0
397            viscA4_DSmg(i,j)=0. _d 0            viscA4_DSmg(i,j)=0. _d 0
398           ENDIF           ENDIF
399    #endif /* ALLOW_AUTODIFF_TAMC */
400    
401  C  Harmonic on Div.u points  C  Harmonic on Div.u points
402           Alin=viscAhD+viscAhGrid*L2rdt           Alin=viscAhD+viscAhGrid*L2rdt
# Line 359  C  BiHarmonic on Div.u points Line 414  C  BiHarmonic on Div.u points
414           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
415           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))
416    
417    #ifdef ALLOW_NONHYDROSTATIC
418    C     /* Pass Viscosities to calc_gw, if constant, not necessary */
419    
420             kp1  = MIN(k+1,Nr)
421    
422             if (k .eq. 1) then
423              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
424              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
425    
426              viscAh_W(i,j,k,bi,bj)=viscAh_D(i,j) /* These values dont get used */
427              viscA4_W(i,j,k,bi,bj)=viscA4_D(i,j)
428             else
429    C Note that previous call of this function has already added half.
430              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
431              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
432    
433              viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)
434              viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)
435             endif
436    #endif /* ALLOW_NONHYDROSTATIC */
437    
438  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
439  C These are (powers of) length scales  C These are (powers of) length scales
440           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
# Line 396  C Velocity Reynolds Scale (Pb here at CS Line 472  C Velocity Reynolds Scale (Pb here at CS
472             U4scl=0.             U4scl=0.
473           ENDIF           ENDIF
474    
475    #ifndef ALLOW_AUTODIFF_TAMC
476  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
477           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
478            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
479       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
480       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
481       &     +((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)  
482    
483  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
484            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 497  C This is the vector magnitude of grad(d
497    
498           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
499  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
500            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)) )
501            grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i,j-1)) )
502       &     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)))  
503    
504            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )
505            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 524  C but this approximation will work on cu
524            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
525            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
526           ENDIF           ENDIF
527    #endif /* ALLOW_AUTODIFF_TAMC */
528    
529  C  Harmonic on Zeta points  C  Harmonic on Zeta points
530           Alin=viscAhZ+viscAhGrid*L2rdt           Alin=viscAhZ+viscAhGrid*L2rdt
# Line 487  C  BiHarmonic on Zeta points Line 560  C  BiHarmonic on Zeta points
560         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
561         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
562         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
563    #ifdef ALLOW_NONHYDROSTATIC
564           CALL DIAGNOSTICS_FILL(viscAh_W,'VISCAHW ',k,1,2,bi,bj,myThid)
565           CALL DIAGNOSTICS_FILL(viscA4_W,'VISCA4W ',k,1,2,bi,bj,myThid)
566    #endif
567    
568         CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
569         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.24
  ViewVC Help
Powered by ViewVC 1.1.22