PR: Fix vars2fortran: No typespec for argument when compiling HELP3O

[jnsebgosselin]

Fixes #24

When compiling the Fortran code for HELP, we get a lot of warnings of the type :

{}
In: :HELP3O:pyhelp/HELP3O.FOR:ahu
vars2fortran: No typespec for argument "m".

This means that the typespec of the function variables need to be defined explicitly in the code. The goal of this PR is to fix that by declaring explicitly the variables in each function that raises a warning. This need to be done very carefully and rigorously though not to introduce unwanted errors in the computation.

Below is a part of the log from AppVeyor during the building process:

[jnsebgosselin]

In function COMPUT, if I use Double Precision instead of Real to declare the variables typespec, then the tests that are run against the RCA test case fail. The RCA test case is an example that was provided with the original Fortran source code of HELP in 1997. We use the results from the RCA case as a reference to test the model each time modifications are done to the code. The input and output files for the RCA test case are available here.

C  *****************************  COMPUT  *****************************

C  COMPUT is used in function AHU, which accumulates heat units and 
C  radiation to calculate potential heat units for each crop before daily
C  simulation begin.

       FUNCTION COMPUT(AC,A,B,I,N)
       Double Precision, INTENT(IN) :: AC
       Double Precision, INTENT(IN) :: A
       Double Precision, INTENT(IN) :: B
       INTEGER, INTENT(IN) :: I
       INTEGER, INTENT(IN) :: N

       AI       =DFLOAT(I)-0.5
       AN       =DFLOAT(N)
       ANG      =6.283185*AI/AN
       COMPUT   =AC+A*COS(ANG)+B*SIN(ANG)
       RETURN
       END

The tests succeed if I declare the variables as REAL in COMPUT.

Below are presented the results from the reference RCA test case, the results obtained with PyHELP when declaring the variables as REAL and those obtained with PyHELP when declaring the variable as DOUBLE PRECISION

ref_evapo_RCA_case : [32.178, 33.420, 36.540] in/yea
evapo_with_real : [32.1836, 33.416042, 36.53447 ] in/year
evapo_with_dreal : [34.042145 35.85652 38.433487] in/year

As we can see, using REAL or DOUBLE PRECISION has a significant impact on the results. Note that the results are in in/year.

[jnsebgosselin]

As explained in the StackOverflow answer (https://stackoverflow.com/a/50918289/4481445):

the double precision definition can change across different platforms and compilers, it is a good practice to explicitly declare the desired kind of the constant using the standard Fortran convention

So here is what I think is probably happening. The RCA test case was run with an executable compiled on a 16bits Windows machine (but I'm not sure this is relevant). They probably compiled with some flags, so that the implicit REAL had a 32bits precision, or maybe the implicit precision of REAL was 32bits already by default on their system.

On my 64bits Windows machine, when compiling without any flags, the implicit precision of Real is 32bits, while the precision of DOUBLE PRECISION is 64bits. I can confirm this is the case, because if I set the kind of REAL explicitly to 32bits and 64bits in COMPUT, I get the exact same results then those obtained when using, respectively, REAL and DOUBLE PRECISION :

C  *****************************  COMPUT  *****************************

       FUNCTION COMPUT(AC,A,B,I,N)
       integer, parameter :: sp = selected_real_kind(6, 37)
       integer, parameter :: dp = selected_real_kind(15, 307)
       integer, parameter :: qp = selected_real_kind(33, 4931) 
       real(kind=qp) :: AC
       real(kind=qp) :: A
       real(kind=qp) :: B
       INTEGER, INTENT(IN) :: I
       INTEGER, INTENT(IN) :: N

       AI       =REAL(I, qp)-0.5
       AN       =REAL(N, qp)
       ANG      =6.283185*AI/AN
       COMPUT   =AC+A*COS(ANG)+B*SIN(ANG)
       RETURN
       END

Evapo 32bits: [32.1836 33.416042 36.53447 ] in/year
Evapo 64bits: [34.042145 35.85652 38.433487] in/year
Evapo 128bits: no supported

The fact that we do not obtain the same results when using 32 and 64 bits precision seems to indicate that using the 32bits precision is not precise enough and introduce round-off errors in the calculation. In this case, using a 32bits precision seems to introduce an underestimation of evapotranspiration by 47.3 mm , 61.9 mm, and 48.1 mm for, respectively, year 1, 2, and 3. This is significant. Note that in other cases, this can cause overestimation of evapotranspiration too.

This will need to be considered very carefully. What should I do with this? Should I use 64 bits precision instead when pertinent? In my opinion, I think it would be better to use 64bits precision, even if this means that we do not obtain the same results as in the RCA test case anymore.

This topic is somewhat related to what was also discussed in PR #1

Here is the list of compiler flags for gfortran:
https://gcc.gnu.org/onlinedocs/gfortran/Fortran-Dialect-Options.html#Fortran-Dialect-Options

[jnsebgosselin]

@Rene-Lefebvre-INRS Bonjour René. Est-ce que tu pourrais regarder ce que j'ai documenté ici quand tu auras une minute stp? J'aimerais avoir ton avis sur ce sujet.

[Rene-Lefebvre-INRS]

C’est plutôt embêtant… Il semble que les tests du code aient effectivement été faits avec 32 bits (ce qui était déjà à l’époque de la double précision)… Mais cela semble insuffisant comme précision puisqu’une changement à 64 bits donne des résultats différents. Je crois aussi qu’il faudrait y aller avec 64 bits. Peux-tu compiler en double précision de façon à éviter de changer explicitement REAL pour « Double Precision » dans le code? Idéalement cela nous prendrait un test dont nous connaissons le résultat pour juger de la meilleure solution. Cela pourrait provenir d’un autre code ou d’une solution analytique…

[jnsebgosselin]

I went back to compile and test the Fortran code directly (without compiling a Python extension). I tested it when (1) compiling without any flags (32 bits precision for REAL) and when (2) compiling it with the flag -fdefault-real-8 (64 bits precision for REAL). Then, for both exe, I ran the RCA test case.

-fdefault-real-8 :

Set the default real type to an 8 byte wide type. This option also affects the kind of non-double real constants like 1.0, and does promote the default width of DOUBLE PRECISION to 16 bytes if possible, unless -fdefault-double-8 is given, too. Unlike -freal-4-real-8, it does not promote variables with explicit kind declaration.

I compared the results and they were about the same in both cases. This suggests that the default 32bits precision for REAL is enough. So I went back to investigate a little bit more to understand where the aforementioned computational errors comes from.

The problem is that when we decide to explicitly set the typespec of the variables, like Double Precision, INTENT(IN) :: AC in the COMPUT function, we need to do it over all the chain of code were the variable AC is used or else, we get this error when compiling:

Warning: Type mismatch in argument 'ac' at (1); passed REAL(4) to REAL(8)

This is what's causing the computational errors. So the solution is simple: fixing all the vars2fortran: No typespec for argument "m" warnings should solve the problem.

I think it is better to use DOUBLE PRECISION only when needed in the code because using the -freal-4-real-8 flag to make all REAL 64bits results in a significant hit on calculation times. I haven't done any benchmark to quantify this, but I can tell this is significant. So in my opinion, this is not an option.

[Rene-Lefebvre-INRS]

OK, these tests seem convincing and I agree with the use of double precision only when needed.

[jnsebgosselin]

https://stackoverflow.com/questions/40748390/gfortran-warning-change-of-value-in-conversion-from-real8-to-integer4