Drying & flooding

[Ulf Gräwe]

Hans and me just realised that severe problems can occur in GETM
simulations when min_depth is not substantially smaller than
crit_depth.

If, e.g., min_depth is 35 cm and crit_depth is 50 cm, then GETM switches on the drying algorithm, no matter how
deep it is.
The problem starts at that line:
d2i = _ONE_/(crit_depth-2*min_depth)
Since d2i can become negative, this will mess up all the following computations. To avoid this, crit_depth must be more than twice as large as min_depth.

Our suggestion is to add to domain.F90 the following lines after
read(NAMLST,domain):

if (crit_depth .lt. 2.5*min_depth) then
stop 'crit_depth must be larger than 2.5 time min_depth'
end if

However, this immediately stops the simulations. Or shall we only output a warning and set
min_depth=crit_depth/2.5

Hans suggested that as a rule of thumb: min_depth=crit_depth/4

What do you think?

--
Ulf Gr�we

Leibniz Institute for Baltic Sea Research Warnemuende
Dept. for Physical Oceanography and Instrumentation
Seestrasse 15
D-18119 Rostock-Warnemuende
Germany

Pho. : +49 (381) 5197-358
Fax : +49 (381) 5197-440
Web : http://www.io-warnemuende.de/ulf-graewe-en.html
Room : 216

[Bjarne Büchmann]

Hi,

> Hans and me just realised that severe problems can occur in
> GETM simulations when min_depth is not substantially smaller
> than crit_depth.

<snip>

Very interesting points...

> Our suggestion is to add to domain.F90 the following lines after
> read(NAMLST,domain):
>
> if (crit_depth .lt. 2.5*min_depth) then
> stop 'crit_depth must be larger than 2.5 time min_depth'
> end if

I think that is a decent solution. I would rather have it stop
than modify my input.

> However, this immediately stops the simulations.

Yes. Well and sound (in my oppinion).

> Or shall we only output a warning and set
> min_depth=crit_depth/2.5

Do everybody read the outputs from parallized GETM?
Stop would force people to understand.

> Hans suggested that as a rule of thumb: min_depth=crit_depth/4

/Bjarne

MED VENLIG HILSEN

Bjarne Büchmann
Civilingeniør, PhD

Overgaden oven vandet 62 B Postboks 1919
DK - 1023 København K
TLF. +45 3268 697
B...@FRV.DK

[Adolf STIPS]

Dear Ulf and Hans,

did you compare 2 test cases about what happens, consequences?

I'm clearly in favor of the "stop" mechanism, users should be forced
to think about what they are doing, not just using GETM as a black box.

Regards, Adolf

On Feb 23, 2012, at 8:08 AM, Ulf Gräwe wrote:

>
> Hans and me just realised that severe problems can occur in GETM
> simulations when min_depth is not substantially smaller than
> crit_depth.
>
> If, e.g., min_depth is 35 cm and crit_depth is 50 cm, then GETM switches on the drying algorithm, no matter how
> deep it is.
> The problem starts at that line:
> d2i = _ONE_/(crit_depth-2*min_depth)
> Since d2i can become negative, this will mess up all the following computations. To avoid this, crit_depth must be more than twice as large as min_depth.
>
> Our suggestion is to add to domain.F90 the following lines after
> read(NAMLST,domain):
>
> if (crit_depth .lt. 2.5*min_depth) then
> stop 'crit_depth must be larger than 2.5 time min_depth'
> end if
>
> However, this immediately stops the simulations. Or shall we only output a warning and set
> min_depth=crit_depth/2.5
>
> Hans suggested that as a rule of thumb: min_depth=crit_depth/4
>
> What do you think?
>
>
> --

> Ulf Gräwe

[Hans Burchard]

Thanks for the comments, it is committed as suggested now.

Cheers, Hans.

> Hi,
>
>> Hans and me just realised that severe problems can occur in
>> GETM simulations when min_depth is not substantially smaller
>> than crit_depth.
>
> <snip>
>
> Very interesting points...
>
>> Our suggestion is to add to domain.F90 the following lines after
>> read(NAMLST,domain):
>>
>> if (crit_depth .lt. 2.5*min_depth) then
>> stop 'crit_depth must be larger than 2.5 time min_depth'
>> end if
>
> I think that is a decent solution. I would rather have it stop
> than modify my input.
>
>> However, this immediately stops the simulations.
>
> Yes. Well and sound (in my oppinion).
>
>> Or shall we only output a warning and set
>> min_depth=crit_depth/2.5
>
> Do everybody read the outputs from parallized GETM?
> Stop would force people to understand.
>
>> Hans suggested that as a rule of thumb: min_depth=crit_depth/4
>
> /Bjarne
>
> MED VENLIG HILSEN
>

> Bjarne B�chmann
> Civilingeni�r, PhD

>
> Overgaden oven vandet 62 B Postboks 1919

> DK - 1023 K�benhavn K

> TLF. +45 3268 697
> B...@FRV.DK
>

--
Hans Burchard

Leibniz Institute for Baltic Sea Research Warnemuende
Dept. for Physical Oceanography and Instrumentation
Seestrasse 15
D-18119 Rostock-Warnemuende
Germany

Tel. : +49 - 381 - 5197 -140
Fax : +49 - 381 - 5197 -114
web : http://www.io-warnemuende.de/hans-burchard-en.html
E-Mail: hans.b...@io-warnemuende.de
skype : hans.burchard

