[MITgcm-devel] lab_sea sensitivity between platforms/compilers

[Martin Losch]

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