1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Uwe Schulzweida
5
6 */
7
8 /*
9 This module contains the following operators:
10
11 Adisit adisit compute insitu from potential temperature
12 Adisit adipot compute potential from insitu temperature
13 */
14
15 #include <cdi.h>
16
17 #include "cdo_options.h"
18 #include "cdo_zaxis.h"
19 #include "process_int.h"
20 #include "cdo_vlist.h"
21 #include "param_conversion.h"
22 #include "util_string.h"
23
24 /*
25 Transformation from potential to in situ temperature according to Bryden, 1973,
26 "New polynomials for thermal expansion, adiabatic temperature gradient and potential temperature of sea water".
27 Deep Sea Research and Oceanographic Abstracts. 20, 401-408 (GILL P.602), which gives the inverse
28 transformation for an approximate value, all terms linear in t are taken after that one newton step.
29 For the check value 8.4678516 the accuracy is 0.2 mikrokelvin.
30 */
31
32 // compute insitu temperature from potential temperature
33 static inline double
adisit(const double tpot,const double sal,const double p)34 adisit(const double tpot, const double sal, const double p)
35 {
36 constexpr double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9, a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
37 a_c1 = 8.9309E-7, a_c2 = 3.1628E-8, a_c3 = 2.1987E-10, a_d = 4.1057E-9, a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
38
39 const double qc = p * (a_a1 + p * (a_c1 - a_e1 * p));
40 const double qv = p * (a_b1 - a_d * p);
41 const double dc = 1.0 + p * (-a_a2 + p * (a_c2 - a_e2 * p));
42 const double dv = a_b2 * p;
43 const double qnq = -p * (-a_a3 + p * a_c3);
44 const double qn3 = -p * a_a4;
45
46 const double tpo = tpot;
47 const double qvs = qv * (sal - 35.0) + qc;
48 const double dvs = dv * (sal - 35.0) + dc;
49 double t = (tpo + qvs) / dvs;
50 const double fne = -qvs + t * (dvs + t * (qnq + t * qn3)) - tpo;
51 const double fst = dvs + t * (2.0 * qnq + 3.0 * qn3 * t);
52 t = t - fne / fst;
53
54 return t;
55 }
56
57 // compute potential temperature from insitu temperature
58 // Ref: Gill, p. 602, Section A3.5:Potential Temperature
59 static inline double
adipot(const double t,const double s,const double p)60 adipot(const double t, const double s, const double p)
61 {
62 constexpr double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9, a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
63 a_c1 = 8.9309E-7, a_c2 = 3.1628E-8, a_c3 = 2.1987E-10, a_d = 4.1057E-9, a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
64
65 const double s_rel = s - 35.0;
66
67 const double aa = (a_a1 + t * (a_a2 - t * (a_a3 - a_a4 * t)));
68 const double bb = s_rel * (a_b1 - a_b2 * t);
69 const double cc = (a_c1 + t * (-a_c2 + a_c3 * t));
70 const double cc1 = a_d * s_rel;
71 const double dd = (-a_e1 + a_e2 * t);
72
73 const double tpot = t - p * (aa + bb + p * (cc - cc1 + p * dd));
74
75 return tpot;
76 }
77
78 static void
calc_adisit(size_t gridsize,size_t nlevel,const Varray<double> & pressure,const FieldVector & tho,const FieldVector & sao,FieldVector & tis)79 calc_adisit(size_t gridsize, size_t nlevel, const Varray<double> &pressure, const FieldVector &tho, const FieldVector &sao,
80 FieldVector &tis)
81 {
82 // pressure units: hPa
83 // tho units: Celsius
84 // sao units: psu
85
86 for (size_t levelID = 0; levelID < nlevel; ++levelID)
87 {
88 const auto &thovec = tho[levelID].vec_d;
89 const auto &saovec = sao[levelID].vec_d;
90 auto &tisvec = tis[levelID].vec_d;
91 const auto thoMissval = tho[levelID].missval;
92 const auto saoMissval = sao[levelID].missval;
93 const auto tisMissval = tis[levelID].missval;
94 for (size_t i = 0; i < gridsize; ++i)
95 {
96 const auto isMissing = (DBL_IS_EQUAL(thovec[i], thoMissval) || DBL_IS_EQUAL(saovec[i], saoMissval));
97 tisvec[i] = isMissing ? tisMissval : adisit(thovec[i], saovec[i], pressure[levelID]);
98 }
99 }
100 }
101
102 static void
calc_adipot(size_t gridsize,size_t nlevel,const Varray<double> & pressure,const FieldVector & t,const FieldVector & s,FieldVector & tpot)103 calc_adipot(size_t gridsize, size_t nlevel, const Varray<double> &pressure, const FieldVector &t, const FieldVector &s,
104 FieldVector &tpot)
105 {
106 // pressure units: hPa
107 // t units: Celsius
108 // s units: psu
109
110 for (size_t levelID = 0; levelID < nlevel; ++levelID)
111 {
112 const auto &tvec = t[levelID].vec_d;
113 const auto &svec = s[levelID].vec_d;
114 auto &tpotvec = tpot[levelID].vec_d;
115 const auto tMissval = t[levelID].missval;
116 const auto sMissval = s[levelID].missval;
117 const auto tpotMissval = tpot[levelID].missval;
118 for (size_t i = 0; i < gridsize; ++i)
119 {
120 const auto isMissing = (DBL_IS_EQUAL(tvec[i], tMissval) || DBL_IS_EQUAL(svec[i], sMissval));
121 tpotvec[i] = isMissing ? tpotMissval : adipot(tvec[i], svec[i], pressure[levelID]);
122 }
123 }
124 }
125
126 void *
Adisit(void * process)127 Adisit(void *process)
128 {
129 int thoID = -1, saoID = -1;
130 char varname[CDI_MAX_NAME], stdname[CDI_MAX_NAME];
131
132 cdo_initialize(process);
133
134 const auto ADISIT = cdo_operator_add("adisit", 0, 0, "");
135 const auto ADIPOT = cdo_operator_add("adipot", 0, 0, "");
136
137 const auto operatorID = cdo_operator_id();
138
139 const auto pin = (cdo_operator_argc() == 1) ? parameter_to_double(cdo_operator_argv(0)) : -1.0;
140
141 const auto streamID1 = cdo_open_read(0);
142 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
143
144 const auto nvars = vlistNvars(vlistID1);
145
146 for (int varID = 0; varID < nvars; varID++)
147 {
148 auto code = vlistInqVarCode(vlistID1, varID);
149 if (code <= 0)
150 {
151 vlistInqVarName(vlistID1, varID, varname);
152 cstr_to_lower_case(varname);
153
154 int length = CDI_MAX_NAME;
155 cdiInqKeyString(vlistID1, varID, CDI_KEY_STDNAME, stdname, &length);
156 cstr_to_lower_case(stdname);
157
158 if (cdo_cmpstr(varname, "s") || cdo_cmpstr(stdname, "sea_water_salinity"))
159 code = 5;
160 else if (cdo_cmpstr(varname, "t"))
161 code = 2;
162
163 const char *cname = (operatorID == ADISIT) ? "sea_water_potential_temperature" : "sea_water_temperature";
164 if (cdo_cmpstr(stdname, cname)) code = 2;
165 }
166
167 if (code == 2 || (code == 20 && operatorID == ADIPOT))
168 thoID = varID;
169 else if (code == 5)
170 saoID = varID;
171 }
172
173 if (saoID == -1) cdo_abort("Sea water salinity not found!");
174 if (thoID == -1) cdo_abort("%s temperature not found!", (operatorID == ADISIT) ? "Potential" : "Insitu");
175
176 char units[CDI_MAX_NAME] = { 0 };
177 int length = CDI_MAX_NAME;
178 cdiInqKeyString(vlistID1, thoID, CDI_KEY_UNITS, units, &length);
179 if (length == 0) strcpy(units, "Celcius");
180
181 const auto gridID = vlistGrid(vlistID1, 0);
182 const auto gridsize = vlist_check_gridsize(vlistID1);
183
184 auto zaxisID = vlistInqVarZaxis(vlistID1, saoID);
185 const auto nlevels1 = zaxisInqSize(zaxisID);
186 zaxisID = vlistInqVarZaxis(vlistID1, thoID);
187 const auto nlevels2 = zaxisInqSize(zaxisID);
188
189 if (nlevels1 != nlevels2) cdo_abort("temperature and salinity have different number of levels!");
190 const auto nlevels = nlevels1;
191
192 Varray<double> pressure(nlevels);
193 cdo_zaxis_inq_levels(zaxisID, pressure.data());
194
195 if (pin >= 0)
196 for (int i = 0; i < nlevels; ++i) pressure[i] = pin;
197 else
198 for (int i = 0; i < nlevels; ++i) pressure[i] /= 10;
199
200 if (Options::cdoVerbose)
201 {
202 cdo_print("Level Pressure");
203 for (int i = 0; i < nlevels; ++i) cdo_print("%5d %g", i + 1, pressure[i]);
204 }
205
206 FieldVector tho(nlevels), sao(nlevels), tis(nlevels);
207 for (int levelID = 0; levelID < nlevels; ++levelID)
208 {
209 tho[levelID].resize(gridsize);
210 sao[levelID].resize(gridsize);
211 tis[levelID].resize(gridsize);
212 tho[levelID].missval = vlistInqVarMissval(vlistID1, thoID);
213 sao[levelID].missval = vlistInqVarMissval(vlistID1, saoID);
214 tis[levelID].missval = tho[levelID].missval;
215 }
216
217 int datatype = CDI_DATATYPE_FLT32;
218 if (vlistInqVarDatatype(vlistID1, thoID) == CDI_DATATYPE_FLT64 && vlistInqVarDatatype(vlistID1, saoID) == CDI_DATATYPE_FLT64)
219 datatype = CDI_DATATYPE_FLT64;
220
221 const auto vlistID2 = vlistCreate();
222 vlistDefNtsteps(vlistID2, vlistNtsteps(vlistID1));
223
224 const auto tisID2 = vlistDefVar(vlistID2, gridID, zaxisID, TIME_VARYING);
225 if (operatorID == ADISIT)
226 {
227 vlistDefVarParam(vlistID2, tisID2, cdiEncodeParam(20, 255, 255));
228 cdiDefKeyString(vlistID2, tisID2, CDI_KEY_NAME, "to");
229 cdiDefKeyString(vlistID2, tisID2, CDI_KEY_LONGNAME, "Sea water temperature");
230 cdiDefKeyString(vlistID2, tisID2, CDI_KEY_STDNAME, "sea_water_temperature");
231 }
232 else
233 {
234 vlistDefVarParam(vlistID2, tisID2, cdiEncodeParam(2, 255, 255));
235 cdiDefKeyString(vlistID2, tisID2, CDI_KEY_NAME, "tho");
236 cdiDefKeyString(vlistID2, tisID2, CDI_KEY_LONGNAME, "Sea water potential temperature");
237 cdiDefKeyString(vlistID2, tisID2, CDI_KEY_STDNAME, "sea_water_potential_temperature");
238 }
239 cdiDefKeyString(vlistID2, tisID2, CDI_KEY_UNITS, units);
240 vlistDefVarMissval(vlistID2, tisID2, tis[0].missval);
241 vlistDefVarDatatype(vlistID2, tisID2, datatype);
242
243 const auto saoID2 = vlistDefVar(vlistID2, gridID, zaxisID, TIME_VARYING);
244 vlistDefVarParam(vlistID2, saoID2, cdiEncodeParam(5, 255, 255));
245 cdiDefKeyString(vlistID2, saoID2, CDI_KEY_NAME, "s");
246 cdiDefKeyString(vlistID2, saoID2, CDI_KEY_LONGNAME, "Sea water salinity");
247 cdiDefKeyString(vlistID2, saoID2, CDI_KEY_STDNAME, "sea_water_salinity");
248 cdiDefKeyString(vlistID2, saoID2, CDI_KEY_UNITS, "psu");
249 vlistDefVarMissval(vlistID2, saoID2, sao[0].missval);
250 vlistDefVarDatatype(vlistID2, saoID2, datatype);
251
252 const auto taxisID1 = vlistInqTaxis(vlistID1);
253 const auto taxisID2 = taxisDuplicate(taxisID1);
254 vlistDefTaxis(vlistID2, taxisID2);
255
256 const auto streamID2 = cdo_open_write(1);
257 cdo_def_vlist(streamID2, vlistID2);
258
259 int tsID = 0;
260 while (true)
261 {
262 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
263 if (nrecs == 0) break;
264
265 cdo_taxis_copy_timestep(taxisID2, taxisID1);
266 cdo_def_timestep(streamID2, tsID);
267
268 for (int recID = 0; recID < nrecs; ++recID)
269 {
270 int varID, levelID;
271 cdo_inq_record(streamID1, &varID, &levelID);
272 if (varID == thoID) cdo_read_record(streamID1, tho[levelID].vec_d.data(), &tho[levelID].nmiss);
273 if (varID == saoID) cdo_read_record(streamID1, sao[levelID].vec_d.data(), &sao[levelID].nmiss);
274
275 if (varID == thoID)
276 {
277 constexpr double MIN_T = -10.0;
278 constexpr double MAX_T = 40.0;
279 const auto mm = field_min_max(tho[levelID]);
280 if (mm.min < MIN_T || mm.max > MAX_T)
281 cdo_warning("Temperature in degree Celsius out of range (min=%g max=%g) [timestep:%d levelIndex:%d]!", mm.min,
282 mm.max, tsID + 1, levelID + 1);
283 }
284 }
285
286 if (operatorID == ADISIT)
287 calc_adisit(gridsize, nlevels, pressure, tho, sao, tis);
288 else
289 calc_adipot(gridsize, nlevels, pressure, tho, sao, tis);
290
291 for (int levelID = 0; levelID < nlevels; ++levelID)
292 {
293 cdo_def_record(streamID2, tisID2, levelID);
294 cdo_write_record(streamID2, tis[levelID].vec_d.data(), field_num_miss(tis[levelID]));
295 }
296
297 for (int levelID = 0; levelID < nlevels; ++levelID)
298 {
299 cdo_def_record(streamID2, saoID2, levelID);
300 cdo_write_record(streamID2, sao[levelID].vec_d.data(), field_num_miss(sao[levelID]));
301 }
302
303 tsID++;
304 }
305
306 cdo_stream_close(streamID2);
307 cdo_stream_close(streamID1);
308
309 vlistDestroy(vlistID2);
310
311 cdo_finish();
312
313 return nullptr;
314 }
315