Vapor Saturation (QSat) table computation Fortran vs C f32 issues

[FlorianDeconinck]

QSat, saturation of vapor, is based on the calculation of Ice and Liquid tables estimated at an incremental temperature values. The table is filled using an exact formulation for both water species described in our code in pyMoist.saturation.qsat_liquid and pyMoist.saturation.qsat_ice.

Code has been heavily debugged to be correct but numeric doesn't validate, which bubbles to qsat itself being wrong and this code is at the base of a lot of moist computation. We side-stepped the validation issue by loading the table from netcdf.

Debugging further we figured out that the issues boils down to actual C (python) vs Fortran numerical differences. Here's two calculation that don't agree with all compilation flags deactivated

C:

#include <stdio.h>

double S16 = 0.516000335e-11 * 100.0;
double S15 = 0.276961083e-8 * 100.0;
double S14 = 0.623439266e-6 * 100.0;
double S13 = 0.754129933e-4 * 100.0;
double S12 = 0.517609116e-2 * 100.0;
double S11 = 0.191372282e0 * 100.0;
double S10 = 0.298152339e1 * 100.0;
double S26 = 0.314296723e-10 * 100.0;
double S25 = 0.132243858e-7 * 100.0;
double S24 = 0.236279781e-5 * 100.0;
double S23 = 0.230325039e-3 * 100.0;
double S22 = 0.129690326e-1 * 100.0;
double S21 = 0.401390832e0 * 100.0;
double S20 = 0.535098336e1 * 100.0;

const float TI = 178.16f;
const float TT = TI - 273.16f;

int main(int argc, void *argv)
{
    float EX = TT * (TT * (TT * (TT * (TT * (TT * S16 + S15) + S14) + S13) + S12) + S11) + S10;
    printf("EX: %f\n", EX);
    return 0;
}

prints:

EX: 0.003803

F90:

program qsat
    real*8 :: S16= 0.516000335E-11*100.0
    real*8 :: S15= 0.276961083E-8 *100.0
    real*8 :: S14= 0.623439266E-6 *100.0
    real*8 :: S13= 0.754129933E-4 *100.0
    real*8 :: S12= 0.517609116E-2 *100.0
    real*8 :: S11= 0.191372282E+0 *100.0
    real*8 :: S10= 0.298152339E+1 *100.0
    real*8 :: S26= 0.314296723E-10*100.0
    real*8 :: S25= 0.132243858E-7 *100.0
    real*8 :: S24= 0.236279781E-5 *100.0
    real*8 :: S23= 0.230325039E-3 *100.0
    real*8 :: S22= 0.129690326E-1 *100.0
    real*8 :: S21= 0.401390832E+0 *100.0
    real*8 :: S20= 0.535098336E+1 *100.0
    
    real*4, parameter :: TI = 178.16
    real*4 :: TT = TI - 273.16
    real*4 :: EX = 0

    EX = TT * (TT * (TT * (TT * (TT * (TT * S16 + S15) + S14) + S13) + S12) + S11) + S10

    print *, EX

end program qsat

prints

3.77704715E-03

Both of those examples are gcc and gfortran with no options.

We need to figure out how we can make those agree, if at all.

Parent: #57

---

• Match C and Fortran
• Deploy "numerical fix" on the DSL stack, even it means a debug mode running slowly
• OR Document findings

---

PS:

The actual errors

For an expected value of

[FlorianDeconinck]

Digging in the assembly, we discovered that the value of the constants are different, e.g.

Language	S16
C	5.16000335000000020638850831546182856834903418530302587897e-10
F90	5.160003535564783305744640529155731201171875e-10

Considering the operations comprised of mult and add, this brings the difference. If we hardcode the Fortran assembly-level value into C, we validated

[FlorianDeconinck]

Same result under icx and ifx, with clang frontend

[FlorianDeconinck]

It seems Fortran is putting a 32--bit representation into a 64-bit registry on literals. E.g. doing the following on C gives the same result

const double S16 = 0.516000335e-11f * 100.f;

[FlorianDeconinck]

Best match

Lang	Ex
C	0.003777
F90	3.77704715E-03
Py	0.0037770470880218454

With code as follows

#include <stdio.h>

const double S16 = 0.516000335e-11f * 100.f;
const double S15 = 0.276961083e-8f * 100.f;
const double S14 = 0.623439266e-6f * 100.f;
const double S13 = 0.754129933e-4f * 100.f;
const double S12 = 0.517609116e-2f * 100.f;
const double S11 = 0.191372282e0f * 100.f;
const double S10 = 0.298152339e1f * 100.f;

const float TI = 178.16f;
const float TT = TI - 273.16f;

int main(int argc, char **argv)
{
    float EX = TT * (TT * (TT * (TT * (TT * (TT * S16 + S15) + S14) + S13) + S12) + S11) + S10;
    printf("EX: %f\n", EX);
    return 0;
}

program hello
    real*8 :: S16= 0.516000335E-11*100.0
    real*8 :: S15= 0.276961083E-8 *100.0
    real*8 :: S14= 0.623439266E-6 *100.0
    real*8 :: S13= 0.754129933E-4 *100.0
    real*8 :: S12= 0.517609116E-2 *100.0
    real*8 :: S11= 0.191372282E+0 *100.0
    real*8 :: S10= 0.298152339E+1 *100.0
    
    real*4, parameter :: TI = 178.16
    real*4, parameter :: TT = TI - 273.16
    real*4 :: EX = 0

    EX = TT * (TT * (TT * (TT * (TT * (TT * S16 + S15) + S14) + S13) + S12) + S11) + S10

    print *, 'EX', EX

end program hello 

import numpy as np

S16 = np.float64(np.float32(0.516000335e-11) * np.float32(100.0))
S15 = np.float64(np.float32(0.276961083e-8) * np.float32(100.0))
S14 = np.float64(np.float32(0.623439266e-6) * np.float32(100.0))
S13 = np.float64(np.float32(0.754129933e-4) * np.float32(100.0))
S12 = np.float64(np.float32(0.517609116e-2) * np.float32(100.0))
S11 = np.float64(np.float32(0.191372282e0) * np.float32(100.0))
S10 = np.float64(np.float32(0.298152339e1) * np.float32(100.0))
S26 = np.float64(np.float32(0.314296723e-10) * np.float32(100.0))
S25 = np.float64(np.float32(0.132243858e-7) * np.float32(100.0))
S24 = np.float64(np.float32(0.236279781e-5) * np.float32(100.0))
S23 = np.float64(np.float32(0.230325039e-3) * np.float32(100.0))
S22 = np.float64(np.float32(0.129690326e-1) * np.float32(100.0))
S21 = np.float64(np.float32(0.401390832e0) * np.float32(100.0))
S20 = np.float64(np.float32(0.535098336e1) * np.float32(100.0))

TI = np.float32(178.16)  # 273.16 - 95.0
TT = TI - np.float32(273.16)

EX = TT * (TT * (TT * (TT * (TT * (TT * S16 + S15) + S14) + S13) + S12) + S11) + S10

print(S16)
print(EX)

[FlorianDeconinck]

This was fixed in GEOS-ESM/GEOSgcm_GridComp#1007

The issues was a mix of 32 representation pushed into a 64-bit float and the fact that in python an i32 * f32 will result in a f64 which then percolate in the calculation and shift the errors into 64-bit land compared to the Fortran.

The actual line for this work was:

-            t = i * DELTA_T + TMINTBL
+            t = Float(i * DELTA_T) + TMINTBL