Arithlat.c 3.76 KB
Newer Older
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

5
  Copyright (C) 2003-2010 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
  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.
*/

/*
   This module contains the following operators:

      Arithlat   mulcoslat       Multiply with cos(lat)
      Arithlat   divcoslat       Divide by cos(lat)
*/


Ralf Mueller's avatar
Ralf Mueller committed
26
#include <cdi.h>
27
28
29
30
31
32
33
34
35
36
37
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"


#ifndef  DEG2RAD
#define  DEG2RAD  (M_PI/180.)   /* conversion for deg to rad */
#endif

void *Arithlat(void *argument)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
  static const char *func = "Arithlat";
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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
  int operatorID;
  int operfunc;
  int streamID1, streamID2;
  int gridsize, gridtype;
  int gridID, gridID0 = -1;
  int nlon = 0, nlat = 0, i, j;
  int nrecs, recID;
  int tsID;
  int varID, levelID;
  int vlistID1, vlistID2;
  int taxisID1, taxisID2;
  int nmiss;
  double *scale = NULL;
  double *array = NULL;

  cdoInitialize(argument);

  cdoOperatorAdd("mulcoslat", func_mul, 0, NULL);
  cdoOperatorAdd("divcoslat", func_div, 0, NULL);

  operatorID = cdoOperatorID();
  operfunc = cdoOperatorFunc(operatorID);

  streamID1 = streamOpenRead(cdoStreamName(0));
  if ( streamID1 < 0 ) cdiError(streamID1, "Open failed on %s", cdoStreamName(0));

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

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

  streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
  if ( streamID2 < 0 ) cdiError(streamID2, "Open failed on %s", cdoStreamName(1));

  streamDefVlist(streamID2, vlistID2);

  gridsize = vlistGridsizeMax(vlistID1);

  array = (double *) malloc(gridsize*sizeof(double));

  tsID = 0;
  while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
    {
      taxisCopyTimestep(taxisID2, taxisID1);

      streamDefTimestep(streamID2, tsID);

      for ( recID = 0; recID < nrecs; recID++ )
	{
	  streamInqRecord(streamID1, &varID, &levelID);
	  streamReadRecord(streamID1, array, &nmiss);
	  
	  gridID = vlistInqVarGrid(vlistID1, varID);

	  if ( gridID != gridID0 )
	    {
	      gridtype = gridInqType(gridID);
	      if ( gridtype != GRID_LONLAT && gridtype != GRID_GAUSSIAN )
		{
		  if ( gridInqType(gridID) == GRID_GAUSSIAN_REDUCED )
		    cdoAbort("Gaussian reduced grid found. Use option -R to convert it to a regular grid!");
		  else
		    cdoAbort("LONLAT or GAUSSIAN grid not found!");
		}

	      gridsize = gridInqSize(gridID);
	      nlon = gridInqXsize(gridID);
	      nlat = gridInqYsize(gridID);

	      scale = (double *) realloc(scale, nlat*sizeof(double));
	      gridInqYvals(gridID, scale);

	      if ( operfunc == func_mul )
		for ( j = 0; j < nlat; j++ ) scale[j] = cos(scale[j]*DEG2RAD);
	      else
		for ( j = 0; j < nlat; j++ ) scale[j] = 1./cos(scale[j]*DEG2RAD);

	      if ( cdoVerbose ) for ( j = 0; j < nlat; j++ ) cdoPrint("coslat  %3d  %g", j+1, scale[j]);
		  
	      gridID0 = gridID;
	    }

	  for ( j = 0; j < nlat; j++ )
	    for ( i = 0; i < nlon; i++ )
	      array[i+j*nlon] *= scale[j];

	  streamDefRecord(streamID2, varID, levelID);
	  streamWriteRecord(streamID2, array, nmiss);
	}
      tsID++;
    }

  streamClose(streamID2);
  streamClose(streamID1);

  if ( array ) free(array);
  if ( scale ) free(scale);

  cdoFinish();

  return (0);
}