1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Copyright (C) 2006 Brockmann Consult
5
6 Author: Ralf Quast
7
8 */
9
10 /*
11 This module contains the following operators:
12
13 Wct wct Compute the windchill temperature (degree C)
14 */
15
16 #include <cdi.h>
17
18 #include "cdo_options.h"
19 #include "process_int.h"
20 #include "cdo_vlist.h"
21 #include "cdo_cdi_wrapper.h"
22
23 static const char WCT_NAME[] = "wind_chill_temperature";
24 static const char WCT_LONGNAME[] = "Windchill temperature describes the fact that low temperatures are felt "
25 "to be even lower in case of wind. It is based on the rate of heat loss "
26 "from exposed skin caused by wind and cold. It is calculated according "
27 "to the empirical formula: 33 + (T - 33) * (0.478 + 0.237 * ( "
28 "SQRT(ff*3.6) - 0.0124 * ff * 3.6)) with T = air temperature in "
29 "degree Celsius, ff = 10 m wind speed in m/s. Windchill temperature is "
30 "only defined for temperatures at or below 33 degree Celsius and wind "
31 "speeds above 1.39 m/s. It is mainly used for freezing temperatures.";
32 static const char WCT_UNITS[] = "Celsius";
33
34 static const int FIRST_VAR = 0;
35
36 static double
windchillTemperature(double t,double ff,double missval)37 windchillTemperature(double t, double ff, double missval)
38 {
39 constexpr double tmax = 33.0;
40 constexpr double vmin = 1.39; /* minimum wind speed (m/s) */
41
42 return ff < vmin || t > tmax ? missval : tmax + (t - tmax) * (0.478 + 0.237 * (std::sqrt(ff * 3.6) - 0.0124 * ff * 3.6));
43 }
44
45 static void
farexpr(Field & field1,Field & field2,double (* expression)(double,double,double))46 farexpr(Field &field1, Field &field2, double (*expression)(double, double, double))
47 {
48 const auto missval1 = field1.missval;
49 const auto missval2 = field2.missval;
50
51 const auto len = field1.size;
52 if (len != field2.size) cdo_abort("Fields have different size (%s)", __func__);
53
54 if (field1.nmiss || field2.nmiss)
55 {
56 for (size_t i = 0; i < len; i++)
57 if (DBL_IS_EQUAL(field1.vec_d[i], missval1) || DBL_IS_EQUAL(field2.vec_d[i], missval2))
58 field1.vec_d[i] = missval1;
59 else
60 field1.vec_d[i] = expression(field1.vec_d[i], field2.vec_d[i], missval1);
61 }
62 else
63 {
64 for (size_t i = 0; i < len; i++) field1.vec_d[i] = expression(field1.vec_d[i], field2.vec_d[i], missval1);
65 }
66
67 field1.nmiss = field_num_miss(field1);
68 }
69
70 void *
Wct(void * process)71 Wct(void *process)
72 {
73 cdo_initialize(process);
74 cdo_operator_add("wct", 0, 0, nullptr);
75
76 operator_check_argc(0);
77
78 const auto streamID1 = cdo_open_read(0);
79 const auto streamID2 = cdo_open_read(1);
80
81 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
82 const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
83
84 const auto taxisID1 = vlistInqTaxis(vlistID1);
85
86 vlist_compare(vlistID1, vlistID2, CMP_DIM);
87
88 const auto gridsizemax = vlistGridsizeMax(vlistID1);
89
90 Field field1, field2;
91 field1.resize(gridsizemax);
92 field2.resize(gridsizemax);
93
94 if (Options::cdoVerbose) cdo_print("Number of timesteps: file1 %d, file2 %d", vlistNtsteps(vlistID1), vlistNtsteps(vlistID2));
95
96 const auto vlistID3 = vlistCreate();
97 const auto gridID = vlistInqVarGrid(vlistID1, FIRST_VAR);
98 const auto zaxisID = vlistInqVarZaxis(vlistID1, FIRST_VAR);
99 const auto varID3 = vlistDefVar(vlistID3, gridID, zaxisID, TIME_VARYING);
100
101 const auto taxisID3 = cdo_taxis_create(TAXIS_RELATIVE);
102 taxisDefTunit(taxisID3, TUNIT_MINUTE);
103 taxisDefCalendar(taxisID3, CALENDAR_STANDARD);
104 taxisDefRdate(taxisID3, 19550101);
105 taxisDefRtime(taxisID3, 0);
106 vlistDefTaxis(vlistID3, taxisID3);
107
108 cdiDefKeyString(vlistID3, varID3, CDI_KEY_NAME, WCT_NAME);
109 cdiDefKeyString(vlistID3, varID3, CDI_KEY_LONGNAME, WCT_LONGNAME);
110 cdiDefKeyString(vlistID3, varID3, CDI_KEY_UNITS, WCT_UNITS);
111
112 const auto streamID3 = cdo_open_write(2);
113 cdo_def_vlist(streamID3, vlistID3);
114
115 int tsID = 0;
116 while (true)
117 {
118 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
119 if (nrecs == 0) break;
120
121 const auto nrecs2 = cdo_stream_inq_timestep(streamID2, tsID);
122 if (nrecs2 == 0) cdo_abort("Input streams have different number of timesteps!");
123
124 cdo_taxis_copy_timestep(taxisID3, taxisID1);
125 cdo_def_timestep(streamID3, tsID);
126
127 for (int recID = 0; recID < nrecs; recID++)
128 {
129 int varID1, levelID1;
130 cdo_inq_record(streamID1, &varID1, &levelID1);
131 cdo_read_record(streamID1, field1.vec_d.data(), &field1.nmiss);
132
133 int varID2, levelID2;
134 cdo_inq_record(streamID2, &varID2, &levelID2);
135 cdo_read_record(streamID2, field2.vec_d.data(), &field2.nmiss);
136
137 if (varID1 != varID2 || levelID1 != levelID2) cdo_abort("Input streams have different structure!");
138
139 if (varID1 != FIRST_VAR) continue;
140
141 field1.missval = vlistInqVarMissval(vlistID1, varID1);
142 field2.missval = vlistInqVarMissval(vlistID2, varID2);
143
144 farexpr(field1, field2, windchillTemperature);
145
146 cdo_def_record(streamID3, varID3, levelID1);
147 cdo_write_record(streamID3, field1.vec_d.data(), field1.nmiss);
148 }
149
150 tsID++;
151 }
152
153 cdo_stream_close(streamID3);
154 cdo_stream_close(streamID2);
155 cdo_stream_close(streamID1);
156
157 cdo_finish();
158
159 return nullptr;
160 }
161