search

TEOS-10 / GSW-C

C implementation of the Thermodynamic Equation Of Seawater - 2010 (TEOS-10)
Other
18 stars 17 forks source link

Bug in gsw_geo_strf_dyn_height_pc.c #4

Closed fnencio closed 7 years ago

fnencio commented 7 years ago

I have recently tested this function using the simple examples provided at http://www.teos-10.org/pubs/gsw/html/gsw_geo_strf_dyn_height_pc.html.

I found that the returned values are not as expected.

Looking at the code, it seems that line 60 should be changed form: prv_dyn_height_deep = (i==0)? -delta_h[0]: dyn_height_deep; to prv_dyn_height_deep = (i==0)? 0: dyn_height_deep;

Currently, the first value of dyn_height_deep used in the computation is -delta_h[0]*2 (see line 61) so that the resulting geo_strf_dyn_height_pc are all offset by -delta_h[0].

PaulMBarker commented 7 years ago

I am not sure why you think it is incorrect. You want the reference level to be the deep water were there tends to be little or no motion, this way you can see the water columns motion, see Eqn. (3.32.2) of the TEOS-10 Manual, http://www.teos-10.org/pubs/TEOS-10_Manual.pdf.

fnencio commented 7 years ago

I just looked at the results I obtained using the minimal example shown at the webpage I pointed above. With the original code the values of dynamic height I obtain are all offset by -delta_h[0] compared to the ones from the webpage.

Looking at the equivalent matlab code, in gsw_geo_strf_dyn_height_pc.m the dyn_height_deep vector is defined on line 136:

dyn_height_deep = -cumsum(delta_h);

and then used in lines 142-143 to compute geo_strf_dyn_height_pc:

geo_strf_dyn_height_pc = gsw_enthalpy_SSO_0(p_mid) + dyn_height_deep + delta_h_half

So that the first element of dyn_height used to compute geo_strf_dyn_height_pc is -delta_h[0].

In the C function however, dyn_height_deep (here a scalar) is defined on lines 60-61:

prv_dyn_height_deep = (i==0)? -delta_h[0]: dyn_height_deep; dyn_height_deep = prv_dyn_height_deep - delta_h[i];

and then used on lines 66-67:

geo_strf_dyn_height_pc[i] = gsw_enthalpy_sso_0(p_mid[i]) + dyn_height_deep + delta_h_half;

So that geo_strf_dyn_height_pc[0] is computed using a value of dyn_height_deep=-delta[0]*2, and the subsequent geo_strf_dyn_height_pc[i] with values of dyn_height_deep offset by -delta[0]

hylandg commented 7 years ago

The calculation of dyn_height_deep on lines 60-61 was incorrect and has now been fixed - commit 8410c41ba4ab6ce44ef0a9944bc50ec45b119643