Thank you Yvan, your information was valuable. I found the error, it ended up being something quite childish from my part. I was thinking all the time in the classic formula for 2D airfoils where instead of using the surface for adimensionalization, we use the chord. However, I forgot that despite I'm trying to solve a 2D case, I must model it as a 3D case, and therefore I do have to take into account the airfoil surface and not only its chord, even if the depth is minimal and with only one element. The depth of my mesh was not 1 but a tenth of the chord, a value chosen so to be sure the nonplanar velocity components would be zero (besides symmetry condition for the lateral walls), but that small depth was not taken into account in my coefficient formula.
As my original chord is 0.1524m, the depth becomes 0.01524m and we can find my error factor dividing by this depth: 0.01524^(-1)=65.617... yes, the scaling factor I was lacking of

. So, here below I write the corrected formula:
DO II=1,NDIM
XCOF(II) = xfor(II) / ((1.0d0/2.0d0) * RO0(1) * UREF(1)**2 * ALMAX(1)**2/10.0d0)
ENDDO
In my case I set the reference length ALMAX(IPHAS) to be my chord length, that's the reason of its inclussion in the formula.
I wasn't having that problem with the pressures for the simple reason the pressure is already divided by the surface, the whole surface and not just the chord.
Thanks to all of you and my apologies for making the thread unnecessarily long, I should have payed more attention.