diff --git a/gsw_check_data.h b/gsw_check_data.h
index 8a288c7..07b6b81 100644
--- a/gsw_check_data.h
+++ b/gsw_check_data.h
@@ -5,7 +5,7 @@
 **  version_date: 15th_May_2011
 **  version_number: 3.06.16
 **  mat_zip_version: 3_06_16
-**  This file created on: 2024-02-12 11:26:21.045021
+**  This file created on: 2024-09-12 12:38:26.783150
 */
 
 /*
@@ -24,6 +24,8 @@
 #define	cast_ice_n	3
 #define	cast_mpres_m	44
 #define	cast_mpres_n	3
+#define	interp_m	51
+#define	interp_n	3
 
 static UNUSED double	ct[135] = {
 27.996436412058213, 27.993857176639086, 27.944017615967979, 27.948372399994319,
@@ -296,6 +298,13 @@ static UNUSED double	p[135] = {
 8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90
 };
 
+static UNUSED double	p_i[51] = {
+0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400,
+1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200, 2300, 2400, 2500, 2600, 2700,
+2800, 2900, 3000, 3100, 3200, 3300, 3400, 3500, 3600, 3700, 3800, 3900, 4000,
+4100, 4200, 4300, 4400, 4500, 4600, 4700, 4800, 4900, 5000
+};
+
 static UNUSED double	delta_p[135] = {
 0, 10, 10, 10, 10, 10, 26, 25, 25, 25, 25, 26, 50, 51, 50, 51, 101, 101, 101,
 101, 101, 101, 101, 102, 101, 102, 101, 254, 254, 254, 255, 255, 256, 255, 256,
@@ -7933,6 +7942,200 @@ static UNUSED double	rsubrho[132] = {
 
 #define	rsubrho_ca	1.690921180852456e-08
 
+static UNUSED double	sai_sactinterp[153] = {
+34.468236430490606, 34.981684011594275, 34.845532135029757, 34.594880552563986,
+34.646910055942165, 34.658726028373025, 34.667610684451823, 34.681241443337626,
+34.693152815925814, 34.704242814066355, 34.715310018759446, 34.727191420710305,
+34.737970181514079, 34.749908756622247, 34.75865652413632, 34.767977853907084,
+34.776667058040957, 34.784334257348284, 34.79178015056219, 34.799429517990411,
+34.806552614204548, 34.812691689382405, 34.818169026804789, 34.822997593767191,
+34.827267561289105, 34.83107566345695, 34.834475961281619, 34.837525533487018,
+34.840341326091348, 34.843005908156698, 34.845452404272699, 34.847569222979239,
+34.849391445147802, 34.851249802860892, 34.853492505612103, 34.855706803475442,
+34.857247772994597, 34.858497276933569, 34.859092057229148, 34.858381605762077,
+34.856823030966915, 34.855679183654935, 34.854739459543588, 34.854308742400669,
+34.855299304684273, 34.857513577115284, 34.85900682590129, 34.859483671997552,
+34.859866852752297, 34.860716166459341, 34.862075375111715, 34.556977983037378,
+35.143552028278961, 34.743955515099316, 34.823482547808247, 34.804942319288521,
+34.760842876649718, 34.73197711758587, 34.714860568930348, 34.713553992644862,
+34.718238136890584, 34.723558346430032, 34.732932331364822, 34.741835167338756,
+34.751656516251295, 34.760987253055077, 34.770634831831927, 34.779831462444207,
+34.787896487807899, 34.795973697286236, 34.804846003583471, 34.812671598162424,
+34.817995350850836, 34.822145849128248, 34.826452359442619, 34.831415297815248,
+34.836065183487634, 34.839513990464418, 34.842264529794129, 34.845044087445629,
+34.848173560330608, 34.851227079351254, 34.853788960522365, 34.856002592079761,
+34.858008149143096, 34.859816041891456, 34.861511533268335, 34.863314423422317,
+34.865265218099864, 34.86708175309419, 34.868590932972623, 34.869910820929661,
+34.871133587707021, 34.872546396459661, 34.873351690456062, 34.872775691608311,
+34.871547239164059, 34.870721688384769, 34.870383374148055, 34.870265396867914,
+34.870338679141682, 34.870546744141642, 6.6699043409245728, 10.346695553446363,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90
+};
+
+#define	sai_sactinterp_ca	1.3001510978938313e-10
+
+static UNUSED double	cti_sactinterp[153] = {
+27.996436412058227, 25.533144569030263, 16.019861520340768, 10.341891502382861,
+8.4456966374532545, 7.2721994214211767, 6.5074004434211012, 5.8840829086292556,
+5.3600631828344829, 4.8832999981632534, 4.4343297157729946, 4.0194186111174233,
+3.6464959786663744, 3.3302175416713258, 3.0360389453491519, 2.7908185312784872,
+2.5904631148025259, 2.423166010494068, 2.2715684781542738, 2.1281996714820086,
+2.0038859800808941, 1.9057066559935123, 1.8235890129026773, 1.7482666532235847,
+1.6762447258413102, 1.6115589460702915, 1.5572356514838051, 1.5101427761616257,
+1.4678690688374059, 1.4291052039518444, 1.3940324397786752, 1.3631384058258678,
+1.3360645039300703, 1.3092326083968704, 1.2796554533485569, 1.2504454968981085,
+1.2262396689330433, 1.2054680503482365, 1.187164480536496, 1.1684829755830797,
+1.1477518433658003, 1.1286972420975527, 1.1093760918020004, 1.0943939612593574,
+1.0873483431052355, 1.0817477224915897, 1.0745478755635225, 1.0634152009122975,
+1.0515515736003886, 1.0426620383447947, 1.0351918856149576, 27.322890736362599,
+22.011254836061603, 11.560227514375722, 9.8757443413275237, 9.0022892545844027,
+8.05506023421108, 7.1412969334874559, 6.2985045169753038, 5.6689611002002005,
+5.1619651842686372, 4.6499362388415042, 4.3058993935480157, 3.9759763533444907,
+3.6420020781766036, 3.3196655649536524, 2.9994634344611657, 2.7820843247989426,
+2.6147405296234192, 2.4435788959532885, 2.2471193199195505, 2.0822572605932823,
+1.9899103600091417, 1.9259139687763396, 1.8589909978970756, 1.7804728271451022,
+1.705937917028272, 1.6475130925660368, 1.5987301233149647, 1.5517208830746407,
+1.5039888179174263, 1.4579677965960154, 1.4151089811281323, 1.3749691776165345,
+1.3359607264588433, 1.2978547591325253, 1.2603848521798859, 1.2218814621548661,
+1.1820088935554813, 1.1431502947536156, 1.1066031829996259, 1.0714954834787538,
+1.0377190422510292, 1.0030329830608418, 0.97360861725787962,
+0.9534590677768886, 0.93592175165578984, 0.92023524972450943,
+0.90565719701788883, 0.8921688460133359, 0.8795567615391916,
+0.86702210345731012, 10.502767735302228, 4.5888811896604285,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90
+};
+
+#define	cti_sactinterp_ca	6.2611249518340628e-10
+
+static UNUSED double	traceri_tracerctinterp[153] = {
+34.468236430490606, 34.981684011594275, 34.845532135029757, 34.594880552563986,
+34.646910055942165, 34.658726028373025, 34.667610684451823, 34.681241443337626,
+34.693152815925814, 34.704242814066355, 34.715310018759446, 34.727191420710305,
+34.737970181514079, 34.749908756622247, 34.75865652413632, 34.767977853907084,
+34.776667058040957, 34.784334257348284, 34.79178015056219, 34.799429517990411,
+34.806552614204548, 34.812691689382405, 34.818169026804789, 34.822997593767191,
+34.827267561289105, 34.83107566345695, 34.834475961281619, 34.837525533487018,
+34.840341326091348, 34.843005908156698, 34.845452404272699, 34.847569222979239,
+34.849391445147802, 34.851249802860892, 34.853492505612103, 34.855706803475442,
+34.857247772994597, 34.858497276933569, 34.859092057229148, 34.858381605762077,
+34.856823030966915, 34.855679183654935, 34.854739459543588, 34.854308742400669,
+34.855299304684273, 34.857513577115284, 34.85900682590129, 34.859483671997552,
+34.859866852752297, 34.860716166459341, 34.862075375111715, 34.556977983037378,
+35.143552028278961, 34.743955515099316, 34.823482547808247, 34.804942319288521,
+34.760842876649718, 34.73197711758587, 34.714860568930348, 34.713553992644862,
+34.718238136890584, 34.723558346430032, 34.732932331364822, 34.741835167338756,
+34.751656516251295, 34.760987253055077, 34.770634831831927, 34.779831462444207,
+34.787896487807899, 34.795973697286236, 34.804846003583471, 34.812671598162424,
+34.817995350850836, 34.822145849128248, 34.826452359442619, 34.831415297815248,
+34.836065183487634, 34.839513990464418, 34.842264529794129, 34.845044087445629,
+34.848173560330608, 34.851227079351254, 34.853788960522365, 34.856002592079761,
+34.858008149143096, 34.859816041891456, 34.861511533268335, 34.863314423422317,
+34.865265218099864, 34.86708175309419, 34.868590932972623, 34.869910820929661,
+34.871133587707021, 34.872546396459661, 34.873351690456062, 34.872775691608311,
+34.871547239164059, 34.870721688384769, 34.870383374148055, 34.870265396867914,
+34.870338679141682, 34.870546744141642, 6.6699043409245728, 10.346695553446363,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90
+};
+
+#define	traceri_tracerctinterp_ca	1.3001510978938313e-10
+
+static UNUSED double	cti_tracerctinterp[153] = {
+27.996436412058227, 25.533144569030263, 16.019861520340768, 10.341891502382861,
+8.4456966374532545, 7.2721994214211767, 6.5074004434211012, 5.8840829086292556,
+5.3600631828344829, 4.8832999981632534, 4.4343297157729946, 4.0194186111174233,
+3.6464959786663744, 3.3302175416713258, 3.0360389453491519, 2.7908185312784872,
+2.5904631148025259, 2.423166010494068, 2.2715684781542738, 2.1281996714820086,
+2.0038859800808941, 1.9057066559935123, 1.8235890129026773, 1.7482666532235847,
+1.6762447258413102, 1.6115589460702915, 1.5572356514838051, 1.5101427761616257,
+1.4678690688374059, 1.4291052039518444, 1.3940324397786752, 1.3631384058258678,
+1.3360645039300703, 1.3092326083968704, 1.2796554533485569, 1.2504454968981085,
+1.2262396689330433, 1.2054680503482365, 1.187164480536496, 1.1684829755830797,
+1.1477518433658003, 1.1286972420975527, 1.1093760918020004, 1.0943939612593574,
+1.0873483431052355, 1.0817477224915897, 1.0745478755635225, 1.0634152009122975,
+1.0515515736003886, 1.0426620383447947, 1.0351918856149576, 27.322890736362599,
+22.011254836061603, 11.560227514375722, 9.8757443413275237, 9.0022892545844027,
+8.05506023421108, 7.1412969334874559, 6.2985045169753038, 5.6689611002002005,
+5.1619651842686372, 4.6499362388415042, 4.3058993935480157, 3.9759763533444907,
+3.6420020781766036, 3.3196655649536524, 2.9994634344611657, 2.7820843247989426,
+2.6147405296234192, 2.4435788959532885, 2.2471193199195505, 2.0822572605932823,
+1.9899103600091417, 1.9259139687763396, 1.8589909978970756, 1.7804728271451022,
+1.705937917028272, 1.6475130925660368, 1.5987301233149647, 1.5517208830746407,
+1.5039888179174263, 1.4579677965960154, 1.4151089811281323, 1.3749691776165345,
+1.3359607264588433, 1.2978547591325253, 1.2603848521798859, 1.2218814621548661,
+1.1820088935554813, 1.1431502947536156, 1.1066031829996259, 1.0714954834787538,
+1.0377190422510292, 1.0030329830608418, 0.97360861725787962,
+0.9534590677768886, 0.93592175165578984, 0.92023524972450943,
+0.90565719701788883, 0.8921688460133359, 0.8795567615391916,
+0.86702210345731012, 10.502767735302228, 4.5888811896604285,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90, 8.9999999999999998e+90, 8.9999999999999998e+90,
+8.9999999999999998e+90
+};
+
+#define	cti_tracerctinterp_ca	6.2611249518340628e-10
+
 static UNUSED double	p_mid_tursr[132] = {
 5, 15, 25, 35, 45, 63, 88.5, 113.5, 138.5, 163.5, 189, 227, 277.5, 328, 378.5,
 454.5, 555.5, 656.5, 757.5, 858.5, 959.5, 1060.5, 1162, 1263.5, 1365, 1466.5,
diff --git a/gsw_check_data.nc b/gsw_check_data.nc
index 75b54eb..e019372 100644
Binary files a/gsw_check_data.nc and b/gsw_check_data.nc differ
diff --git a/gsw_check_functions.c b/gsw_check_functions.c
index a396a4e..b8bab4b 100644
--- a/gsw_check_functions.c
+++ b/gsw_check_functions.c
@@ -130,7 +130,8 @@ main(int argc, char **argv)
         int     count = cast_m*cast_n, i, j, k, l, n;
         double  saturation_fraction, value[cast_m*cast_n], lat[cast_m*cast_n],
                 lon[cast_m*cast_n], val1[cast_m*cast_n], val2[cast_m*cast_n], val3[cast_m*cast_n],
-                val4[cast_m*cast_n], val5[cast_m*cast_n];
+                val4[cast_m*cast_n], val5[cast_m*cast_n], val6[interp_m*interp_n],
+                val7[interp_m*interp_n];
 
         if (argc==2 && !strcmp(argv[1],"-debug"))
             debug       = 1;
@@ -583,6 +584,38 @@ main(int argc, char **argv)
         test_func(o2sol, (sa[i],ct[i],p[i],lon[i],lat[i]), value, o2sol);
         test_func(o2sol_sp_pt, (sp[i],pt[i]), value, o2sol_sp_pt);
 
+        section_title("Vertical Interpolation");
+
+        for (j = 0; j<cast_n; j++) {
+            k = j*cast_m;
+            l = j*interp_m;
+            for (n=0; n<cast_m; n++)
+                if (isnan(sa[k+n]) || fabs(sa[k+n]) >= GSW_ERROR_LIMIT)
+                    break;
+            if (gsw_sa_ct_interp(&sa[k],&ct[k],&p[k],n,
+                p_i,interp_m,&val6[l],&val7[l]) == 1)
+                printf("gsw_sa_ct_interp returned error.\n");
+        }
+        check_accuracy("gsw_sa_ct_interp",sai_sactinterp_ca,
+                "sai_sactinterp",interp_n*interp_m, val6, sai_sactinterp);
+        check_accuracy("gsw_sa_ct_interp",cti_sactinterp_ca,
+                "cti_sactinterp",interp_n*interp_m, val7, cti_sactinterp);
+
+        for (j = 0; j<cast_n; j++) {
+            k = j*cast_m;
+            l = j*interp_m;
+            for (n=0; n<cast_m; n++)
+                if (isnan(sa[k+n]) || fabs(sa[k+n]) >= GSW_ERROR_LIMIT)
+                    break;
+            if (gsw_tracer_ct_interp(&sa[k],&ct[k],&p[k],n,
+                p_i,interp_m,9.,&val6[l],&val7[l]) == 1)
+                printf("gsw_tracer_ct_interp returned error.\n");
+        }
+        check_accuracy("gsw_tracer_ct_interp",traceri_tracerctinterp_ca,
+                "traceri_tracerctinterp",interp_n*interp_m, val6, traceri_tracerctinterp);
+        check_accuracy("gsw_tracer_ct_interp",cti_tracerctinterp_ca,
+                "cti_tracerctinterp",interp_n*interp_m, val7, cti_tracerctinterp);
+
         test_infunnel();
 
         if (gsw_error_flag)
diff --git a/gsw_data_ncdump.txt b/gsw_data_ncdump.txt
index f840ab1..4b6b4b2 100644
--- a/gsw_data_ncdump.txt
+++ b/gsw_data_ncdump.txt
@@ -15,6 +15,8 @@ dimensions:
 	test_cast_midpressure_number = 3 ;
 	test_cast_midlocation_length = 45 ;
 	test_cast_midlocation_number = 2 ;
+	test_interp_length = 51 ;
+	test_interp_number = 3 ;
 variables:
 	double p_ref(nz) ;
 		p_ref:standard_name = "pressure_reference_atlas" ;
@@ -85,6 +87,7 @@ variables:
 	double SP_chck_cast(test_cast_number, test_cast_length) ;
 	double t_chck_cast(test_cast_number, test_cast_length) ;
 	double p_chck_cast(test_cast_number, test_cast_length) ;
+	double p_i(test_interp_length) ;
 	double lat_chck_cast(test_cast_number) ;
 	double long_chck_cast(test_cast_number) ;
 	double pr(value_test_cast) ;
@@ -320,6 +323,14 @@ variables:
 		Tu:computation_accuracy = 2.19097096021414e-08 ;
 	double Rsubrho(test_cast_number, test_cast_midpressure_length) ;
 		Rsubrho:computation_accuracy = 1.70965819279445e-08 ;
+	double SAi_SACTinterp(test_interp_number, test_interp_length) ;
+		SAi_SACTinterp:computation_accuracy = 1.300151097893831e-10 ;
+	double CTi_SACTinterp(test_interp_number, test_interp_length) ;
+		CTi_SACTinterp:computation_accuracy = 6.261124951834063e-10 ;
+	double traceri_tracerCTinterp(test_interp_number, test_interp_length) ;
+		traceri_tracerCTinterp:computation_accuracy = 1.300151097893831e-10 ;
+	double CTi_tracerCTinterp(test_interp_number, test_interp_length) ;
+		CTi_tracerCTinterp:computation_accuracy = 6.261124951834063e-10 ;
 	double p_mid_TuRsr(test_cast_number, test_cast_midpressure_length) ;
 		p_mid_TuRsr:computation_accuracy = 2.30002115131356e-08 ;
 	double n2(test_cast_number, test_cast_midpressure_length) ;
diff --git a/gsw_oceanographic_toolbox.c b/gsw_oceanographic_toolbox.c
index 371bce2..1bda1f3 100644
--- a/gsw_oceanographic_toolbox.c
+++ b/gsw_oceanographic_toolbox.c
@@ -8589,6 +8589,406 @@ rr68_interp_section(int sectnum, double *sa, double *ct, double *p, int mp,
 }
 /*
 !==========================================================================
+subroutine gsw_SA_CT_interp (sa,ct,p,p_i)
+!==========================================================================
+!
+! SA and CT interpolation to p_i on a cast
+!
+!  SA   =  Absolute Salinity                                       [ g/kg ]
+!  CT   =  Conservative Temperature (ITS-90)                      [ deg C ]
+!  p    =  sea pressure                                            [ dbar ]
+!           ( i.e. absolute pressure - 10.1325 dbar )
+!  p_i  =  specific query points at which the interpolated SA_i and CT_i
+!            are required                                          [ dbar ]
+
+!
+!  SA_i  =  interpolated SA values at pressures p_i                [ g/kg ]
+!  CT_i  =  interpolated CT values at pressures p_i               [ deg C ]
+!
+!--------------------------------------------------------------------------
+*/
+int
+gsw_sa_ct_interp(double *sa, double *ct, double *p, int m,
+        double *p_i, int m_i, double *sa_i, double *ct_i)
+{
+        double  factor      = 9.,
+                rec_factor  = 1./factor,
+
+                sin_kpi_on_16[] = {
+                    1.950903220161283e-1,
+                    3.826834323650898e-1,
+                    5.555702330196022e-1,
+                    7.071067811865475e-1,
+                    8.314696123025452e-1,
+                    9.238795325112867e-1,
+                    9.807852804032304e-1
+                },
+                cos_kpi_on_16[] = {
+                    9.807852804032304e-1,
+                    9.238795325112867e-1,
+                    8.314696123025452e-1,
+                    7.071067811865476e-1,
+                    5.555702330196023e-1,
+                    3.826834323650898e-1,
+                    1.950903220161283e-1
+                };
+
+        int i, j, k, prof_len,
+            not_monotonic, unique_count, new_len, p_all_len,
+            i_min_p_obs, i_obs_plus_interp_len, i_surf_and_obs_plus_interp_len,
+            i_out_len, i_2_len, i_frozen, i_shallower, i_above, i_above_i, i_below_i;
+        int *p_idx, *p_all_idx, *i_obs_plus_interp, *i_surf_and_obs_plus_interp,
+            *i_out, *i_1, *i_2, *i_3;
+        double d, ct_f, unique_p, ct_sum, sa_sum, min_p_obs, max_p_obs;
+        double *sa_obs, *ct_obs, *p_obs, *p_i_tmp,
+               *p_sort, *sa_sort, *ct_sort, *p_all, *p_all_sort,
+               *p_obs_plus_interp, *p_surf_and_obs_plus_interp,
+               *independent_variable, *independent_variable_obs_plus_interp,
+               *scaled_sa_obs, *v_tmp, *q_tmp, *v_i, *q_i,
+               *sa_i_obs_plus_interp, *ct_i_obs_plus_interp,
+               *sa_i_limiting_obs_plus_interp, *ct_i_limiting_obs_plus_interp,
+               *ctf_i_tointerp, *sa_i_tooutput, *ct_i_tooutput;
+
+        if (m < 4)
+            return 1; // There must be at least 4 bottles'
+
+        // Check if interpolating pressure is monotonic
+        for (i=0; i<m_i-1; ++i) {
+            d = p_i[i + 1] - p_i[i];
+            if (d < 0.) {
+                return 1;
+            }
+        }
+
+        for (i=0; i<m_i; i++) {
+            sa_i[i] = NAN;
+            ct_i[i] = NAN;
+        }
+
+        // Find NaNs in profile
+        prof_len = 0;
+        for (i=0; i<m; ++i) {
+            d = sa[i] + ct[i] + p[i];
+            if (!isnan(d)) {
+                ++prof_len;
+            }
+        }
+
+        if (prof_len < 2) {
+            return 1;
+        }
+
+        sa_obs = (double *) malloc(prof_len*sizeof (double));
+        ct_obs = (double *) malloc(prof_len*sizeof (double));
+        p_obs = (double *) malloc(prof_len*sizeof (double));
+
+        prof_len = 0;
+        for (i=0; i<m; ++i) {
+            d = sa[i] + ct[i] + p[i];
+            if (!isnan(d)) {
+                sa_obs[prof_len] = sa[i];
+                ct_obs[prof_len] = ct[i];
+
+                p_obs[prof_len] = 1e-3 * round(1e3 * p[i]);
+
+                ct_f = gsw_ct_freezing_poly(sa_obs[prof_len],
+                                            p_obs[prof_len],
+                                            0.);
+                if (ct_obs[prof_len] < (ct_f - 0.1)) {
+                    ct_obs[prof_len] = ct_f;
+                }
+
+                ++prof_len;
+            }
+        }
+
+        p_i_tmp = (double *) malloc(m_i*sizeof (double));
+        for (i=0; i<m_i; ++i) {
+            p_i_tmp[i] = 1e-3 * round(1e3 * p_i[i]);
+        }
+
+        // Check if profile pressure values are monotonic
+        // If they are not, then sort pressure values along with SA and CT
+        not_monotonic = 0;
+        for (i=0; i<prof_len-1; ++i) {
+            d = p_obs[i+1] - p_obs[i];
+            if (d < 0.) {
+                ++not_monotonic;
+            }
+        }
+
+        p_idx = (int *) malloc(prof_len*sizeof (int));
+        p_sort = (double *) malloc(prof_len*sizeof (double));
+        sa_sort = (double *) malloc(prof_len*sizeof (double));
+        ct_sort = (double *) malloc(prof_len*sizeof (double));
+
+        if (not_monotonic > 0) {
+            gsw_util_sort_real(p_obs, prof_len, p_idx);
+            for (i=0; i<prof_len; ++i) {
+                p_sort[i] = p_obs[p_idx[i]];
+                sa_sort[i] = sa_obs[p_idx[i]];
+                ct_sort[i] = ct_obs[p_idx[i]];
+            }
+
+            // Once sorted, only save unique pressure values.
+            // SA and CT observations with the same pressure
+            // will be averaged.
+            new_len = 0;
+            unique_count = 1;
+            unique_p = p_sort[0];
+            sa_sum = sa_sort[0];
+            ct_sum = ct_sort[0];
+            for (i=1; i<prof_len; ++i) {
+                if (unique_p == p_sort[i]) {
+                    sa_sum += sa_sort[i];
+                    ct_sum += ct_sort[i];
+                    ++unique_count;
+                } else {
+                    p_obs[new_len] = unique_p;
+                    sa_obs[new_len] = sa_sum/(double)unique_count;
+                    ct_obs[new_len] = ct_sum/(double)unique_count;
+                    ++new_len;
+                    unique_p = p_sort[i];
+                    sa_sum = sa_sort[i];
+                    ct_sum = ct_sort[i];
+                    unique_count = 1;
+                }
+            }
+            p_obs[new_len] = unique_p;
+            sa_obs[new_len] = sa_sum/(double)unique_count;
+            ct_obs[new_len] = ct_sum/(double)unique_count;
+            ++new_len;
+            prof_len = new_len;
+        }
+
+        // Combine pressure values of observed and interpolated SA and CT
+        p_all_len = prof_len + m_i;
+
+        p_all = (double *) malloc(p_all_len*sizeof (double));
+        p_all_sort = (double *) malloc(p_all_len*sizeof (double));
+        p_all_idx = (int *) malloc(p_all_len*sizeof (int));
+
+        memcpy(p_all, p_obs, prof_len*sizeof (double));
+        memcpy(&p_all[prof_len], p_i_tmp, m_i*sizeof (double));
+        gsw_util_sort_real(p_all, p_all_len, p_all_idx);
+        for (i=0; i<p_all_len; ++i) {
+            p_all_sort[i] = p_all[p_all_idx[i]];
+        }
+        new_len = 0;
+        unique_p = p_all_sort[0];
+        for (i=1; i<p_all_len; ++i) {
+            if (unique_p != p_all_sort[i]) {
+                p_all[new_len] = unique_p;
+                ++new_len;
+                unique_p = p_all_sort[i];
+            }
+        }
+        p_all[new_len] = unique_p;
+        ++new_len;
+        p_all_len = new_len;
+
+        i_min_p_obs = 0;
+        min_p_obs = p_obs[0];
+        max_p_obs = p_obs[0];
+        for (i=1; i<prof_len; ++i) {
+            if (p_obs[i] < min_p_obs) {
+                i_min_p_obs = i;
+                min_p_obs = p_obs[i];
+            }
+            max_p_obs = max(max_p_obs, p_obs[i]);
+        }
+
+        i_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
+        i_surf_and_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
+        p_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
+        p_surf_and_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
+        i_out = (int *) malloc(p_all_len*sizeof (int));
+        i_1 = (int *) malloc(p_all_len*sizeof (int));
+        i_2 = (int *) malloc(p_all_len*sizeof (int));
+        i_3 = (int *) malloc(p_all_len*sizeof (int));
+
+        i_obs_plus_interp_len = 0;
+        i_surf_and_obs_plus_interp_len = 0;
+        for (i=0; i<p_all_len; ++i) {
+            if (p_all[i] >= min_p_obs && p_all[i] <= max_p_obs) {
+                i_obs_plus_interp[i_obs_plus_interp_len] = i;
+                p_obs_plus_interp[i_obs_plus_interp_len] = p_all[i];
+                i_obs_plus_interp_len++;
+            }
+            if (p_all[i] <= max_p_obs) {
+                i_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = i;
+                p_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = p_all[i];
+                i_surf_and_obs_plus_interp_len++;
+            }
+        }
+        i_out_len = gsw_util_intersect(p_i_tmp, m_i, p_surf_and_obs_plus_interp, i_surf_and_obs_plus_interp_len, i_out, i_1);
+        i_2_len = gsw_util_intersect(p_obs, prof_len, p_obs_plus_interp, i_obs_plus_interp_len, i_2, i_3);
+
+        independent_variable = (double *) malloc(prof_len*sizeof (double));
+        independent_variable_obs_plus_interp = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+
+        for(i=0; i<prof_len; ++i) {
+            independent_variable[i] = (float)i;
+        }
+        gsw_util_pchip_interp(p_obs, independent_variable, prof_len,
+                                p_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+
+        scaled_sa_obs = (double *) malloc(prof_len*sizeof (double));
+        v_tmp = (double *) malloc(prof_len*sizeof (double));
+        q_tmp = (double *) malloc(prof_len*sizeof (double));
+        v_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        q_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        sa_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        ct_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        sa_i_limiting_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        ct_i_limiting_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+
+        for(i=0; i<prof_len; ++i) {
+            scaled_sa_obs[i] = factor * sa_obs[i];
+        }
+
+        gsw_util_pchip_interp(independent_variable, ct_obs, prof_len,
+                                independent_variable_obs_plus_interp, ct_i_obs_plus_interp, i_obs_plus_interp_len);
+        gsw_util_pchip_interp(independent_variable, scaled_sa_obs, prof_len,
+                                independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
+
+        for(i=0; i<i_obs_plus_interp_len; ++i) {
+            sa_i_obs_plus_interp[i] = rec_factor * q_i[i];
+        }
+
+        for(k=0; k<7; ++k) {
+            for(i=0; i<prof_len; ++i) {
+                v_tmp[i] = scaled_sa_obs[i] * sin_kpi_on_16[k] + ct_obs[i] * cos_kpi_on_16[k];
+                q_tmp[i] = scaled_sa_obs[i] * cos_kpi_on_16[k] - ct_obs[i] * sin_kpi_on_16[k];
+            }
+            gsw_util_pchip_interp(independent_variable, v_tmp, prof_len,
+                                    independent_variable_obs_plus_interp, v_i, i_obs_plus_interp_len);
+            gsw_util_pchip_interp(independent_variable, q_tmp, prof_len,
+                                    independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
+            for(i=0; i<i_obs_plus_interp_len; ++i) {
+                ct_i_obs_plus_interp[i] += -q_i[i] * sin_kpi_on_16[k] + v_i[i] * cos_kpi_on_16[k];
+                sa_i_obs_plus_interp[i] += rec_factor * (q_i[i] * cos_kpi_on_16[k] + v_i[i] * sin_kpi_on_16[k]);
+            }
+        }
+
+        for(i=0; i<i_obs_plus_interp_len; ++i) {
+            sa_i_obs_plus_interp[i] *= 0.125;
+            ct_i_obs_plus_interp[i] *= 0.125;
+        }
+
+        gsw_linear_interp_sa_ct(sa_obs, ct_obs, independent_variable, prof_len,
+                                independent_variable_obs_plus_interp, i_obs_plus_interp_len,
+                                sa_i_limiting_obs_plus_interp, ct_i_limiting_obs_plus_interp);
+
+        ctf_i_tointerp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+
+        while (1) {
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                d = gsw_ct_freezing_poly(sa_i_obs_plus_interp[i], p_obs_plus_interp[i], 0.);
+                if (ct_i_limiting_obs_plus_interp[i] < (d - 0.1)) {
+                    ctf_i_tointerp[i] = ct_i_limiting_obs_plus_interp[i];
+                } else {
+                    ctf_i_tointerp[i] = d;
+                }
+            }
+
+            for (i_frozen=0; i_frozen<i_obs_plus_interp_len; ++i_frozen) {
+                if (ct_i_obs_plus_interp[i_frozen] < (ctf_i_tointerp[i_frozen] - 0.1)) {
+                    break;
+                }
+            }
+
+            if (i_frozen == i_obs_plus_interp_len) {
+                break;
+            }
+
+            for (i=0; i < i_2_len; ++i) {
+                if ((i_3[i] - i_frozen) <= 0) {
+                    i_shallower = i;
+                }
+            }
+            i_above = i_2[i_shallower];
+            i_above_i = i_3[i_shallower];
+            if ((i_above + 1) > i_3[i_2_len - 1]) {
+                i_below_i = i_3[i_2_len - 1];
+            } else {
+                i_below_i = i_3[i_above + 1];
+            }
+
+            for (i=i_above_i; i<=i_below_i; ++i) {
+                sa_i_obs_plus_interp[i] = sa_i_limiting_obs_plus_interp[i];
+                ct_i_obs_plus_interp[i] = ct_i_limiting_obs_plus_interp[i];
+                ctf_i_tointerp[i] = gsw_ct_freezing_poly(sa_i_obs_plus_interp[i], p_all[i], 0.);
+            }
+        }
+
+        sa_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+        ct_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+
+        if (min_p_obs != 0.) {
+            for (i=0; i<i_surf_and_obs_plus_interp_len; ++i) {
+                sa_i_tooutput[i] = NAN;
+                ct_i_tooutput[i] = NAN;
+                if (p_i_tmp[i] < min_p_obs) {
+                    sa_i_tooutput[i] = sa_i_obs_plus_interp[i_3[i_min_p_obs]];
+                    ct_i_tooutput[i] = ct_i_obs_plus_interp[i_3[i_min_p_obs]];
+                }
+            }
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                sa_i_tooutput[i_obs_plus_interp[i]] = sa_i_obs_plus_interp[i];
+                ct_i_tooutput[i_obs_plus_interp[i]] = ct_i_obs_plus_interp[i];
+            }
+        } else {
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                sa_i_tooutput[i] = sa_i_obs_plus_interp[i];
+                ct_i_tooutput[i] = ct_i_obs_plus_interp[i];
+            }
+
+        }
+
+        for(i=0; i<i_out_len; ++i) {
+            sa_i[i_out[i]] = sa_i_tooutput[i_1[i]];
+            ct_i[i_out[i]] = ct_i_tooutput[i_1[i]];
+        }
+
+        free(independent_variable);
+        free(independent_variable_obs_plus_interp);
+        free(scaled_sa_obs);
+        free(v_tmp);
+        free(q_tmp);
+        free(v_i);
+        free(q_i);
+        free(sa_i_obs_plus_interp);
+        free(ct_i_obs_plus_interp);
+        free(sa_i_limiting_obs_plus_interp);
+        free(ct_i_limiting_obs_plus_interp);
+        free(ctf_i_tointerp);
+        free(sa_i_tooutput);
+        free(ct_i_tooutput);
+        free(sa_obs);
+        free(ct_obs);
+        free(p_obs);
+        free(p_i_tmp);
+        free(p_idx);
+        free(p_sort);
+        free(sa_sort);
+        free(ct_sort);
+        free(p_all);
+        free(p_all_sort);
+        free(p_all_idx);
+        free(i_obs_plus_interp);
+        free(i_surf_and_obs_plus_interp);
+        free(p_obs_plus_interp);
+        free(p_surf_and_obs_plus_interp);
+        free(i_out);
+        free(i_1);
+        free(i_2);
+        free(i_3);
+
+        return 0;
+}
+/*
+!==========================================================================
 elemental function gsw_sa_freezing_estimate (p, saturation_fraction, ct, t)
 !==========================================================================
 !
@@ -10918,6 +11318,348 @@ gsw_thermobaric(double sa, double ct, double p)
 }
 /*
 !==========================================================================
+subroutine gsw_tracer_CT_interp (tracer,ct,p,p_i,factor)
+!==========================================================================
+!
+!  Tracer and CT interpolation to p_i on a cast
+!
+!  tracer   =  tracer                                                 [ ? ]
+!  CT   =  Conservative Temperature (ITS-90)                      [ deg C ]
+!  p    =  sea pressure                                            [ dbar ]
+!           ( i.e. absolute pressure - 10.1325 dbar )
+!  p_i  =  specific query points at which the interpolated tracer_i
+!          and CT_i are required                                   [ dbar ]
+!  factor  =  ratio between the ranges of Conservative Temperature
+!             and tracer in the world ocean.                          [ ? ]
+
+!
+!  tracer_i  =  interpolated tracer values at pressures p_i        [ g/kg ]
+!  CT_i  =  interpolated CT values at pressures p_i               [ deg C ]
+!
+!--------------------------------------------------------------------------
+*/
+int
+gsw_tracer_ct_interp(double *tracer, double *ct, double *p, int m,
+        double *p_i, int m_i, double factor, double *tracer_i, double *ct_i)
+{
+        double  rec_factor  = 1./factor,
+
+                sin_kpi_on_16[] = {
+                    1.950903220161283e-1,
+                    3.826834323650898e-1,
+                    5.555702330196022e-1,
+                    7.071067811865475e-1,
+                    8.314696123025452e-1,
+                    9.238795325112867e-1,
+                    9.807852804032304e-1
+                },
+                cos_kpi_on_16[] = {
+                    9.807852804032304e-1,
+                    9.238795325112867e-1,
+                    8.314696123025452e-1,
+                    7.071067811865476e-1,
+                    5.555702330196023e-1,
+                    3.826834323650898e-1,
+                    1.950903220161283e-1
+                };
+
+        int i, j, k, prof_len,
+            not_monotonic, unique_count, new_len, p_all_len,
+            i_min_p_obs, i_obs_plus_interp_len, i_surf_and_obs_plus_interp_len,
+            i_out_len, i_2_len;
+        int *p_idx, *p_all_idx, *i_obs_plus_interp, *i_surf_and_obs_plus_interp,
+            *i_out, *i_1, *i_2, *i_3;
+        double d, ct_f, unique_p, ct_sum, tracer_sum, min_p_obs, max_p_obs;
+        double *tracer_obs, *ct_obs, *p_obs, *p_i_tmp,
+               *p_sort, *tracer_sort, *ct_sort, *p_all, *p_all_sort,
+               *p_obs_plus_interp, *p_surf_and_obs_plus_interp,
+               *independent_variable, *independent_variable_obs_plus_interp,
+               *scaled_tracer_obs, *v_tmp, *q_tmp, *v_i, *q_i,
+               *tracer_i_obs_plus_interp, *ct_i_obs_plus_interp,
+               *tracer_i_tooutput, *ct_i_tooutput;
+
+        if (m < 4)
+            return 1; // There must be at least 4 bottles'
+
+        // Check if interpolating pressure is monotonic
+        for (i=0; i<m_i-1; ++i) {
+            d = p_i[i + 1] - p_i[i];
+            if (d < 0.) {
+                return 1;
+            }
+        }
+
+        for (i=0; i<m_i; i++) {
+            tracer_i[i] = NAN;
+            ct_i[i] = NAN;
+        }
+
+        // Find NaNs in profile
+        prof_len = 0;
+        for (i=0; i<m; ++i) {
+            d = tracer[i] + ct[i] + p[i];
+            if (!isnan(d)) {
+                ++prof_len;
+            }
+        }
+
+        if (prof_len < 2) {
+            return 1;
+        }
+
+        tracer_obs = (double *) malloc(prof_len*sizeof (double));
+        ct_obs = (double *) malloc(prof_len*sizeof (double));
+        p_obs = (double *) malloc(prof_len*sizeof (double));
+
+        prof_len = 0;
+        for (i=0; i<m; ++i) {
+            d = tracer[i] + ct[i] + p[i];
+            if (!isnan(d)) {
+                tracer_obs[prof_len] = tracer[i];
+                ct_obs[prof_len] = ct[i];
+
+                p_obs[prof_len] = 1e-3 * round(1e3 * p[i]);
+
+                ++prof_len;
+            }
+        }
+
+        p_i_tmp = (double *) malloc(m_i*sizeof (double));
+        for (i=0; i<m_i; ++i) {
+            p_i_tmp[i] = 1e-3 * round(1e3 * p_i[i]);
+        }
+
+        // Check if profile pressure values are monotonic
+        // If they are not, then sort pressure values along with tracer and CT
+        not_monotonic = 0;
+        for (i=0; i<prof_len-1; ++i) {
+            d = p_obs[i+1] - p_obs[i];
+            if (d < 0.) {
+                ++not_monotonic;
+            }
+        }
+
+        p_idx = (int *) malloc(prof_len*sizeof (int));
+        p_sort = (double *) malloc(prof_len*sizeof (double));
+        tracer_sort = (double *) malloc(prof_len*sizeof (double));
+        ct_sort = (double *) malloc(prof_len*sizeof (double));
+
+        if (not_monotonic > 0) {
+            gsw_util_sort_real(p_obs, prof_len, p_idx);
+            for (i=0; i<prof_len; ++i) {
+                p_sort[i] = p_obs[p_idx[i]];
+                tracer_sort[i] = tracer_obs[p_idx[i]];
+                ct_sort[i] = ct_obs[p_idx[i]];
+            }
+
+            // Once sorted, only save unique pressure values.
+            // tracer and CT observations with the same pressure
+            // will be averaged.
+            new_len = 0;
+            unique_count = 1;
+            unique_p = p_sort[0];
+            tracer_sum = tracer_sort[0];
+            ct_sum = ct_sort[0];
+            for (i=1; i<prof_len; ++i) {
+                if (unique_p == p_sort[i]) {
+                    tracer_sum += tracer_sort[i];
+                    ct_sum += ct_sort[i];
+                    ++unique_count;
+                } else {
+                    p_obs[new_len] = unique_p;
+                    tracer_obs[new_len] = tracer_sum/(double)unique_count;
+                    ct_obs[new_len] = ct_sum/(double)unique_count;
+                    ++new_len;
+                    unique_p = p_sort[i];
+                    tracer_sum = tracer_sort[i];
+                    ct_sum = ct_sort[i];
+                    unique_count = 1;
+                }
+            }
+            p_obs[new_len] = unique_p;
+            tracer_obs[new_len] = tracer_sum/(double)unique_count;
+            ct_obs[new_len] = ct_sum/(double)unique_count;
+            ++new_len;
+            prof_len = new_len;
+        }
+
+        // Combine pressure values of observed and interpolated tracer and CT
+        p_all_len = prof_len + m_i;
+
+        p_all = (double *) malloc(p_all_len*sizeof (double));
+        p_all_sort = (double *) malloc(p_all_len*sizeof (double));
+        p_all_idx = (int *) malloc(p_all_len*sizeof (int));
+
+        memcpy(p_all, p_obs, prof_len*sizeof (double));
+        memcpy(&p_all[prof_len], p_i_tmp, m_i*sizeof (double));
+        gsw_util_sort_real(p_all, p_all_len, p_all_idx);
+        for (i=0; i<p_all_len; ++i) {
+            p_all_sort[i] = p_all[p_all_idx[i]];
+        }
+        new_len = 0;
+        unique_p = p_all_sort[0];
+        for (i=1; i<p_all_len; ++i) {
+            if (unique_p != p_all_sort[i]) {
+                p_all[new_len] = unique_p;
+                ++new_len;
+                unique_p = p_all_sort[i];
+            }
+        }
+        p_all[new_len] = unique_p;
+        ++new_len;
+        p_all_len = new_len;
+
+        i_min_p_obs = 0;
+        min_p_obs = p_obs[0];
+        max_p_obs = p_obs[0];
+        for (i=1; i<prof_len; ++i) {
+            if (p_obs[i] < min_p_obs) {
+                i_min_p_obs = i;
+                min_p_obs = p_obs[i];
+            }
+            max_p_obs = max(max_p_obs, p_obs[i]);
+        }
+
+        i_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
+        i_surf_and_obs_plus_interp = (int *) malloc(p_all_len*sizeof (int));
+        p_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
+        p_surf_and_obs_plus_interp = (double *) malloc(p_all_len*sizeof (double));
+        i_out = (int *) malloc(p_all_len*sizeof (int));
+        i_1 = (int *) malloc(p_all_len*sizeof (int));
+        i_2 = (int *) malloc(p_all_len*sizeof (int));
+        i_3 = (int *) malloc(p_all_len*sizeof (int));
+
+        i_obs_plus_interp_len = 0;
+        i_surf_and_obs_plus_interp_len = 0;
+        for (i=0; i<p_all_len; ++i) {
+            if (p_all[i] >= min_p_obs && p_all[i] <= max_p_obs) {
+                i_obs_plus_interp[i_obs_plus_interp_len] = i;
+                p_obs_plus_interp[i_obs_plus_interp_len] = p_all[i];
+                i_obs_plus_interp_len++;
+            }
+            if (p_all[i] <= max_p_obs) {
+                i_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = i;
+                p_surf_and_obs_plus_interp[i_surf_and_obs_plus_interp_len] = p_all[i];
+                i_surf_and_obs_plus_interp_len++;
+            }
+        }
+        i_out_len = gsw_util_intersect(p_i_tmp, m_i, p_surf_and_obs_plus_interp, i_surf_and_obs_plus_interp_len, i_out, i_1);
+        i_2_len = gsw_util_intersect(p_obs, prof_len, p_obs_plus_interp, i_obs_plus_interp_len, i_2, i_3);
+
+        independent_variable = (double *) malloc(prof_len*sizeof (double));
+        independent_variable_obs_plus_interp = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+
+        for(i=0; i<prof_len; ++i) {
+            independent_variable[i] = (float)i;
+        }
+        gsw_util_pchip_interp(p_obs, independent_variable, prof_len,
+                                p_obs_plus_interp, independent_variable_obs_plus_interp, i_obs_plus_interp_len);
+
+        scaled_tracer_obs = (double *) malloc(prof_len*sizeof (double));
+        v_tmp = (double *) malloc(prof_len*sizeof (double));
+        q_tmp = (double *) malloc(prof_len*sizeof (double));
+        v_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        q_i = (double *) malloc(i_obs_plus_interp_len*sizeof (double));
+        tracer_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+        ct_i_obs_plus_interp = (double *)malloc(i_obs_plus_interp_len*sizeof (double));
+
+        for(i=0; i<prof_len; ++i) {
+            scaled_tracer_obs[i] = factor * tracer_obs[i];
+        }
+
+        gsw_util_pchip_interp(independent_variable, ct_obs, prof_len,
+                                independent_variable_obs_plus_interp, ct_i_obs_plus_interp, i_obs_plus_interp_len);
+        gsw_util_pchip_interp(independent_variable, scaled_tracer_obs, prof_len,
+                                independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
+
+        for(i=0; i<i_obs_plus_interp_len; ++i) {
+            tracer_i_obs_plus_interp[i] = rec_factor * q_i[i];
+        }
+
+        for(k=0; k<7; ++k) {
+            for(i=0; i<prof_len; ++i) {
+                v_tmp[i] = scaled_tracer_obs[i] * sin_kpi_on_16[k] + ct_obs[i] * cos_kpi_on_16[k];
+                q_tmp[i] = scaled_tracer_obs[i] * cos_kpi_on_16[k] - ct_obs[i] * sin_kpi_on_16[k];
+            }
+            gsw_util_pchip_interp(independent_variable, v_tmp, prof_len,
+                                    independent_variable_obs_plus_interp, v_i, i_obs_plus_interp_len);
+            gsw_util_pchip_interp(independent_variable, q_tmp, prof_len,
+                                    independent_variable_obs_plus_interp, q_i, i_obs_plus_interp_len);
+            for(i=0; i<i_obs_plus_interp_len; ++i) {
+                ct_i_obs_plus_interp[i] += -q_i[i] * sin_kpi_on_16[k] + v_i[i] * cos_kpi_on_16[k];
+                tracer_i_obs_plus_interp[i] += rec_factor * (q_i[i] * cos_kpi_on_16[k] + v_i[i] * sin_kpi_on_16[k]);
+            }
+        }
+
+        for(i=0; i<i_obs_plus_interp_len; ++i) {
+            tracer_i_obs_plus_interp[i] *= 0.125;
+            ct_i_obs_plus_interp[i] *= 0.125;
+        }
+
+        tracer_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+        ct_i_tooutput = (double *)malloc(i_surf_and_obs_plus_interp_len*sizeof (double));
+
+        if (min_p_obs != 0.) {
+            for (i=0; i<i_surf_and_obs_plus_interp_len; ++i) {
+                tracer_i_tooutput[i] = NAN;
+                ct_i_tooutput[i] = NAN;
+                if (p_i_tmp[i] < min_p_obs) {
+                    tracer_i_tooutput[i] = tracer_i_obs_plus_interp[i_3[i_min_p_obs]];
+                    ct_i_tooutput[i] = ct_i_obs_plus_interp[i_3[i_min_p_obs]];
+                }
+            }
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                tracer_i_tooutput[i_obs_plus_interp[i]] = tracer_i_obs_plus_interp[i];
+                ct_i_tooutput[i_obs_plus_interp[i]] = ct_i_obs_plus_interp[i];
+            }
+        } else {
+            for (i=0; i<i_obs_plus_interp_len; ++i) {
+                tracer_i_tooutput[i] = tracer_i_obs_plus_interp[i];
+                ct_i_tooutput[i] = ct_i_obs_plus_interp[i];
+            }
+
+        }
+
+        for(i=0; i<i_out_len; ++i) {
+            tracer_i[i_out[i]] = tracer_i_tooutput[i_1[i]];
+            ct_i[i_out[i]] = ct_i_tooutput[i_1[i]];
+        }
+
+        free(independent_variable);
+        free(independent_variable_obs_plus_interp);
+        free(scaled_tracer_obs);
+        free(v_tmp);
+        free(q_tmp);
+        free(v_i);
+        free(q_i);
+        free(tracer_i_obs_plus_interp);
+        free(ct_i_obs_plus_interp);
+        free(tracer_i_tooutput);
+        free(ct_i_tooutput);
+        free(tracer_obs);
+        free(ct_obs);
+        free(p_obs);
+        free(p_i_tmp);
+        free(p_idx);
+        free(p_sort);
+        free(tracer_sort);
+        free(ct_sort);
+        free(p_all);
+        free(p_all_sort);
+        free(p_all_idx);
+        free(i_obs_plus_interp);
+        free(i_surf_and_obs_plus_interp);
+        free(p_obs_plus_interp);
+        free(p_surf_and_obs_plus_interp);
+        free(i_out);
+        free(i_1);
+        free(i_2);
+        free(i_3);
+
+        return 0;
+}
+/*
+!==========================================================================
 subroutine gsw_turner_rsubrho(sa,ct,p,nz,tu,rsubrho,p_mid)
 !==========================================================================
 
@@ -11489,6 +12231,96 @@ gsw_z_from_p(double p, double lat, double geo_strf_dyn_height,
         return (-2.0*c/(b + sqrt(b*b - 4.0*a*c)));
 }
 
+/*
+!==========================================================================
+subroutine gsw_util_intersect (x, nx, y, ny, ix, iy)
+!==========================================================================
+!
+! Find unique values common to both x and y arrays.
+!
+!  x      :  x array of values
+!  nx     :  length of x array
+!  y      :  y array of values
+!  ny     :  length of y array
+!  xi     :  indexes of x of values common to y
+!  yi     :  indexes of y of values common to x
+!
+!  ni     :  length of ix and iy
+!
+*/
+int
+gsw_util_intersect(double *x, int nx, double *y, int ny, int *ix, int *iy)
+{
+        int *sort_ix, *sort_iy, *unique_ix, *unique_iy;
+        int i, uix, uiy, nxs, nys, itx, ity, ni;
+
+        double xx, yy;
+
+        if (nx <= 0 || ny <= 0) {
+            return 0;
+        }
+
+        sort_ix = (int *) malloc(nx * sizeof(int));
+        sort_iy = (int *) malloc(ny * sizeof(int));
+        unique_ix = (int *) malloc(nx * sizeof(int));
+        unique_iy = (int *) malloc(ny * sizeof(int));
+
+        gsw_util_sort_real(x, nx, sort_ix);
+        gsw_util_sort_real(y, ny, sort_iy);
+
+        nxs = 0;
+        uix = sort_ix[0];
+        for (i=1; i<nx; ++i) {
+            if (x[uix] != x[sort_ix[i]]) {
+                unique_ix[nxs] = uix;
+                ++nxs;
+                uix = sort_ix[i];
+            } else if (sort_ix[i] < uix) {
+                uix = sort_ix[i];
+            }
+        }
+        unique_ix[nxs] = uix;
+        ++nxs;
+
+        nys = 0;
+        uiy = sort_iy[0];
+        for (i=1; i<ny; ++i) {
+            if (y[uiy] != y[sort_iy[i]]) {
+                unique_iy[nys] = uiy;
+                ++nys;
+                uiy = sort_iy[i];
+            } else if (sort_iy[i] < uiy) {
+                uiy = sort_iy[i];
+            }
+        }
+        unique_iy[nys] = uiy;
+        ++nys;
+
+        itx = ity = ni = 0;
+        while(itx < nxs && ity < nys) {
+            xx = x[unique_ix[itx]];
+            yy = y[unique_iy[ity]];
+            if (xx < yy) {
+                ++itx;
+            } else if (xx > yy) {
+                ++ity;
+            } else {
+                ix[ni] = unique_ix[itx];
+                iy[ni] = unique_iy[ity];
+                ++itx;
+                ++ity;
+                ++ni;
+            }
+        }
+
+        free(sort_ix);
+        free(sort_iy);
+        free(unique_ix);
+        free(unique_iy);
+
+        return ni;
+}
+
 /*
 **  The End
 **!==========================================================================
diff --git a/gswteos-10.h b/gswteos-10.h
index 0d43caa..e273cad 100644
--- a/gswteos-10.h
+++ b/gswteos-10.h
@@ -228,6 +228,8 @@ DECLSPEC extern double gsw_rho_t_exact(double sa, double t, double p);
 DECLSPEC extern void   gsw_rr68_interp_sa_ct(double *sa, double *ct, double *p, int mp,
                 double *p_i, int mp_i, double *sa_i, double *ct_i);
 DECLSPEC extern double gsw_saar(double p, double lon, double lat);
+DECLSPEC extern int    gsw_sa_ct_interp(double *sa, double *ct, double *p,
+                int m, double *p_i, int m_i, double *sa_i, double *ct_i);
 DECLSPEC extern double gsw_sa_freezing_estimate(double p, double saturation_fraction,
                 double *ct, double *t);
 DECLSPEC extern double gsw_sa_freezing_from_ct(double ct, double p,
@@ -297,11 +299,14 @@ DECLSPEC extern double gsw_t_freezing_poly(double sa, double p,
 DECLSPEC extern double gsw_t_from_ct(double sa, double ct, double p);
 DECLSPEC extern double gsw_t_from_pt0_ice(double pt0_ice, double p);
 DECLSPEC extern double gsw_thermobaric(double sa, double ct, double p);
+DECLSPEC extern int    gsw_tracer_ct_interp(double *tracer, double *ct, double *p,
+                int m, double *p_i, int m_i, double factor, double *sa_i, double *ct_i);
 DECLSPEC extern void   gsw_turner_rsubrho(double *sa, double *ct, double *p, int nz,
                 double *tu, double *rsubrho, double *p_mid);
 DECLSPEC extern int    gsw_util_indx(double *x, int n, double z);
 DECLSPEC extern double *gsw_util_interp1q_int(int nx, double *x, int *iy, int nxi,
                 double *x_i, double *y_i);
+DECLSPEC extern int    gsw_util_intersect(double *x, int nx, double *y, int ny, int *ix, int *iy);
 DECLSPEC extern double *gsw_util_linear_interp(int nx, double *x, int ny, double *y,
                 int nxi, double *x_i, double *y_i);
 DECLSPEC extern void   gsw_util_sort_real(double *rarray, int nx, int *iarray);
diff --git a/make_data_from_mat.py b/make_data_from_mat.py
index dede788..6dc8b7a 100644
--- a/make_data_from_mat.py
+++ b/make_data_from_mat.py
@@ -97,6 +97,7 @@ def get_arraydims():
     ["sp", "SP_chck_cast"],
     ["t", "t_chck_cast"],
     ["p", "p_chck_cast"],
+    ["p_i", "p_i"],
     ["delta_p", "delta_p_chck_cast"],
     ["p_shallow", "p_chck_cast_shallow"],
     ["p_deep", "p_chck_cast_deep"],
@@ -271,6 +272,10 @@ def get_arraydims():
     ['p_mid_n2', ""],
     ['Tu', ""],
     ['Rsubrho', ""],
+    ['SAi_SACTinterp', ""],
+    ['CTi_SACTinterp', ""],
+    ['traceri_tracerCTinterp', ""],
+    ['CTi_tracerCTinterp', ""],
     ['p_mid_TuRsr', ""],
     ['IPVfN2', ""],
     ['p_mid_IPVfN2', ""],
@@ -380,6 +385,7 @@ def get_arraydims():
 work_dims["cast_m"], work_dims["cast_n"] = cv['p_chck_cast'].shape
 work_dims["cast_ice_m"], work_dims["cast_ice_n"] = cv["p_Arctic"].shape
 work_dims["cast_mpres_m"], work_dims["cast_mpres_n"] = cv["p_mid_n2"].shape
+work_dims["interp_m"], work_dims["interp_n"] = cv['SAi_SACTinterp'].shape
 
 for key, val in work_dims.items():
     out.write("#define\t%s\t%d\n" % (key, val))
@@ -410,6 +416,8 @@ def get_arraydims():
     test_cast_midpressure_number = 3,
     test_cast_midlocation_length = 45,
     test_cast_midlocation_number = 2,
+    test_interp_length = 51,
+    test_interp_number = 3
 )
 arraydims = get_arraydims()