[MITgcm-devel] lab_sea sensitivity between platforms/compilers

[Martin Losch]

Hi Jean-Michel,

I might be chasing a Friday afternoon blackout, but I ran lab_sea  
without the dynamics solver, and got very similar numbers on for my  
different platforms.
Turning off the thermodynamics did not have the same effect. But it  
may very well be a coincidence for my particular combination of  
platforms/compilers.

How do you normally test these things?

Martin

On Jul 17, 2009, at 11:48 PM, Jean-Michel Campin wrote:

> Hi Martin,
>
> Are you sure that this lab_sea "sensitivity" to compiler
> is not comming from the thermodynamics part ?
> If I look at the results of ifort/ifc (or pgi) test, compared
> with g77/gfortran, I still get a good results for
> global_ocean.cs32x15.icedyn (10 or 11 digits for cg2d)
> but all the lab_sea are failing with only 5 digits.
>
> I don't feel enough courage now to look in details in seaice_growth.F
> but do we still apply instantaneous melting when SST > t_freeze ?
> and do we still start from 1.m ice thickness everywhere
> in all those lab_sea test ?
>
> Cheers,
> Jean-Michel
>
> On Fri, Jul 17, 2009 at 10:00:16PM +0200, Martin Losch wrote:
>> Hi David (and others),
>>
>> your are right, SQRT(deltaC)/(SQRT(deltaC)+SEAICE_EPS) is order one  
>> or
>> zero (if deltaC=0), the correct expression would have been
>>>> deltaCreg = SQRT(deltaC)+SEAICE_EPS
>> to avoid divisions by zero.
>> And consequently, when I do that, then the sensitivity returns.
>> I found that when I fix zetaC to a small value (say 1.), the  
>> sensitivity
>> goes away, but since zetaC is a factor that scales everything, it  
>> just
>> means that the entire stress term is scaled down by orders of  
>> magnitude
>> (normally zetaC is order 1e10 and more), and the dynamics solver  
>> gives
>> nothing much. All that means, is that the dynamics solver is very
>> sensitive; and I still cannot figure out why (the same is true for  
>> LSR of
>> course), rats.
>>
>> Martin
>>
>> On Jul 17, 2009, at 4:32 PM, David Ferreira wrote:
>>
>>> Hi Martin,
>>>> In the EVP code there is a line
>>>>           deltaCreg = SQRT(MAX(deltaC,SEAICE_EPS_SQ))
>>>> If you replace this line with
>>>>           deltaCreg = SQRT(deltaC)/(SQRT(deltaC)+SEAICE_EPS)
>>> It doesn't look like these 2 are equivalent.
>>> Also, what about something like :
>>>
>>> IF deltaC >= SEAICE_EPS_SQ
>>>  deltaCreg =SQRT(deltaC)
>>> else
>>>  deltaCreg = SEAICE_EPS
>>> end
>>>
>>> it adds a "IF" but remove a "MAX" ... Not sure which one costs more
>>> david
>>>
>>>
>>>
>>>> the number of matching digits increase from 5 to 12 (for cg2d) for
>>>> the lab_sea/input experiment, between a
>>>> Linux csysl18 2.6.22.19-0.2-bigsmp #1 SMP 2008-12-18 10:17:03 +0100
>>>> i686 i686 i386
>>>> gfortran (gcc version 4.2.1)
>>>> and a
>>>> Darwin csysm15.dmawi.de 9.7.0 Darwin Kernel Version 9.7.0: Tue Mar
>>>> 31 22:52:17 PDT 2009; gfortran (gcc version 4.3.1)
>>>>
>>>> Both expressions give about the same results for large DeltaC, but
>>>> still I would like to find out what would be the correct  
>>>> replacement
>>>> of the SQRT(MAX(deltaC,SEAICE_EPS_SQ)). Here's what did *NOT* work:
>>>>           deltaCreg = max(sqrt(deltaC),seaice_eps)
>>>>           deltaCreg = sqrt(dmax1(deltaC,seaice_eps_sq))
>>>>           deltaCreg = SQRT(MAX(DBLE(deltaC),DBLE(SEAICE_EPS_SQ)))
>>>>           deltaCreg = SQRT(MAX(SEAICE_EPS_SQ,deltaC))
>>>>
>>>> Any idea, what might be going on? If I cannot solve this problem I
>>>> would suggest replacing this expression with the more stable one,
>>>> also in seaice_calc_viscosities (although it will change all
>>>> results).
>>>>
>>>> Martin
>>>> _______________________________________________
>>>> 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