[FLASH-BUGS] "minor" glitch in flatten.F90 (part 2)

[Tomasz Plewa]

Markus -

I think any bias-free formula is OK as long as the threshold for
pressure is reasonable. This is, of course, problem dependent.

I am generally against using floor-values in the code as they are (1)
problem-dependent, and (2) mask some real problems. One can almost
always find a reasonable (unbiased, right order of magnitued) estimate
for the limiting value based on physical or numerical
arguments. 1.e-10 or 1.e+10 does not really mean much when taken out
of context. Hardwiring -10/+10 and -2/+2, as you suggested, may not be
the best solution. Tying those constant to epsiln may work better.

Best wishes -

Tomek
--
On Wed, Jul 16, 2003 at 01:33:44PM +0100, Markus Gross wrote:
> Hi again!
> 
> Realizing that my fix was not much better than the original for any non 
> constructed case, i.e. when rounding error may be present, may I suggest the 
> following:
> 
> !this is the original:
> !!$       if (dp(i) .LT. 0.e0) then
> !!$          scrch2(i) = scrch3(i+1)
> !!$       else
> !!$          scrch2(i) = scrch3(i-1)
> !!$       endif
> 
> !my "new" suggestion:
>        if (dp(i) .LT. -10.e0) THEN
>           scrch2(i) = scrch3(i+1)
>        ELSE
>           if (dp(i) .GT. 10.e0) THEN
>              scrch2(i) = scrch3(i-1)
>           ELSE
>              weight1 = 0.5d0*dabs(derf(dp(i)+2.d0)-1.d0)
>              weight2 = 0.5d0*dabs(derf(dp(i)-2.d0)+1.d0)
>              scrch2(i) = scrch3(i+1) * weight1 +&
>                          scrch3(i-1) * weight2 +&
>                          scrch3(i  ) * (1.d0-(weight1+weight2))
>           END IF
>        END IF
> 
> This of course assumes that a -2.d0< dp(i) < 2.d0 is effectively a dp(i)=0 
> which for my case is true. Modifications for other "small" values is of course 
> straight forward.
> 
> Again, what do you think? A weibull may of course be faster than the erf...
> 
> Regards,
> 
> Markus.
> 
> 

--