inconsitency/bug in verttransform_gfs.f90

[shachinger]

Dear all,

in FLEXPART 9.2 I stumbled upon what is most probably a bug, and it is still present in the FLEXPART git dev branch.

Comparing verttransform.f90 (now seems to be verttransform_ecmwf.f90) to verttransform_gfs.f90, around line 500, one notices pieces of code which probably should be identical, but they aren't:

for ECMWF it looks like:

if (vv(nx/2-1,0,iz,n).lt.0.) then

ddpol=atan(uu(nx/2-1,0,iz,n)/ &

vv(nx/2-1,0,iz,n))+xlonr

else if (vv(nx/2-1,0,iz,n).gt.0.) then

ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &

vv(nx/2-1,0,iz,n))+xlonr

(+ + xlonr)

whereas for GFS:

if(vv(nx/2-1,0,iz,n).lt.0.) then

ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr

elseif (vv(nx/2-1,0,iz,n).gt.0.) then

ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))-xlonr

(+ - xlonr)

I would thus propose changes as you can see in the attached picture (i.e. changing to "+ +" in the GFS version, and modifying the infinity case to be consistent with both approaches to the limit).

Unfortunately, I am neither an atmospheric transport expert nor do I have an idea anymore why I fixed it this way (I looked into that like several months ago, and I think I did a comparison with corresponding code pieces in FLEXTRA).

However, I think I write this here now as "as is" and some more expert people maybe check and fix this. Thanks very much!

Stephan

---

contact: ​http://shaching.userweb.mwn.de