← Back to branch summary

~javier-lopez/ubuntu/quantal/cdo/sru-1023329

« back to all changes in this revision

Viewing changes to src/Gridcell.c

  • Committer: Bazaar Package Importer
  • Author(s): Alastair McKinstry
  • Date: 2010-09-22 15:58:09 UTC
  • mfrom: (1.2.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20100922155809-u1d3vlmlqj02uxjt
Tags: 1.4.6.dfsg.1-1
New upstream release. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/Gridcell.c
2
2
  This file is part of CDO. CDO is a collection of Operators to
3
3
  manipulate and analyse Climate model Data.
4
4
 
5
 
  Copyright (C) 2003-2009 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
 
5
  Copyright (C) 2003-2010 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
6
6
  See COPYING file for copying and redistribution conditions.
7
7
 
8
8
  This program is free software; you can redistribute it and/or modify
20
20
 
21
21
      Gridcell   gridarea        Grid cell area in m^2
22
22
      Gridcell   gridweights     Grid cell weights
 
23
      Gridcell   gridmask        Grid mask
23
24
*/
24
25
 
25
26
 
26
 
#include "cdi.h"
 
27
#include <cdi.h>
27
28
#include "cdo.h"
28
29
#include "cdo_int.h"
29
30
#include "pstream.h"
31
32
 
32
33
void *Gridcell(void *argument)
33
34
{
34
 
  static char func[] = "Gridcell";
35
 
  int GRIDAREA, GRIDWGTS;
 
35
  static const char *func = "Gridcell";
 
36
  int GRIDAREA, GRIDWGTS, GRIDMASK;
36
37
  int operatorID;
37
38
  int streamID1, streamID2;
38
39
  int vlistID1, vlistID2;
39
40
  int gridID, zaxisID;
40
 
  int gridsize, gridtype;
 
41
  int gridtype;
41
42
  int status;
42
43
  int ngrids;
43
44
  int tsID, varID, levelID, taxisID;
44
 
  double *grid_area = NULL;
45
 
  double *grid_wgts = NULL;
46
 
  double *pdata;
 
45
  long i, gridsize;
47
46
  char *envstr;
 
47
  double *array = NULL;
48
48
  double  EarthRadius = 6371000; /* default radius of the earth in m */
49
49
  double PlanetRadius = EarthRadius;
50
50
 
51
51
  cdoInitialize(argument);
52
52
 
53
 
  envstr = getenv("PLANET_RADIUS");
54
 
  if ( envstr )
55
 
    {
56
 
      double fval;
57
 
      fval = atof(envstr);
58
 
      if ( fval > 0 )
59
 
        {
60
 
          PlanetRadius = fval;
61
 
          if ( cdoVerbose )
62
 
            cdoPrint("Set PlanetRadius to %g", PlanetRadius);
63
 
        }
64
 
    }
65
 
 
66
53
  GRIDAREA = cdoOperatorAdd("gridarea",     0,  0, NULL);
67
54
  GRIDWGTS = cdoOperatorAdd("gridweights",  0,  0, NULL);
 
55
  GRIDMASK = cdoOperatorAdd("gridmask",     0,  0, NULL);
68
56
 
69
57
  operatorID = cdoOperatorID();
70
58
 
 
59
  if ( operatorID == GRIDAREA || operatorID == GRIDWGTS )
 
60
    {
 
61
      envstr = getenv("PLANET_RADIUS");
 
62
      if ( envstr )
 
63
        {
 
64
          double fval;
 
65
          fval = atof(envstr);
 
66
          if ( fval > 0 )
 
67
            {
 
68
              PlanetRadius = fval;
 
69
              if ( cdoVerbose )
 
70
                cdoPrint("Set PlanetRadius to %g", PlanetRadius);
 
71
            }
 
72
        }
 
73
    }
 
74
 
71
75
  streamID1 = streamOpenRead(cdoStreamName(0));
72
76
  if ( streamID1 < 0 ) cdiError(streamID1, "Open failed on %s", cdoStreamName(0));
73
77
 
93
97
      vlistDefVarUnits(vlistID2, varID, "m2");
94
98
    }
95
99
  else if ( operatorID == GRIDWGTS )
96
 
    vlistDefVarName(vlistID2, varID, "cell_weights");
 
100
    {
 
101
      vlistDefVarName(vlistID2, varID, "cell_weights");
 
102
    }
 
103
  else if ( operatorID == GRIDMASK )
 
104
    {
 
105
      vlistDefVarName(vlistID2, varID, "grid_mask");
 
106
      vlistDefVarDatatype(vlistID2, varID, DATATYPE_UINT8);
 
107
    }
97
108
 
98
109
  taxisID = taxisCreate(TAXIS_ABSOLUTE);
99
110
  vlistDefTaxis(vlistID2, taxisID);
105
116
 
106
117
 
107
118
  gridsize = gridInqSize(gridID);
108
 
  grid_area = (double *) malloc(gridsize*sizeof(double));
109
 
  grid_wgts = (double *) malloc(gridsize*sizeof(double));
 
119
  array = (double *) malloc(gridsize*sizeof(double));
110
120
 
111
121
 
112
122
  tsID = 0;
120
130
  if ( operatorID == GRIDAREA )
121
131
    {
122
132
      gridtype = gridInqType(gridID);
123
 
      if ( gridtype != GRID_LONLAT      &&
124
 
           gridtype != GRID_GAUSSIAN    &&
125
 
           gridtype != GRID_LCC     &&
126
 
           gridtype != GRID_GME         &&
127
 
           gridtype != GRID_CURVILINEAR &&
128
 
           gridtype != GRID_CELL )
129
 
        {
130
 
          if ( gridInqType(gridID) == GRID_GAUSSIAN_REDUCED )
131
 
            cdoAbort("Unsupported grid type: %s, use CDO option -R to convert reduced to regular grid!",
132
 
                     gridNamePtr(gridtype));
133
 
          else
134
 
            cdoAbort("Unsupported grid type: %s", gridNamePtr(gridtype));
135
 
        }
136
 
      else
 
133
      if ( gridtype == GRID_LONLAT      ||
 
134
           gridtype == GRID_GAUSSIAN    ||
 
135
           gridtype == GRID_LCC         ||
 
136
           gridtype == GRID_GME         ||
 
137
           gridtype == GRID_CURVILINEAR ||
 
138
           gridtype == GRID_CELL )
137
139
        {
138
140
          if ( gridHasArea(gridID) )
139
141
            {
140
142
              if ( cdoVerbose ) cdoPrint("Using existing grid cell area!");
141
 
              gridInqArea(gridID, grid_area);
 
143
              gridInqArea(gridID, array);
142
144
            }
143
145
          else
144
146
            {
145
 
              int i;
146
 
 
147
 
              status = gridGenArea(gridID, grid_area);
 
147
              status = gridGenArea(gridID, array);
148
148
              if ( status == 1 )
149
149
                cdoAbort("Grid corner missing!");
150
150
              else if ( status == 2 )
151
151
                cdoAbort("Can't compute grid cell areas for this grid!");
152
152
 
153
153
              for ( i = 0; i < gridsize; ++i )
154
 
                grid_area[i] *= PlanetRadius*PlanetRadius;
 
154
                array[i] *= PlanetRadius*PlanetRadius;
155
155
            }
156
156
        }
157
 
 
158
 
      pdata = grid_area;
 
157
      else
 
158
        {
 
159
          if ( gridtype == GRID_GAUSSIAN_REDUCED )
 
160
            cdoAbort("Unsupported grid type: %s, use CDO option -R to convert reduced to regular grid!",
 
161
                     gridNamePtr(gridtype));
 
162
          else
 
163
            cdoAbort("Unsupported grid type: %s", gridNamePtr(gridtype));
 
164
        }
159
165
    }
160
 
  else
 
166
  else if ( operatorID == GRIDWGTS )
161
167
    {
162
 
      status = gridWeights(gridID, grid_wgts);
 
168
      status = gridWeights(gridID, array);
163
169
      if ( status != 0 )
164
170
          cdoWarning("Using constant grid cell area weights!");
165
 
 
166
 
      pdata = grid_wgts;
167
 
    }
168
 
 
169
 
  streamWriteRecord(streamID2, pdata, 0);
 
171
    }
 
172
  else if ( operatorID == GRIDMASK )
 
173
    {
 
174
      int *mask;
 
175
      mask = (int *) malloc(gridsize*sizeof(int));
 
176
      if ( gridInqMask(gridID, NULL) )
 
177
        {
 
178
          gridInqMask(gridID, mask);
 
179
        }
 
180
      else
 
181
        {
 
182
          for ( i = 0; i < gridsize; ++i ) mask[i] = 1;
 
183
        }
 
184
 
 
185
      for ( i = 0; i < gridsize; ++i ) array[i] = mask[i];
 
186
      free(mask);
 
187
    }
 
188
 
 
189
  streamWriteRecord(streamID2, array, 0);
170
190
 
171
191
 
172
192
  streamClose(streamID2);
173
193
  streamClose(streamID1);
174
194
 
175
 
  if ( grid_area ) free(grid_area);
176
 
  if ( grid_wgts ) free(grid_wgts);
 
195
  if ( array ) free(array);
177
196
 
178
197
  cdoFinish();
179
198