1 /*
2   This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3 
4   Author: Cedrick Ansorge
5 
6 */
7 
8 /*
9   This module contains the following operators:
10 
11   Eofcoeff             eofcoeff             process eof coefficients
12 */
13 
14 #include <cdi.h>
15 
16 #include "cdo_options.h"
17 #include "process_int.h"
18 #include "cdo_vlist.h"
19 #include "util_files.h"
20 
21 // NO MISSING VALUE SUPPORT ADDED SO FAR
22 
23 void *
Eofcoeff(void * process)24 Eofcoeff(void *process)
25 {
26   char eof_name[16], oname[1024];
27   double missval1 = -999, missval2;
28   int varID, levelID;
29   int nrecs;
30   size_t nmiss;
31 
32   cdo_initialize(process);
33 
34   if (process_self().m_ID != 0) cdo_abort("This operator can't be combined with other operators!");
35 
36   cdo_operator_add("eofcoeff", 0, 0, nullptr);
37 
38   const auto streamID1 = cdo_open_read(0);
39   const auto streamID2 = cdo_open_read(1);
40 
41   const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
42   const auto vlistID2 = cdo_stream_inq_vlist(streamID2);
43   const auto vlistID3 = vlistDuplicate(vlistID2);
44 
45   // taxisID1 = vlistInqTaxis(vlistID1);
46   const auto taxisID2 = vlistInqTaxis(vlistID2);
47   const auto taxisID3 = taxisDuplicate(taxisID2);
48 
49   const auto gridID1 = vlistInqVarGrid(vlistID1, 0);
50   const auto gridID2 = vlistInqVarGrid(vlistID2, 0);
51 
52   const auto gridsize = vlistGridsizeMax(vlistID1);
53   if (gridsize != vlistGridsizeMax(vlistID2))
54     cdo_abort("Gridsize of input files does not match! %zu and %zu", gridsize, vlistGridsizeMax(vlistID2));
55 
56   if (vlistNgrids(vlistID2) > 1 || vlistNgrids(vlistID1) > 1) cdo_abort("Too many different grids in input!");
57 
58   const auto nvars = vlistNvars(vlistID1) == vlistNvars(vlistID2) ? vlistNvars(vlistID1) : -1;
59   const auto nlevs = zaxisInqSize(vlistInqVarZaxis(vlistID1, 0));
60 
61   if (gridID1 != gridID2) cdo_compare_grids(gridID1, gridID2);
62 
63   strcpy(oname, cdo_get_obase().c_str());
64   int nchars = strlen(oname);
65 
66   auto refname = cdo_get_stream_name(0);
67   char filesuffix[32] = { 0 };
68   FileUtils::gen_suffix(filesuffix, sizeof(filesuffix), cdo_inq_filetype(streamID1), vlistID1, refname);
69 
70   FieldVector3D eof(nvars);
71   for (varID = 0; varID < nvars; varID++) eof[varID].resize(nlevs);
72 
73   int eofID = 0;
74   while (1)
75     {
76       nrecs = cdo_stream_inq_timestep(streamID1, eofID);
77       if (nrecs == 0) break;
78 
79       for (int recID = 0; recID < nrecs; recID++)
80         {
81           cdo_inq_record(streamID1, &varID, &levelID);
82           missval1 = vlistInqVarMissval(vlistID1, varID);
83           if (eofID == 0)
84             eof[varID][levelID].resize(1);
85           else
86             eof[varID][levelID].resize(eofID + 1);
87           eof[varID][levelID][eofID].grid = gridID1;
88           eof[varID][levelID][eofID].missval = missval1;
89           eof[varID][levelID][eofID].resize(gridsize);
90           varray_fill(gridsize, eof[varID][levelID][eofID].vec_d, missval1);
91 
92           if (varID >= nvars) cdo_abort("Internal error - varID >= nvars");
93           if (levelID >= nlevs) cdo_abort("Internal error - levelID >= nlevs");
94 
95           cdo_read_record(streamID1, eof[varID][levelID][eofID].vec_d.data(), &nmiss);
96           eof[varID][levelID][eofID].nmiss = nmiss;
97         }
98       eofID++;
99     }
100 
101   const int neof = eofID;
102 
103   if (Options::cdoVerbose) cdo_print("%s contains %i eof's", cdo_get_stream_name(0), neof);
104   // Create 1x1 Grid for output
105   const int gridID3 = gridCreate(GRID_LONLAT, 1);
106   gridDefXsize(gridID3, 1);
107   gridDefYsize(gridID3, 1);
108   double xvals = 0, yvals = 0;
109   gridDefXvals(gridID3, &xvals);
110   gridDefYvals(gridID3, &yvals);
111 
112   // Create var-list and time-axis for output
113 
114   const auto ngrids = vlistNgrids(vlistID3);
115   for (int i = 0; i < ngrids; i++) vlistChangeGridIndex(vlistID3, i, gridID3);
116 
117   vlistDefTaxis(vlistID3, taxisID3);
118   for (varID = 0; varID < nvars; varID++) vlistDefVarTimetype(vlistID3, varID, TIME_VARYING);
119 
120   // open streams for eofcoeff output
121   std::vector<CdoStreamID> streamIDs(neof);
122   for (eofID = 0; eofID < neof; eofID++)
123     {
124       oname[nchars] = '\0';
125 
126       sprintf(eof_name, "%5.5i", eofID);
127       strcat(oname, eof_name);
128       if (filesuffix[0]) strcat(oname, filesuffix);
129 
130       streamIDs[eofID] = cdo_open_write(oname);
131 
132       if (Options::cdoVerbose) cdo_print("opened %s ('w')  as stream%i for %i. eof", oname, streamIDs[eofID]->get_id(), eofID + 1);
133 
134       cdo_def_vlist(streamIDs[eofID], vlistID3);
135     }
136 
137   // ALLOCATE temporary fields for data read and write
138   Field in, out;
139   in.resize(gridsize);
140   in.grid = gridID1;
141   out.missval = missval1;
142   out.nmiss = 0;
143   out.resize(1);
144 
145   int tsID = 0;
146   while (1)
147     {
148       nrecs = cdo_stream_inq_timestep(streamID2, tsID);
149       if (nrecs == 0) break;
150 
151       cdo_taxis_copy_timestep(taxisID3, taxisID2);
152       /*for ( eofID=0; eofID<neof; eofID++)
153         {
154           fprintf(stderr, "defining ts %i\n", tsID);
155           cdo_def_timestep(streamIDs[eofID],tsID);
156         }
157       */
158       for (int recID = 0; recID < nrecs; recID++)
159         {
160           cdo_inq_record(streamID2, &varID, &levelID);
161           missval2 = vlistInqVarMissval(vlistID2, varID);
162           cdo_read_record(streamID2, in.vec_d.data(), &in.nmiss);
163 
164           for (eofID = 0; eofID < neof; eofID++)
165             {
166               if (recID == 0) cdo_def_timestep(streamIDs[eofID], tsID);
167               // if ( recID == 0 ) fprintf(stderr, "ts%i rec%i eof%i\n", tsID, recID, eofID);
168               out.vec_d[0] = 0;
169               out.grid = gridID3;
170               out.missval = missval2;
171               for (size_t i = 0; i < gridsize; i++)
172                 {
173                   if (!DBL_IS_EQUAL(in.vec_d[i], missval2) && !DBL_IS_EQUAL(eof[varID][levelID][eofID].vec_d[i], missval1))
174                     {
175                       // out.vec_d[0] += w[i]*in.vec_d[i]*eof[varID][levelID][eofID].vec_d[i];
176                       out.vec_d[0] += in.vec_d[i] * eof[varID][levelID][eofID].vec_d[i];
177                     }
178                 }
179               if (!DBL_IS_EQUAL(out.vec_d[0], 0.))
180                 nmiss = 0;
181               else
182                 {
183                   nmiss = 1;
184                   out.vec_d[0] = missval2;
185                 }
186 
187               cdo_def_record(streamIDs[eofID], varID, levelID);
188               // fprintf(stderr, "%d %d %d %d %d %g\n", streamIDs[eofID],tsID, recID, varID, levelID,*out.vec_d.data());
189               cdo_write_record(streamIDs[eofID], out.vec_d.data(), nmiss);
190             }
191           if (varID >= nvars) cdo_abort("Internal error - varID >= nvars");
192           if (levelID >= nlevs) cdo_abort("Internal error - levelID >= nlevs");
193         }
194 
195       tsID++;
196     }
197 
198   for (eofID = 0; eofID < neof; eofID++) cdo_stream_close(streamIDs[eofID]);
199 
200   cdo_stream_close(streamID2);
201   cdo_stream_close(streamID1);
202 
203   cdo_finish();
204 
205   return nullptr;
206 }
207