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 #include <cdi.h>
9
10 #include "cdo_options.h"
11 #include "cdo_cdi_wrapper.h"
12 #include "cdo_zaxis.h"
13 #include "process_int.h"
14 #include "param_conversion.h"
15 #include <mpim_grid.h>
16 #include "griddes.h"
17 #include "timer.h"
18 #include "util_files.h"
19 #include "timebase.h"
20
21 static void
print_stat(const char * sinfo,MemType memtype,int datatype,int filetype,off_t nvalues,double data_size,double file_size,double tw)22 print_stat(const char *sinfo, MemType memtype, int datatype, int filetype, off_t nvalues, double data_size, double file_size,
23 double tw)
24 {
25 nvalues /= 1000000;
26 data_size /= 1024. * 1024. * 1024.;
27
28 double rout = (tw > 0) ? nvalues / tw : -1;
29 cdo_print("%s Wrote %.1f GB of %d bit floats to %s %s, %.1f MVal/s", sinfo, data_size, memtype == MemType::Float ? 32 : 64,
30 cdi_datatype_to_str(datatype), cdi_filetype_to_str(filetype), rout);
31
32 file_size /= 1024. * 1024. * 1024.;
33
34 rout = (tw > 0) ? 1024 * file_size / tw : -1;
35 cdo_print("%s Wrote %.1f GB in %.1f seconds, total %.1f MB/s", sinfo, file_size, tw, rout);
36 }
37
38 void *
CDIwrite(void * process)39 CDIwrite(void *process)
40 {
41 const auto memtype = Options::CDO_Memtype;
42 int nvars = 10, nlevs = 0, ntimesteps = 30;
43 int zaxisID;
44 int filetype = -1, datatype = -1;
45 int nruns = 1;
46 char sinfo[64] = { 0 };
47 off_t nvalues = 0;
48 double file_size = 0, data_size = 0;
49 double tw, twsum = 0;
50
51 cdo_initialize(process);
52
53 if (Options::cdoVerbose) cdo_print("parameter: <nruns, <grid, <nlevs, <ntimesteps, <nvars>>>>>");
54
55 if (cdo_operator_argc() > 5) cdo_abort("Too many arguments!");
56
57 std::string defaultgrid = "global_.2";
58 auto gridfile = defaultgrid;
59 if (cdo_operator_argc() > 0) nruns = parameter_to_int(cdo_operator_argv(0));
60 if (cdo_operator_argc() > 1) gridfile = cdo_operator_argv(1);
61 if (cdo_operator_argc() > 2) nlevs = parameter_to_int(cdo_operator_argv(2));
62 if (cdo_operator_argc() > 3) ntimesteps = parameter_to_int(cdo_operator_argv(3));
63 if (cdo_operator_argc() > 4) nvars = parameter_to_int(cdo_operator_argv(4));
64
65 nruns = std::min(std::max(nruns, 0), 9999);
66 nlevs = std::min(std::max(nlevs, 1), 255);
67 nvars = std::max(nvars, 1);
68 ntimesteps = std::max(ntimesteps, 1);
69
70 const auto gridID = cdo_define_grid(gridfile);
71 auto gridsize = gridInqSize(gridID);
72
73 if (nlevs == 1)
74 zaxisID = zaxis_from_name("surface");
75 else
76 {
77 Varray<double> levels(nlevs);
78 for (int i = 0; i < nlevs; ++i) levels[i] = 100 * i;
79 zaxisID = zaxisCreate(ZAXIS_HEIGHT, nlevs);
80 zaxisDefLevels(zaxisID, &levels[0]);
81 }
82
83 if (Options::cdoVerbose)
84 {
85 cdo_print("nruns : %d", nruns);
86 cdo_print("gridsize : %zu", gridsize);
87 cdo_print("nlevs : %d", nlevs);
88 cdo_print("ntimesteps : %d", ntimesteps);
89 cdo_print("nvars : %d", nvars);
90 }
91
92 Varray<double> array(gridsize), xvals(gridsize), yvals(gridsize);
93
94 auto gridID2 = gridID;
95 if (gridInqType(gridID) == GRID_GME) gridID2 = gridToUnstructured(gridID, 0);
96
97 if (gridInqType(gridID) != GRID_UNSTRUCTURED && gridInqType(gridID) != GRID_CURVILINEAR) gridID2 = gridToCurvilinear(gridID, 0);
98
99 gridInqXvals(gridID2, &xvals[0]);
100 gridInqYvals(gridID2, &yvals[0]);
101
102 // Convert lat/lon units if required
103 cdo_grid_to_radian(gridID2, CDI_XAXIS, gridsize, &xvals[0], "grid center lon");
104 cdo_grid_to_radian(gridID2, CDI_YAXIS, gridsize, &yvals[0], "grid center lat");
105
106 for (size_t i = 0; i < gridsize; i++) array[i] = 2 - std::cos(std::acos(std::cos(xvals[i]) * std::cos(yvals[i])) / 1.2);
107
108 Varray3D<double> vars(nvars);
109 for (int varID = 0; varID < nvars; varID++)
110 {
111 vars[varID].resize(nlevs);
112 for (int levelID = 0; levelID < nlevs; levelID++)
113 {
114 vars[varID][levelID].resize(gridsize);
115 for (size_t i = 0; i < gridsize; ++i) vars[varID][levelID][i] = varID + array[i] * (levelID + 1);
116 }
117 }
118
119 std::vector<float> farray;
120 if (memtype == MemType::Float) farray.resize(gridsize);
121
122 const auto vlistID = vlistCreate();
123
124 for (int i = 0; i < nvars; ++i)
125 {
126 const auto varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_VARYING);
127 vlistDefVarParam(vlistID, varID, cdiEncodeParam(varID + 1, 255, 255));
128 }
129
130 const auto taxisID = cdo_taxis_create(TAXIS_RELATIVE);
131 vlistDefTaxis(vlistID, taxisID);
132
133 // vlistDefNtsteps(vlistID, 1);
134
135 for (int irun = 0; irun < nruns; ++irun)
136 {
137 double tw0 = timer_val(timer_write);
138 data_size = 0;
139 nvalues = 0;
140
141 const auto streamID = cdo_open_write(0);
142 cdo_def_vlist(streamID, vlistID);
143
144 filetype = cdo_inq_filetype(streamID);
145 datatype = vlistInqVarDatatype(vlistID, 0);
146 if (datatype == CDI_UNDEFID) datatype = CDI_DATATYPE_FLT32;
147
148 const auto julday = date_to_julday(CALENDAR_PROLEPTIC, 19870101);
149
150 double t0 = timer_val(timer_write);
151
152 for (int tsID = 0; tsID < ntimesteps; tsID++)
153 {
154 const auto vdate = julday_to_date(CALENDAR_PROLEPTIC, julday + tsID);
155 const auto vtime = 0;
156 taxisDefVdate(taxisID, vdate);
157 taxisDefVtime(taxisID, vtime);
158 cdo_def_timestep(streamID, tsID);
159
160 for (int varID = 0; varID < nvars; varID++)
161 {
162 for (int levelID = 0; levelID < nlevs; levelID++)
163 {
164 nvalues += gridsize;
165 cdo_def_record(streamID, varID, levelID);
166 if (memtype == MemType::Float)
167 {
168 for (size_t i = 0; i < gridsize; ++i) farray[i] = vars[varID][levelID][i];
169 cdo_write_record_f(streamID, &farray[0], 0);
170 data_size += gridsize * 4;
171 }
172 else
173 {
174 cdo_write_record(streamID, &vars[varID][levelID][0], 0);
175 data_size += gridsize * 8;
176 }
177 }
178 }
179
180 if (Options::cdoVerbose)
181 {
182 tw = timer_val(timer_write) - t0;
183 t0 = timer_val(timer_write);
184 cdo_print("Timestep %d: %.2f seconds", tsID + 1, tw);
185 }
186 }
187
188 cdo_stream_close(streamID);
189
190 tw = timer_val(timer_write) - tw0;
191 twsum += tw;
192
193 file_size = (double) FileUtils::size(cdo_get_stream_name(0));
194
195 if (nruns > 1) snprintf(sinfo, sizeof(sinfo), "(run %d)", irun + 1);
196
197 print_stat(sinfo, memtype, datatype, filetype, nvalues, data_size, file_size, tw);
198 }
199
200 if (nruns > 1) print_stat("(mean)", memtype, datatype, filetype, nvalues, data_size, file_size, twsum / nruns);
201
202 vlistDestroy(vlistID);
203
204 cdo_finish();
205
206 return nullptr;
207 }
208