/[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.22 by heimbach, Wed May 31 19:53:15 2006 UTC
# Line 99  C     == Local variables == Line 99  C     == Local variables ==
99        _RL Uscl,U4scl        _RL Uscl,U4scl
100        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102          _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103          _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106        _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 206  C     == Local variables ==
206          leithD4fac=0. _d 0          leithD4fac=0. _d 0
207        ENDIF        ENDIF
208    
209    #ifdef ALLOW_AUTODIFF_TAMC
210          IF ( calcLeith .OR. calcSmag ) THEN
211           STOP 'calcLeith or calcSmag not implemented for ADJOINT'
212          ENDIF
213          DO j=1-Oly,sNy+Oly
214            DO i=1-Olx,sNx+Olx
215              viscAh_D(i,j)=viscAhD
216              viscAh_Z(i,j)=viscAhZ
217              viscA4_D(i,j)=viscA4D
218              viscA4_Z(i,j)=viscA4Z
219    c
220              visca4_zsmg(i,j) = 0. _d 0
221              viscah_zsmg(i,j) = 0. _d 0
222    c
223              viscAh_Dlth(i,j) = 0. _d 0
224              viscA4_Dlth(i,j) = 0. _d 0
225              viscAh_DlthD(i,j)= 0. _d 0
226              viscA4_DlthD(i,j)= 0. _d 0
227    c
228              viscAh_DSmg(i,j) = 0. _d 0
229              viscA4_DSmg(i,j) = 0. _d 0
230    c
231              viscAh_ZLth(i,j) = 0. _d 0
232              viscA4_ZLth(i,j) = 0. _d 0
233              viscAh_ZLthD(i,j)= 0. _d 0
234              viscA4_ZLthD(i,j)= 0. _d 0
235            ENDDO
236          ENDDO
237    #endif
238    
239    
240    
241  C     - Viscosity  C     - Viscosity
242        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
243    
244  C      horizontal gradient of horizontal divergence:  C-     Initialise to zero gradient of vorticity & divergence:
245         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
246          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
247            divDx(i,j) = 0.            divDx(i,j) = 0.
248            divDy(i,j) = 0.            divDy(i,j) = 0.
249              vrtDx(i,j) = 0.
250              vrtDy(i,j) = 0.
251          ENDDO          ENDDO
252         ENDDO         ENDDO
253    
254         IF (calcleith) THEN         IF (calcleith) THEN
255    C      horizontal gradient of horizontal divergence:
256    
257  C-       gradient in x direction:  C-       gradient in x direction:
258  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
259           IF (useCubedSphereExchange) THEN           IF (useCubedSphereExchange) THEN
# Line 240  C        to compute d/dy(hDiv), fill cor Line 279  C        to compute d/dy(hDiv), fill cor
279              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)
280            ENDDO            ENDDO
281           ENDDO           ENDDO
282    
283    C      horizontal gradient of vertical vorticity:
284    C-       gradient in x direction:
285             DO j=2-Oly,sNy+Oly
286              DO i=2-Olx,sNx+Olx-1
287                vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
288         &                  *recip_DXG(i,j,bi,bj)
289         &                  *maskS(i,j,k,bi,bj)
290              ENDDO
291             ENDDO
292    C-       gradient in y direction:
293             DO j=2-Oly,sNy+Oly-1
294              DO i=2-Olx,sNx+Olx
295                vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
296         &                  *recip_DYG(i,j,bi,bj)
297         &                  *maskW(i,j,k,bi,bj)
298              ENDDO
299             ENDDO
300    
301         ENDIF         ENDIF
302    
303         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
# Line 279  C Velocity Reynolds Scale Line 337  C Velocity Reynolds Scale
337             U4scl=0.             U4scl=0.
338           ENDIF           ENDIF
339    
340    #ifndef ALLOW_AUTODIFF_TAMC
341           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
342  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
343            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
344       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
345       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
346       &     +((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)  
347    
348  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
349  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 363  C Using it in Leith serves to damp insta
363           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
364  C but this approximation will work on cube  C but this approximation will work on cube
365  c (and differs by as much as 4X)  c (and differs by as much as 4X)
366            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)) )
367            grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i+1,j)) )
368       &     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)))  
369    
370    c This approximation is good to the same order as above...
371            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )
372            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )
373            grdDiv=max( grdDiv, abs(divDy(i,j))   )            grdDiv=max( grdDiv, abs(divDy(i,j))   )
374    
 c This approximation is good to the same order as above...  
375            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
376            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
377            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 394  c This approximation is good to the same
394            viscAh_DSmg(i,j)=0. _d 0            viscAh_DSmg(i,j)=0. _d 0
395            viscA4_DSmg(i,j)=0. _d 0            viscA4_DSmg(i,j)=0. _d 0
396           ENDIF           ENDIF
397    #endif /* ALLOW_AUTODIFF_TAMC */
398    
399  C  Harmonic on Div.u points  C  Harmonic on Div.u points
400           Alin=viscAhD+viscAhGrid*L2rdt           Alin=viscAhD+viscAhGrid*L2rdt
# Line 396  C Velocity Reynolds Scale (Pb here at CS Line 449  C Velocity Reynolds Scale (Pb here at CS
449             U4scl=0.             U4scl=0.
450           ENDIF           ENDIF
451    
452    #ifndef ALLOW_AUTODIFF_TAMC
453  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
454           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
455            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
456       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
457       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
458       &     +((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)  
459    
460  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
461            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 474  C This is the vector magnitude of grad(d
474    
475           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
476  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
477            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)) )
478            grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i,j-1)) )
479       &     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)))  
480    
481            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )
482            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 501  C but this approximation will work on cu
501            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
502            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
503           ENDIF           ENDIF
504    #endif /* ALLOW_AUTODIFF_TAMC */
505    
506  C  Harmonic on Zeta points  C  Harmonic on Zeta points
507           Alin=viscAhZ+viscAhGrid*L2rdt           Alin=viscAhZ+viscAhGrid*L2rdt
# Line 481  C  BiHarmonic on Zeta points Line 531  C  BiHarmonic on Zeta points
531         ENDDO         ENDDO
532        ENDIF        ENDIF
533    
534    
535  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
536        IF (useDiagnostics) THEN        IF (useDiagnostics) THEN
537         CALL DIAGNOSTICS_FILL(viscAh_D,'VISCAHD ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_D,'VISCAHD ',k,1,2,bi,bj,myThid)
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22