1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Fabian Wachsmann
5
6 */
7
8 #include <cdi.h>
9
10 #include "process_int.h"
11 #include "cdo_options.h"
12
13 void *
EstFreq(void * process)14 EstFreq(void *process)
15 {
16
17 cdo_initialize(process);
18
19 operator_check_argc(0);
20
21 const auto dataIsUnchanged = data_is_unchanged();
22
23 const auto streamID1 = cdo_open_read(0);
24
25 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
26 const auto vlistID2 = vlistDuplicate(vlistID1);
27
28 const auto taxisID1 = vlistInqTaxis(vlistID1);
29 const auto taxisID2 = taxisDuplicate(taxisID1);
30 vlistDefTaxis(vlistID2, taxisID2);
31
32 const auto streamID2 = cdo_open_write(1);
33
34 cdo_def_vlist(streamID2, vlistID2);
35
36 const auto ntsteps = vlistNtsteps(vlistID1);
37
38 Field field;
39
40 VarList varList1;
41 varListInit(varList1, vlistID1);
42
43 int tsID = 0;
44 int fyear = 0, lyear, fmonth = 0, lmonth, dummy;
45 int step_per_year = 0, currentyear, currentmon;
46
47 while (true)
48 {
49 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
50 if (nrecs == 0) break;
51
52 if (tsID == 0)
53 {
54 cdiDecodeDate(taxisInqVdate(taxisID1), &fyear, &fmonth, &dummy);
55 currentyear = fyear;
56 currentmon = fmonth;
57 }
58 else
59 cdiDecodeDate(taxisInqVdate(taxisID1), ¤tyear, ¤tmon, &dummy);
60 if (currentyear == fyear)
61 {
62 // lymonth = currentmon;
63 step_per_year++;
64 }
65 cdo_taxis_copy_timestep(taxisID2, taxisID1);
66 cdo_def_timestep(streamID2, tsID);
67
68 for (int recID = 0; recID < nrecs; recID++)
69 {
70 int varID, levelID;
71 cdo_inq_record(streamID1, &varID, &levelID);
72 cdo_def_record(streamID2, varID, levelID);
73
74 if (dataIsUnchanged)
75 {
76 cdo_copy_record(streamID2, streamID1);
77 }
78 else
79 {
80 field.init(varList1[varID]);
81 cdo_read_record(streamID1, field);
82 cdo_write_record(streamID2, field);
83 }
84 }
85
86 tsID++;
87 }
88
89 const char *frequency = "";
90
91 if (ntsteps > 2)
92 {
93 cdo_stream_inq_timestep(streamID2, ntsteps);
94 cdiDecodeDate(taxisInqVdate(taxisID2), &lyear, &lmonth, &dummy);
95 /* First, estimation by maximal number of time steps divided by covered years between last and first time step */
96 if (Options::cdoVerbose)
97 printf("Frequency is calculated by dividing the number of time steps '%d' included in the time axis by the covered years"
98 " of the time axis\ncomputed by the difference of the year of the last time stamp '%d' and the year of the first"
99 " time stamp '%d'.\n",
100 ntsteps, lyear, fyear);
101 double covered_years = lyear - fyear + 1.0;
102 double freq = ntsteps / covered_years;
103 if (DBL_IS_EQUAL(freq, 1.))
104 frequency = "yr";
105 else if (DBL_IS_EQUAL(freq, 12.))
106 frequency = "mon";
107 else if (IS_EQUAL(freq, 365.) || IS_EQUAL(freq, 365.25) || IS_EQUAL(freq, 366.))
108 frequency = "day";
109 else if (IS_EQUAL(freq, 365. * 4) || IS_EQUAL(freq, 365.25 * 4) || IS_EQUAL(freq, 366. * 4))
110 frequency = "6hr";
111 else if (IS_EQUAL(freq, 365. * 8) || IS_EQUAL(freq, 365.25 * 8) || IS_EQUAL(freq, 366. * 8))
112 frequency = "3hr";
113 else
114 {
115 int covered_months = lmonth - fmonth + 1;
116 if (Options::cdoVerbose)
117 printf("The fraction ntsteps / covered_years = '%f' is neither 1, 12, 365, 365.25, 366 nor a multiple of 365 which "
118 "would correspond to frequencies yearly, monthly, daily or subdaily respectively.\n Next try:\n\nFrequency is "
119 "calculated by dividing the number of time steps '%d' in year '%d' by the covered months in that year '%d'.\n",
120 ntsteps / covered_years, step_per_year, fyear, covered_months);
121 if (step_per_year > 366 * 8)
122 cdo_abort("Step per year '%d' in year '%d' is bigger than 366*8 which corresponds to a frequency of sub-3hourly!"
123 " This is not yet enabled.",
124 step_per_year, fyear);
125 else
126 {
127 freq = (double) step_per_year / (double) covered_months;
128 // clang-format off
129 if (freq > 31 * 8) cdo_abort("Frequency is sub-3hourly! Not yet enabled.");
130 else if (freq > 31 * 4) frequency = "3hr";
131 else if (freq > 31) frequency = "6hr";
132 else if (freq > 1) frequency = "day";
133 else frequency = "mon";
134 // clang-format on
135 }
136 }
137 }
138 else
139 cdo_abort("For %d found timesteps no frequency can be computed - at least 3 timesteps are required.", ntsteps);
140
141 if (Options::cdoVerbose) printf("Your file indicates a frequency of '%s'.\n", frequency);
142 cdiDefAttTxt(vlistID2, CDI_GLOBAL, "frequency", 3, frequency);
143
144 cdo_stream_close(streamID1);
145 cdo_stream_close(streamID2);
146
147 vlistDestroy(vlistID2);
148
149 cdo_finish();
150
151 return nullptr;
152 }
153