[MITgcm-devel] mom_calc_visc.F

[Hong Zhang]

Hi David,
We use useAreaViscLength=.FALSE. in CS510 adjoint.
The forward run is OK but we did have blowup at corner in backward 
integration.
Have no idea it's related to this bug.
But we capped the adjoint fields in the neighborhood of corners
to get a stable adjoint run.

cheers,
hong

David Ferreira wrote:
> 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
>>
>>
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel