[MITgcm-devel] mom_calc_visc.F

[David Ferreira]

Since nobody answers my question, I'm doing my own
follow-up: Even more bugs in mom_calc_visc.F ?

My first below concerned the computation of the viscosities
at the div points. There is an equivalent for at vorticity points:

CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
C These are (powers of) length scales
         IF (useAreaViscLength) THEN
          L2=rAz(i,j,bi,bj)
          L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2
         ELSE
          L2=2. _d 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))
          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

As for the divergence, the 2 methods for a carsetian square grid
converge for L2 but not L4rdt. Even weirder, the coefficients
are now 1/8 and 1/20 (compared to 1/32 and 1/20 at the div points).

In the best case scenario, that's only 2 correct coeff out of 4, no ?
Comments ?
david



David Ferreira wrote:
> Hi,
> I've got a question about the variable bi-harmonic viscosities. In the
> MITgcm, there is the possibility to specify a non-dimensional horizontal
> bi-harmonic viscosity through viscA4Grid. The model then computes
> the actual viscosity taking into account the local size of the 
> grid-point.
> This local grid-scale can be computed 2 different ways (through the
> flag useAreaViscLength):
>
> C These are (powers of) length scales
>         IF (useAreaViscLength) THEN
>          L2=rA(i,j,bi,bj)
>          L4rdt=0.03125 _d 0*recip_dt*L2**2
>         ELSE
>          L2=2. _d 0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))
>          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
>
> The default (useAreaViscLength=.FALSE.) compute the grid
> factor directly from the length scales, otherwise it is done from
> the area of the grid-point (for the cube-sphere for example).
>
> I was expecting that the 2 formula above would give the same results
> for  cartesian square grid-box (dx=dy). This is the case for L2 but 
> not for
> L4rdt.
> If useAreaViscLength is TRUE, then L4rdt = dx^4/(32*deltaT),
> otherwise it is L4rdt = dx^4/(20*deltaT).
>
> Does anyone know which coefficient is the correct one ?
> Thanks,
> david
>
>
>
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel