Arithdays.c 4.48 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2014 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
10
11
12
13
14
15
16
17
18
  See COPYING file for copying and redistribution conditions.

  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.

  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.
*/

/*
19
   This module contains the following operators:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
20

21
22
23
24
      Arithdays  muldpm          Multiply with days per month
      Arithdays  divdpm          Divide by days per month
      Arithdays  muldpy          Multiply with days per year
      Arithdays  divdpy          Divide by days per year
25
      Arithdays  muldoy          Multiply with day of year
26
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
28


Ralf Mueller's avatar
Ralf Mueller committed
29
#include <cdi.h>
30
31
32
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
33
34


35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
static
double dayofyear(int calendar, int vdate, int vtime)
{
  int month_360[12] = {30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30};
  int month_365[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
  int month_366[12] = {31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
  int *dpm;
  int im, dpy;
  int year, month, day;
  int hour, minute, second;
  double doy = 0;

  cdiDecodeDate(vdate, &year, &month, &day);
  cdiDecodeTime(vtime, &hour, &minute, &second);

  dpy = days_per_year(calendar, year);

  for ( im = 1; im < month; ++im )
    {
      if      ( dpy == 360 ) dpm = month_360;
      else if ( dpy == 365 ) dpm = month_365;
      else                   dpm = month_366;

      if ( im >= 1 && im <= 12 ) doy += dpm[im-1];
    }

  doy += (day-1);
  doy += (second+minute*60+hour*3600)/86400.;

  if ( cdoVerbose )
65
    cdoPrint("%d %d %d %g\n", vdate, vtime, dpy, doy);
66
67
68
69
70

  return (doy);
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
72
void *Arithdays(void *argument)
{
73
  int MULDOY;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
75
76
77
78
79
80
81
82
  int operatorID;
  int operfunc, operfunc2;
  int streamID1, streamID2;
  int gridsize;
  int nrecs, recID;
  int tsID;
  int varID, levelID;
  int vlistID1, vlistID2;
  int taxisID1, taxisID2;
83
  int vdate, vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
84
85
86
  int year, month, day;
  int calendar;
  double rconst;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87
  field_t field;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
88
89
90

  cdoInitialize(argument);

91
92
93
94
95
           cdoOperatorAdd("muldpm", func_mul, func_month, NULL);
           cdoOperatorAdd("divdpm", func_div, func_month, NULL);
           cdoOperatorAdd("muldpy", func_mul, func_year,  NULL);
           cdoOperatorAdd("divdpy", func_div, func_year,  NULL);
  MULDOY = cdoOperatorAdd("muldoy", func_mul,         0,  NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96
97

  operatorID = cdoOperatorID();
98
99
  operfunc = cdoOperatorF1(operatorID);
  operfunc2 = cdoOperatorF2(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

  streamID1 = streamOpenRead(cdoStreamName(0));

  vlistID1 = streamInqVlist(streamID1);
  vlistID2 = vlistDuplicate(vlistID1);

  taxisID1 = vlistInqTaxis(vlistID1);
  taxisID2 = taxisDuplicate(taxisID1);
  vlistDefTaxis(vlistID2, taxisID2);

  calendar = taxisInqCalendar(taxisID1);

  streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());

  streamDefVlist(streamID2, vlistID2);

  gridsize = vlistGridsizeMax(vlistID1);

118
  field_init(&field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
120
121
122
123
124
125
  field.ptr    = (double *) malloc(gridsize*sizeof(double));
  field.weight = NULL;

  tsID = 0;
  while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
    {
      vdate = taxisInqVdate(taxisID1);
126
      vtime = taxisInqVtime(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
128
129
130
131

      taxisCopyTimestep(taxisID2, taxisID1);

      streamDefTimestep(streamID2, tsID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
132
      cdiDecodeDate(vdate, &year, &month, &day);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
133

134
135
136
137
      if ( operatorID == MULDOY )
	{
	  rconst = dayofyear(calendar, vdate, vtime);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
      else
139
140
141
142
143
144
	{
	  if ( operfunc2 == func_month )
	    rconst = days_per_month(calendar, year, month);
	  else
	    rconst = days_per_year(calendar, year);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
145
146

      if ( cdoVerbose )
147
	cdoPrint("calendar %d  year %d  month %d  result %g", calendar, year, month, rconst);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173

      for ( recID = 0; recID < nrecs; recID++ )
	{
	  streamInqRecord(streamID1, &varID, &levelID);
	  streamReadRecord(streamID1, field.ptr, &field.nmiss);

	  field.grid    = vlistInqVarGrid(vlistID1, varID);
	  field.missval = vlistInqVarMissval(vlistID1, varID);

	  farcfun(&field, rconst, operfunc);

	  streamDefRecord(streamID2, varID, levelID);
	  streamWriteRecord(streamID2, field.ptr, field.nmiss);
	}
      tsID++;
    }

  streamClose(streamID2);
  streamClose(streamID1);

  if ( field.ptr ) free(field.ptr);

  cdoFinish();

  return (0);
}