[MITgcm-devel] mom_calc_visc.F

[David Ferreira]

Again I reply to myself. But this time, I have a conclusion.
There are 2 bugs in mom_calc_visc.F:

1) with useAreaViscLength=.TRUE., the biharmonic viscosities
 computed at the corner and center points differ by a factor 4 (the 
corner ones
are wrong ones)

2) with useAreaViscLength=.FALSE., the computation of L4rdt is wrong at 
both
positions (a 1 dropped in the formula and a 16 turned into a 6).

I'm going to check in fixes for those 2 bugs. Both will affect the runs with
variables viscosities, 1) probably more than 2). If you have a 
complaint, now
is good time (cs510 people ?).
Cheers,
david


David Ferreira wrote:
> It probably means that the coefficient for the computation
> at the vorticity points should be changed similarly (from 1/8 to 1/32).
> We are left with the "1/32 or 1/20" problem depending on 
> useAreaViscLength
> david
>
>
>
> chris hill wrote:
>> Darn - does that mean cnh has to sort it out! I'll take a look.
>>
>> As I said to David (when I did respond to his original question) - 
>> there is no 100% _correct_ answer here. The optimal scaling of 
>> effective dissipation/diffusion with grid spacing does not need to 
>> follow the spacing exactly .e.g on a lat-lon grid typically the 
>> meridional component of biharmonic mixing is scaled with a l^~2.5 not 
>> l^4 times coefficient factor. Nevertheless we should have a 
>> consistent transfer function and be able to say what it is.
>>
>> cnh
>>
>> Jean-Michel Campin wrote:
>>> Hi David,
>>> In the CVS log-message, I found, from cnh (revision 1.25):
>>>> Change from 1/8 to 1/32 scaling from Baylor.
>>>> Need to check if there are any verif expts affected.
>>> May be this could help ?
>>>
>>> Jean-Michel
>>>
>>> On Mon, Mar 09, 2009 at 10:18:34PM -0400, David Ferreira wrote:
>>>> 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
>>>>
>>>> _______________________________________________
>>>> MITgcm-devel mailing list
>>>> MITgcm-devel at mitgcm.org
>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>> _______________________________________________
>>> MITgcm-devel mailing list
>>> MITgcm-devel at mitgcm.org
>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>
>>
>> _______________________________________________
>> MITgcm-devel mailing list
>> MITgcm-devel at mitgcm.org
>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>
>