erf.f90 type mismatch with newer gfortran compilers

[donaldjmorton]

Hello, I'm pasting in an error report that I emailed to a few people on 11 February 2013. I see that the error is still present in the most recent v9.2 downloaded from flexpart.eu.

I listed priority as "major" because if you're using the newer gfortran, well, it's not going to compile. Other than that, though, it's probably not a big deal.

--- Email from 11 Feb 2013 ---

I've been installing and running the latest FLEXPART v9.02, compiled with gfortran 4.7.1 and 4.7.2.

As I mentioned to Delia and Jerome recently, the newer gfortran compilers are very strict (I think starting somewhere around 4.6 or so). I spent a lot of time last summer and autumn working with NCEP on some porting to gfortran. Stuff that would compile and run fine with Intel, Portland Group, and even older versions of gfortran had to be modified to compile and run correctly with the newest gfortran compilers. In my humble opinion, this is generally a good thing - it means that if you can get your code through the gfortran 4.7.1+ gauntlet, you might have some pretty solid code. Always better to have problems caught by the compiler than to have to spend lots of time debugging elusive runtime problems.

So, in compiling FLEXPART v9.02 tonight, I ran into some type mismatch issues in erf.f90 where double precision values were being stored in single precision variables. I know, it's all petty, but the newer gfortran wants you to be aware of these. My "quick fix" was just to declare an offending external function as type double precision and convert an argument to double precision before passing into a function.

diff erf.f90 erf.f90.ORIG
106c106
< double precision, external :: gammln
---

real, external

gammln

111c111
< gln=gammln(DBLE(a))
---

gln=gammln(a)

140c140
< double precision, external :: gammln
---

real, external

gammln

145c145
< gln=gammln(DBLE(a))
---

gln=gammln(a)

I don't feel like the changes I made adhered to your paradigm of using a "dp" KIND as defined in par_mod.f90, but I didn't want to make assumptions. I just wanted a quick fix, which I would notify you of, and then let you decide if/how you want to modify things.

[saeck]

Thanks for spotting this Don, we will add it to the new version!

[ignacio]

The fix was added to the development version. It will be available in the version 9.2.

[pesei]

The fix is not in 9.2beta.
However, for met it compiles without errors using gfortran from gcc version 4.8.2 (Debian 4.8.2-1)

Is the issue still relevant?

[cillianjoy]

Replying to pesei:

The fix is not in 9.2beta.
However, for met it compiles without errors using gfortran from gcc version 4.8.2 (Debian 4.8.2-1)

Is the issue still relevant?

I recently had the same issue. Running GNU Fortran (Ubuntu 4.8.2-19ubuntu1) 4.8.2.
I edited erf.f90 as above.

[cillianjoy]

I recently had the same issue. Running GNU Fortran (Ubuntu 4.8.2-19ubuntu1) 4.8.2.
I edited erf.f90 as above.

With FLEXPART 9.0.2.

[pesei]

Ignacio, please post your fix here, or mail it to me, I'll integrate it into v9.2!

[tjaya]

kindly add (kind=8) at two lines in erf.90, where the below string appears:

real, external

gammln

It should become

real(kind=8), external

gammln

And then compile.
Good luck.

Replying to donaldjmorton:

Hello, I'm pasting in an error report that I emailed to a few people on 11 February 2013. I see that the error is still present in the most recent v9.2 downloaded from flexpart.eu.

I listed priority as "major" because if you're using the newer gfortran, well, it's not going to compile. Other than that, though, it's probably not a big deal.

--- Email from 11 Feb 2013 ---

I've been installing and running the latest FLEXPART v9.02, compiled with gfortran 4.7.1 and 4.7.2.

As I mentioned to Delia and Jerome recently, the newer gfortran compilers are very strict (I think starting somewhere around 4.6 or so). I spent a lot of time last summer and autumn working with NCEP on some porting to gfortran. Stuff that would compile and run fine with Intel, Portland Group, and even older versions of gfortran had to be modified to compile and run correctly with the newest gfortran compilers. In my humble opinion, this is generally a good thing - it means that if you can get your code through the gfortran 4.7.1+ gauntlet, you might have some pretty solid code. Always better to have problems caught by the compiler than to have to spend lots of time debugging elusive runtime problems.

