Eofcoeff.cc 7.29 KB
Newer Older
Cedrick Ansorge's avatar
Cedrick Ansorge committed
1
2
3
/*
 This file is part of CDO. CDO is a collection of Operators to
 manipulate and analyse Climate model Data.
4

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
 Copyright (C) 2003-2018 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
Cedrick Ansorge's avatar
Cedrick Ansorge committed
6
 See COPYING file for copying and redistribution conditions.
7

Cedrick Ansorge's avatar
Cedrick Ansorge committed
8
9
10
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation; version 2 of the License.
11

Cedrick Ansorge's avatar
Cedrick Ansorge committed
12
13
14
15
16
17
18
19
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.
 */

/*
 This module contains the following operators:
20

Cedrick Ansorge's avatar
Cedrick Ansorge committed
21
22
23
 Eofcoeff             eofcoeff             process eof coefficients
*/

Ralf Mueller's avatar
Ralf Mueller committed
24
#include <cdi.h>
Oliver Heidmann's avatar
Oliver Heidmann committed
25

Cedrick Ansorge's avatar
Cedrick Ansorge committed
26
#include "cdo_int.h"
27
#include "pstream_int.h"
Cedrick Ansorge's avatar
Cedrick Ansorge committed
28
29
30

// NO MISSING VALUE SUPPORT ADDED SO FAR

