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       Hi      hi           Compute the humidity index
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 HI_NAME[] = "hum_index";
24 static const char HI_LONGNAME[]
25     = "Humindex describes empirically in units of temperature how the temperature and humidity influence the wellness of a human "
26       "being. HI = T + 5/9 * (A - 10) with A = e * (6.112 * 10 ** ((7.5 * T)/(237.7 + T)) * R), T  = air temperature in degree "
27       "Celsius, R = relative humidity, e = vapour pressure. Humindex is only defined for temperatures of at least 26 degree "
28       "Celsius and relative humidity of at least 40 percent.";
29 static const char HI_UNITS[] = "Celsius";
30 
31 static const int FIRST_VAR = 0;
32 
33 static double
humidityIndex(double t,double e,double r,double missval)34 humidityIndex(double t, double e, double r, double missval)
35 {
36   static const double tmin = 26.0;
37   static const double rmin = 40.0;
38 
39   if (t < tmin || r < rmin) return missval;
40 
41   return t + (5.0 / 9.0) * ((0.01 * r * e * 6.112 * std::pow(10.0, (7.5 * t) / (237.7 + t))) - 10.0);
42 }
43 
44 static void
farexpr(Field & field1,const Field & field2,const Field & field3,double (* expression)(double,double,double,double))45 farexpr(Field &field1, const Field &field2, const Field &field3, double (*expression)(double, double, double, double))
46 {
47   const auto grid1 = field1.grid;
48   const auto grid2 = field2.grid;
49   const auto grid3 = field3.grid;
50   const auto missval1 = field1.missval;
51   const auto missval2 = field2.missval;
52   const auto missval3 = field3.missval;
53 
54   const auto len = gridInqSize(grid1);
55   if (len != gridInqSize(grid2) || len != gridInqSize(grid3)) cdo_abort("Fields have different gridsize (%s)", __func__);
56 
57   if (field1.nmiss || field2.nmiss || field3.nmiss)
58     {
59       for (size_t i = 0; i < len; i++)
60         if (DBL_IS_EQUAL(field1.vec_d[i], missval1) || DBL_IS_EQUAL(field2.vec_d[i], missval2)
61             || DBL_IS_EQUAL(field3.vec_d[i], missval3))
62           field1.vec_d[i] = missval1;
63         else
64           field1.vec_d[i] = expression(field1.vec_d[i], field2.vec_d[i], field3.vec_d[i], missval1);
65     }
66   else
67     {
68       for (size_t i = 0; i < len; i++) field1.vec_d[i] = expression(field1.vec_d[i], field2.vec_d[i], field3.vec_d[i], missval1);
69     }
70 
71   field1.nmiss = field_num_miss(field1);
72 }
73 
74 void *
Hi(void * process)75 Hi(void *process)
76 {
77   cdo_initialize(process);
78 
79   cdo_operator_add("hi", 0, 0, nullptr);
80 
81   operator_check_argc(0);
82 
83   const auto streamID1 = cdo_open_read(0);
84   const auto streamID2 = cdo_open_read(1);
85   const auto streamID3 = cdo_open_read(2);
86 
87   const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
88   const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
89   const auto vlistID3 = cdo_stream_inq_vlist(streamID3);
90 
91   const auto taxisID1 = vlistInqTaxis(vlistID1);
92   // taxisID2 = vlistInqTaxis(vlistID2);
93   // taxisID3 = vlistInqTaxis(vlistID3);
94 
95   vlist_compare(vlistID1, vlistID2, CMP_DIM);
96   vlist_compare(vlistID1, vlistID3, CMP_DIM);
97 
98   const auto gridsizemax = vlistGridsizeMax(vlistID1);
99 
100   Field field1, field2, field3;
101   field1.resize(gridsizemax);
102   field2.resize(gridsizemax);
103   field3.resize(gridsizemax);
104 
105   if (Options::cdoVerbose)
106     cdo_print("Number of timesteps: file1 %d, file2 %d, file3 %d", vlistNtsteps(vlistID1), vlistNtsteps(vlistID2),
107               vlistNtsteps(vlistID3));
108 
109   const auto vlistID4 = vlistCreate();
110   const auto gridID = vlistInqVarGrid(vlistID1, FIRST_VAR);
111   const auto zaxisID = vlistInqVarZaxis(vlistID1, FIRST_VAR);
112   const auto varID4 = vlistDefVar(vlistID4, gridID, zaxisID, TIME_VARYING);
113 
114   const auto taxisID4 = cdo_taxis_create(TAXIS_RELATIVE);
115   taxisDefTunit(taxisID4, TUNIT_MINUTE);
116   taxisDefCalendar(taxisID4, CALENDAR_STANDARD);
117   taxisDefRdate(taxisID4, 19550101);
118   taxisDefRtime(taxisID4, 0);
119   vlistDefTaxis(vlistID4, taxisID4);
120 
121   cdiDefKeyString(vlistID4, varID4, CDI_KEY_NAME, HI_NAME);
122   cdiDefKeyString(vlistID4, varID4, CDI_KEY_LONGNAME, HI_LONGNAME);
123   cdiDefKeyString(vlistID4, varID4, CDI_KEY_UNITS, HI_UNITS);
124 
125   const auto streamID4 = cdo_open_write(3);
126 
127   cdo_def_vlist(streamID4, vlistID4);
128 
129   int tsID = 0;
130   while (true)
131     {
132       const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
133       if (nrecs == 0) break;
134 
135       const auto nrecs2 = cdo_stream_inq_timestep(streamID2, tsID);
136       const auto nrecs3 = cdo_stream_inq_timestep(streamID3, tsID);
137       if (nrecs2 == 0 || nrecs3 == 0) cdo_abort("Input streams have different number of timesteps!");
138 
139       cdo_taxis_copy_timestep(taxisID4, taxisID1);
140       cdo_def_timestep(streamID4, tsID);
141 
142       for (int recID = 0; recID < nrecs; recID++)
143         {
144           int varID1, levelID1;
145           cdo_inq_record(streamID1, &varID1, &levelID1);
146           cdo_read_record(streamID1, field1.vec_d.data(), &field1.nmiss);
147 
148           int varID2, levelID2;
149           cdo_inq_record(streamID2, &varID2, &levelID2);
150           cdo_read_record(streamID2, field2.vec_d.data(), &field2.nmiss);
151 
152           int varID3, levelID3;
153           cdo_inq_record(streamID3, &varID3, &levelID3);
154           cdo_read_record(streamID3, field3.vec_d.data(), &field3.nmiss);
155 
156           if (varID1 != varID2 || varID1 != varID3 || levelID1 != levelID2 || levelID1 != levelID3)
157             cdo_abort("Input streams have different structure!");
158 
159           if (varID1 != FIRST_VAR) continue;
160 
161           field1.grid = vlistInqVarGrid(vlistID1, varID1);
162           field1.missval = vlistInqVarMissval(vlistID1, varID1);
163 
164           field2.grid = vlistInqVarGrid(vlistID2, varID2);
165           field2.missval = vlistInqVarMissval(vlistID2, varID2);
166 
167           field3.grid = vlistInqVarGrid(vlistID3, varID3);
168           field3.missval = vlistInqVarMissval(vlistID3, varID3);
169 
170           farexpr(field1, field2, field3, humidityIndex);
171 
172           cdo_def_record(streamID4, varID4, levelID1);
173           cdo_write_record(streamID4, field1.vec_d.data(), field1.nmiss);
174         }
175 
176       tsID++;
177     }
178 
179   cdo_stream_close(streamID4);
180   cdo_stream_close(streamID3);
181   cdo_stream_close(streamID2);
182   cdo_stream_close(streamID1);
183 
184   cdo_finish();
185 
186   return nullptr;
187 }
188