So, in compiling FLEXPART v9.02 tonight, I ran into some type mismatch issues in erf.f90 where double precision values were being stored in single precision variables. I know, it's all petty, but the newer gfortran wants you to be aware of these. My "quick fix" was just to declare an offending external function as type double precision and convert an argument to double precision before passing into a function.

diff erf.f90 erf.f90.ORIG
106c106
< double precision, external :: gammln
---

real, external

gammln

111c111
< gln=gammln(DBLE(a))
---

gln=gammln(a)

140c140
< double precision, external :: gammln
---

real, external

gammln

145c145
< gln=gammln(DBLE(a))
---

gln=gammln(a)

I don't feel like the changes I made adhered to your paradigm of using a "dp" KIND as defined in par_mod.f90, but I didn't want to make assumptions. I just wanted a quick fix, which I would notify you of, and then let you decide if/how you want to modify things.

[ahilboll]

Here's the patch (just a pretty-printed copy from the original ticket) to erf.f90 which makes the code compile on my machine (gfortran 4.9.3):

106c106
<   real, external :: gammln
---
>   double precision, external :: gammln
111c111
<   gln=gammln(a)
---
>   gln=gammln(DBLE(a))
140c140
<   real, external :: gammln
---
>   double precision, external :: gammln
145c145
<   gln=gammln(a)
---
>   gln=gammln(DBLE(a))

[pesei]

To finish this old problem, I have done the following changes in erf.f90 which implment kind=dp where needed, and also make some parts of the code a bit nicer to read. The diff is

45c45
<   real(KIND=dp) :: cof(6) = (/ &
---
>   real(kind=dp) :: cof(6) = (/ &
48,49c48,49
<   real(KIND=dp) :: stp = 2.50662827465_dp
<   real(KIND=dp) :: half = 0.5_dp, one = 1.0_dp, fpf = 5.5_dp
---
>   real(kind=dp) :: stp = 2.50662827465_dp
>   real(kind=dp) :: half = 0.5_dp, one = 1.0_dp, fpf = 5.5_dp
58a59
>   
59a61
>   
68c70
<   if(x .lt. 0. .or. a .le. 0.) then
---
>   if (x .lt. 0. .or. a .le. 0.) then
72c74
<   if(x.lt.a+1.)then
---
>   if (x .lt. a+1.) then
78a81
>   
87c90
<   if(x.lt.0..or.a.le.0.) then
---
>   if (x .lt. 0. .or. a .le. 0.) then
91c94
<   if(x.lt.a+1.)then
---
>   if (x .lt. a+1.) then
97a101
>   
101a106,107
>   use par_mod, only: dp
> 
106c112
<   real, external :: gammln
---
>   real(kind=dp), external :: gammln
111,113c117,120
<   gln=gammln(a)
<   if(x.le.0.)then
<     if(x.lt.0.) then
---
> 
>   gln=gammln(real(a,dp))
>   if (x .le. 0.) then
>     if (x .lt. 0.) then
127c134
<     if(abs(del).lt.abs(summ)*eps)go to 1
---
>     if (abs(del) .lt. abs(summ)*eps) goto 1
131c138,140
< 1   gamser=summ*exp(-x+a*log(x)-gln)
---
>   
> 1 gamser=summ*exp(-x+a*log(x)-gln)
> 
135a145,146
>   use par_mod, only: dp
> 
140c151
<   real, external :: gammln
---
>   real(kind=dp), external :: gammln
145c156
<   gln=gammln(a)
---
>   gln=gammln(real(a,kind=dp))
160c171
<     if(a1.ne.0.)then
---
>     if (a1 .ne. 0.) then
163c174
<       if(abs((g-gold)/g).lt.eps)go to 1
---
>       if (abs((g-gold)/g) .lt. eps) goto 1
169c180,182
< 1   gammcf=exp(-x+a*alog(x)-gln)*g
---
> 
> 1 gammcf=exp(-x+a*alog(x)-gln)*g
> 
193,194c206,207
<   if(x.lt.0.)then
<     erfc=1.+gammp(.5,x**2)
---
>   if (x .lt. 0.) then
>     erfc=1.+gammp(0.5,x**2)
196c209
<     erfc=gammq(.5,x**2)
---
>     erfc=gammq(0.5,x**2)
197a211
> 
211c225,226
<   if (x.lt.0.) erfcc=2.-erfcc
---
>   if (x .lt. 0.) erfcc=2.-erfcc
>   

fixed in flexpart_9.0.3

[pesei]

Ignacio, this is not in v10! Please correct.