31
32
void *
Eofcoeff(void *process)
Cedrick Ansorge's avatar
Cedrick Ansorge committed
33
{
34
  char eof_name[16], oname[1024], filesuffix[32];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
  double missval1 = -999, missval2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
37
  Field in;
  Field out;
38
  int i, varID, levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
  int nrecs;
40
41
  size_t nmiss;

42
  cdoInitialize(process);
43

Uwe Schulzweida's avatar
Uwe Schulzweida committed
44
  if (processSelf().m_ID != 0) cdoAbort("This operator can't be combined with other operators!");
45
46

  cdoOperatorAdd("eofcoeff", 0, 0, NULL);
47

48
49
  int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
  int streamID2 = cdoStreamOpenRead(cdoStreamName(1));
50

51
52
  int vlistID1 = cdoStreamInqVlist(streamID1);
  int vlistID2 = cdoStreamInqVlist(streamID2);
53
54
55
56
  int vlistID3 = vlistDuplicate(vlistID2);

  // taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = vlistInqTaxis(vlistID2);
57
  int taxisID3 = taxisDuplicate(taxisID2);
58

59
  int gridID1 = vlistInqVarGrid(vlistID1, 0);
60
  int gridID2 = vlistInqVarGrid(vlistID2, 0);
61
62
63

  size_t gridsize = vlistGridsizeMax(vlistID1);
  if (gridsize != vlistGridsizeMax(vlistID2))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
64
    cdoAbort("Gridsize of input files does not match! %zu and %zu", gridsize, vlistGridsizeMax(vlistID2));
65

Uwe Schulzweida's avatar
Uwe Schulzweida committed
66
  if (vlistNgrids(vlistID2) > 1 || vlistNgrids(vlistID1) > 1) cdoAbort("Too many different grids in input!");
67

Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
  int nvars = vlistNvars(vlistID1) == vlistNvars(vlistID2) ? vlistNvars(vlistID1) : -1;
69
70
71
72
  int nlevs = zaxisInqSize(vlistInqVarZaxis(vlistID1, 0));

  if (gridID1 != gridID2) cdoCompareGrids(gridID1, gridID2);

73
  strcpy(oname, cdoGetObase());
74
  int nchars = strlen(oname);
75

76
  const char *refname = cdoGetObase();
Cedrick Ansorge's avatar
Cedrick Ansorge committed
77
  filesuffix[0] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
  cdoGenFileSuffix(filesuffix, sizeof(filesuffix), pstreamInqFiletype(streamID1), vlistID1, refname);
79

Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
81
  Field ***eof = (Field ***) Malloc(nvars * sizeof(Field **));
  for (varID = 0; varID < nvars; varID++) eof[varID] = (Field **) Malloc(nlevs * sizeof(Field *));
82

83
  int eofID = 0;
84
85
86
87
88
89
90
91
92
93
  while (1)
    {
      nrecs = cdoStreamInqTimestep(streamID1, eofID);
      if (nrecs == 0) break;

      for (int recID = 0; recID < nrecs; recID++)
        {
          pstreamInqRecord(streamID1, &varID, &levelID);
          missval1 = vlistInqVarMissval(vlistID1, varID);
          if (eofID == 0)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94
            eof[varID][levelID] = (Field *) Malloc(1 * sizeof(Field));
95
          else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96
            eof[varID][levelID] = (Field *) Realloc(eof[varID][levelID], (eofID + 1) * sizeof(Field));
97
98
99
          eof[varID][levelID][eofID].grid = gridID1;
          eof[varID][levelID][eofID].nmiss = 0;
          eof[varID][levelID][eofID].missval = missval1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
          eof[varID][levelID][eofID].ptr = (double *) Malloc(gridsize * sizeof(double));
101
102
103
104
105
106
107
108
109
110
111
112
113
          arrayFill(gridsize, &eof[varID][levelID][eofID].ptr[0], missval1);

          if (varID >= nvars) cdoAbort("Internal error - too high varID");
          if (levelID >= nlevs) cdoAbort("Internal error - too high levelID");

          pstreamReadRecord(streamID1, eof[varID][levelID][eofID].ptr, &nmiss);
          eof[varID][levelID][eofID].nmiss = nmiss;
        }
      eofID++;
    }

  int neof = eofID;

114
  if (cdoVerbose) cdoPrint("%s contains %i eof's", cdoGetStreamName(0), neof);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
115
  // Create 1x1 Grid for output
116
  int gridID3 = gridCreate(GRID_LONLAT, 1);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
117
118
  gridDefXsize(gridID3, 1);
  gridDefYsize(gridID3, 1);
119
120
121
122
  double xvals = 0;
  double yvals = 0;
  gridDefXvals(gridID3, &xvals);
  gridDefYvals(gridID3, &yvals);
123

Cedrick Ansorge's avatar
Cedrick Ansorge committed
124
  // Create var-list and time-axis for output
125

126
  int ngrids = vlistNgrids(vlistID3);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
  for (i = 0; i < ngrids; i++) vlistChangeGridIndex(vlistID3, i, gridID3);
128

Cedrick Ansorge's avatar
Cedrick Ansorge committed
129
  vlistDefTaxis(vlistID3, taxisID3);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
  for (varID = 0; varID < nvars; varID++) vlistDefVarTimetype(vlistID3, varID, TIME_VARYING);
131

Cedrick Ansorge's avatar
Cedrick Ansorge committed
132
  // open streams for eofcoeff output
133
  std::vector<int> streamIDs(neof);
134
  for (eofID = 0; eofID < neof; eofID++)
Cedrick Ansorge's avatar
Cedrick Ansorge committed
135
    {
136
      oname[nchars] = '\0';
137

Cedrick Ansorge's avatar
Cedrick Ansorge committed
138
139
      sprintf(eof_name, "%5.5i", eofID);
      strcat(oname, eof_name);
140
141
      if (filesuffix[0]) strcat(oname, filesuffix);

142
      streamIDs[eofID] = cdoStreamOpenWrite(oname, cdoFiletype());
143

Uwe Schulzweida's avatar
Uwe Schulzweida committed
144
      if (cdoVerbose) cdoPrint("opened %s ('w')  as stream%i for %i. eof", oname, streamIDs[eofID], eofID + 1);
145

146
      pstreamDefVlist(streamIDs[eofID], vlistID3);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
147
    }
148

Cedrick Ansorge's avatar
Cedrick Ansorge committed
149
  // ALLOCATE temporary fields for data read and write
150
151
  in.ptr = (double *) Malloc(gridsize * sizeof(double));
  in.grid = gridID1;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
152
153
  out.missval = missval1;
  out.nmiss = 0;
154
155
  out.ptr = (double *) Malloc(1 * sizeof(double));

156
  int tsID = 0;
157
158
  while (1)
    {
159
      nrecs = cdoStreamInqTimestep(streamID2, tsID);
160
161
      if (nrecs == 0) break;

Cedrick Ansorge's avatar
Cedrick Ansorge committed
162
163
164
165
      taxisCopyTimestep(taxisID3, taxisID2);
      /*for ( eofID=0; eofID<neof; eofID++)
        {
          fprintf(stderr, "defining ts %i\n", tsID);
166
          pstreamDefTimestep(streamIDs[eofID],tsID);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
167
168
        }
      */
169
      for (int recID = 0; recID < nrecs; recID++)
Cedrick Ansorge's avatar
Cedrick Ansorge committed
170
        {
171
          pstreamInqRecord(streamID2, &varID, &levelID);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
172
          missval2 = vlistInqVarMissval(vlistID2, varID);
173
          pstreamReadRecord(streamID2, in.ptr, &nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
174
          in.nmiss = nmiss;
175
176

          for (eofID = 0; eofID < neof; eofID++)
Cedrick Ansorge's avatar
Cedrick Ansorge committed
177
            {
178
179
180
181
182
183
184
185
              if (recID == 0) pstreamDefTimestep(streamIDs[eofID], tsID);
              // if ( recID == 0 ) fprintf(stderr, "ts%i rec%i eof%i\n", tsID,
              // recID, eofID);
              out.ptr[0] = 0;
              out.grid = gridID3;
              out.missval = missval2;
              for (size_t i = 0; i < gridsize; i++)
                {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
                  if (!DBL_IS_EQUAL(in.ptr[i], missval2) && !DBL_IS_EQUAL(eof[varID][levelID][eofID].ptr[i], missval1))
Cedrick Ansorge's avatar
Cedrick Ansorge committed
187
                    {
188
189
                      // double tmp =
                      // w[i]*in.ptr[i]*eof[varID][levelID][eofID].ptr[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
                      double tmp = in.ptr[i] * eof[varID][levelID][eofID].ptr[i];
191
                      out.ptr[0] += tmp;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
192
                    }
193
194
195
196
197
198
199
200
201
                }
              if (!DBL_IS_EQUAL(out.ptr[0], 0.))
                nmiss = 0;
              else
                {
                  nmiss = 1;
                  out.ptr[0] = missval2;
                }

202
              pstreamDefRecord(streamIDs[eofID], varID, levelID);
203
204
205
              // fprintf(stderr, "%d %d %d %d %d %g\n", streamIDs[eofID],tsID,
              // recID, varID, levelID,*out.ptr);
              pstreamWriteRecord(streamIDs[eofID], out.ptr, nmiss);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
206
            }
207
208
          if (varID >= nvars) cdoAbort("Internal error - too high varID");
          if (levelID >= nlevs) cdoAbort("Internal error - too high levelID");
Cedrick Ansorge's avatar
Cedrick Ansorge committed
209
        }
210

Cedrick Ansorge's avatar
Cedrick Ansorge committed
211
212
      tsID++;
    }
213

Uwe Schulzweida's avatar
Uwe Schulzweida committed
214
  for (eofID = 0; eofID < neof; eofID++) pstreamClose(streamIDs[eofID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
215

216
217
  pstreamClose(streamID2);
  pstreamClose(streamID1);
218

Cedrick Ansorge's avatar
Cedrick Ansorge committed
219
  cdoFinish();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
220

Uwe Schulzweida's avatar
Uwe Schulzweida committed
221
  return 0;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
222
}