open TELEMAC-MASCARET The mathematically superior suite of solvers

Hi all,

I believe that there may be an error in the calculation in astro.f which is used to calculate the tide generating force. In my large-scale model (where the tide-generating force is important), I noticed that the model always started to behave differently on 2013-11-30 or 2013-12-01, roughly speaking. I figured out that there is a discontinuity (jump) in the calculation of the tide generating force (FU1,FV1 in marast.f), which is caused by a discontinuity in the moon right ascension AL (calculated in astro.f).

As far as I know, there should not be any discontinuities in FU1, FV1, or in AL (except for jumps by 2*pi). I am trying to solve this problem and I will post a fix as soon as I figure it out, but I was wondering if somebody else has experienced the same problem?

Attached is a fortran code martest.f which calculates FU1,FV1 for the timeseries around 1 December 2013 and prints them to a file, to demonstrate that there is that discontinuity. It should be compiled along with astro.f, jultim.f, dmo.f, atanc.f and tsloc.f from the Telemac source code.

Hello,

Very interesting, I just mentionned recently that this tide generating force was difficult to validate. I you do not find I would be interested in having the whole case and look.

Regards,

Jean-Michel Hervouet

Can you clarify what you mean by your last sentence? What do you mean by 'the whole case'?

I just mean all your files so that I can immediately run it.

JMH

Attached is a zip containing all the files you need. Just compile them and run the resulting executable.

When running the compiled program, 2 text output files will be generated:
astroOut.txt contains the following colums: AT,T,AL,DL,X,NU,O,OM,I0
marOut.txt contains the following colums: AT,AL,AS,FU1,FV1

see the files astro.f and martest.f for a description of the variables.

The lunar positions are calculated with 25 november, 2013 as a start date.
Around 30 November-1 December 2014, something strange happens: there is a jump
in AL (right ascension of the moon). There is also a jump (by pi/2) in the angle
nu (see the book Hydrodynamics of Free Surface Flows p. 314, eq. A.60).

This angle nu is calculated line 206 of astro.f . What I do not understand is why it is calculated with atanc so that nu remains within +/- pi/4 from the angle O.

I don't think I will be able to resolve this problem so any help would be greatly appreciated.

JM,

I believe that the issue can be resolved by replacing the lines
byin the file astro.f.

I compared the computations of astro.f for the right ascension and the declination of the sun and moon with the predictions from pyEphem (rhodesmill.org/pyephem/), a python implementation of the XEphem program. I found good agreement (to within about 1 percent) over the period 1990-2035.

However, this would only mean that the position of the sun and moon are calculated correctly; this does not guarantee that the resulting tidal body forces are calculated correctly.

I agree that this is a difficult problem to validate since the tidal body force is often a second-order effect. One option could be to run the model for an idealized case (a rectangular domain with constant depth and wall boundaries) with only the tidal body force, and compare it to a similar run from another model?

Hello,

Hum, this does not seem to solve the matter here, on my side I find that the problem is function ATANC that gives NU and the discontinuity is in NU. So probably function ATANC which is not a Fortran intrinsic function but some homemade function is guilty, as it is discontinuous even with continuous arguments.

Unfortunately I will be able to look at this only next Wednesday, be patient...

With best regards,

Jean-Michel Hervouet

I think you understood my message incorrectly: the old code contains the ATANC function, the new code does not contain it.
I agree, the ATANC function looks strange, I did not understand the logic behind it.

Hello,

Yes, I missed the fact that you changed ATANC into ATAN in the computation of NU. But did you change it also in the computation of X ?

Regards,

JMH

Hello,

After thinking a while, I propose to replace the use of ATANC by intrinsic Fortran function ATAN2. See enclosed subroutine astro.f.
If you find it is OK I'll put it in the next release, and ATANC may now be removed.

With best regards,

Jean-Michel Hervouet