Integrated profile gives wrong values
The pressure at layer centers should not directly used. Instead, the differential pressure dp for each layer should be used in the equation below.
vcd[:, :, level] = (10.0*Ra*primary[:, :, level]*dp)/(g*dryair)
Requires to calculate the a and b coefficients at the interfaces to obtain dp.