[Knut]

Hi,

> Hans and me just realised that severe problems can occur in GETM
> simulations when min_depth is not substantially smaller than
> crit_depth.
>
>

> The problem starts at that line:

> d2i = _ONE_/(crit_depth-2*min_depth)
> Since d2i can become negative, this will mess up all the following computations. To avoid this, crit_depth must be more than twice as large as min_depth.
>
> Our suggestion is to add to domain.F90 the following lines after
> read(NAMLST,domain):
>
> if (crit_depth .lt. 2.5*min_depth) then
> stop 'crit_depth must be larger than 2.5 time min_depth'
> end if
>
>

Related to this missing check is the issue asked by me in April last year:

"""

To learn what these quantities are, I had a look into depth_update and
wondered about the different handling for center and velocity points:

d1 = 2*min_depth
d2i = _ONE_/(*crit_depth*-2*min_depth)
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
if (az(i,j) .gt. 0) then
dry_z(i,j)=max(_ZERO_,min(_ONE_,(D(i,j)-_HALF_*d1)*d2i))
end if
if (au(i,j) .gt. 0) then
dry_u(i,j) = max(_ZERO_,min(_ONE_,(DU(i,j)-d1)*d2i))

This means:
if (D .le. min_depth) dry_z = 0 => OK!
if (D .ge. *crit_depth* - min_depth) dry_z = 1 => ???
if (DU .le. 2*min_depth) dry_u = 0 => ???
if (DU .ge. *crit_depth*) dry_u = 1 => OK!

So why not:

dry_?(i,j) = max(_ZERO_,min(_ONE_,(D?(i,j) - d1)*d2i))
with d1 = min_depth
and d2i = _ONE_/(*crit_depth* - min_depth)

This would cause:
if (D? .le. min_depth) dry_? = 0 => OK!
if (D? .ge. *crit_depth*) dry_? = 1 => OK!

...in agreement with the documentation

"""

Of course the check whether crit_depth and min_depth will guarantee a
positive d2i needs to be adapted.

Apart from the different handling of center and velocity points,
especially this result seems to be very strange to me:
if (D .ge. *crit_depth* - min_depth) dry_z = 1

So far only Richard replied to my question:
http://groups.google.com/group/getm-devel/browse_thread/thread/d2fc0ed95885da59/8c40a947b7554791?#8c40a947b7554791

Regards,
Knut

[Johan Van Der Molen (Cefas)]

Hi,

You would have to have info on why it was implemented the way it was
implemented (I don't know).

Knowing what you're doing and why may be essential, because generally,
there can be (numerical) stability considerations or wave-propagation
speed considerations (in particular for flooding) for treating points
differently, or having different rules/thresholds for flooding as
opposed to drying. I don't know about many publications dealing with
these issues specifically, but it's probably worth having a look to see
if there are any. I know that G.S. Stelling published some work on
numerical techniques for over-land flows and flooding of land
reclamation areas in the Netherlands in the early 2000's that is worth
having a look at: see the attached paper, but there may be subsequent
more applied publications as well.

Hope this helps,

Johan

> Regards,
> Knut

This email and any attachments are intended for the named recipient only. Its unauthorised use, distribution, disclosure, storage or copying is not permitted. If you have received it in error, please destroy all copies and notify the sender. In messages of a non-business nature, the views and opinions expressed are the author's own and do not necessarily reflect those of Cefas. Communications on Cefas’ computer systems may be monitored and/or recorded to secure the effective operation of the system and for other lawful purposes.

[Hans Burchard]

<snip>

> You would have to have info on why it was implemented the way it was
> implemented (I don't know).
>
> Knowing what you're doing and why may be essential, because generally,
> there can be (numerical) stability considerations or wave-propagation
> speed considerations (in particular for flooding) for treating points
> differently, or having different rules/thresholds for flooding as
> opposed to drying. I don't know about many publications dealing with
> these issues specifically, but it's probably worth having a look to see
> if there are any. I know that G.S. Stelling published some work on
> numerical techniques for over-land flows and flooding of land
> reclamation areas in the Netherlands in the early 2000's that is worth
> having a look at: see the attached paper, but there may be subsequent
> more applied publications as well.
>

Hi Johan,

thanks a lot for the papr. I assume that Stelling's method works well for
vertically integrated flow, but for fully baroclinic flow things could be
much more difficult.

As far as I remember, Peter Holtermann once implemented one of Stellings
methods (I don't know if it was this one) and it worked fine for 2D flow
and crashed for 3D flow.

Therefore, for the moment, I think we leave it as it is, unless somebody
has a much better idea and implements it.

All the best,

Hans.

> Hope this helps,
>
> Johan
>
>> Regards,
>> Knut
> This email and any attachments are intended for the named recipient only.
> Its unauthorised use, distribution, disclosure, storage or copying is not
> permitted.
> If you have received it in error, please destroy all copies and notify the
> sender. In messages of a non-business nature, the views and opinions
> expressed are the author's own
> and do not necessarily reflect those of Cefas.
> Communications on Cefas� computer systems may be monitored and/or recorded
> to secure the effective operation of the system and for other lawful
> purposes.

--
Hans Burchard

Leibniz Institute for Baltic Sea Research Warnemuende
Dept. for Physical Oceanography and Instrumentation
Seestrasse 15
D-18119 Rostock-Warnemuende
Germany

Tel. : +49 - 381 - 5197 -140