Adisit.cc 10.3 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-2018 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
  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:

21
22
      Adisit      adisit          compute insitu from potential temperature
      Adisit      adipot          compute potential from insitu temperature
Uwe Schulzweida's avatar
Uwe Schulzweida committed
23
24
25
*/

#include <cdi.h>
Oliver Heidmann's avatar
Oliver Heidmann committed
26

Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
#include "cdo_int.h"
28
#include "pstream_int.h"
29
#include "util_string.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

/*
!>
!! transformation from potential to in situ temperature
!! according to Bryden, 1973, "New polynomials for thermal expansion,
!! adiabatic temperature gradient and potential temperature of sea
!! water". Deep Sea Research and Oceanographic Abstracts. 20, 401-408
!! (GILL P.602), which gives the inverse transformation for an
!! approximate value, all terms linear in t are taken after that one
!! newton step.  for the check value 8.4678516 the accuracy is 0.2
!! mikrokelvin.
!!
*/

/* compute insitu temperature from potential temperature */
45
46
static double
adisit_1(double tpot, double sal, double p)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
48
49
  double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9, a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
         a_c1 = 8.9309E-7, a_c2 = 3.1628E-8, a_c3 = 2.1987E-10, a_d = 4.1057E-9, a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50

51
52
53
54
  double qc = p * (a_a1 + p * (a_c1 - a_e1 * p));
  double qv = p * (a_b1 - a_d * p);
  double dc = 1. + p * (-a_a2 + p * (a_c2 - a_e2 * p));
  double dv = a_b2 * p;
55
56
  double qnq = -p * (-a_a3 + p * a_c3);
  double qn3 = -p * a_a4;
57
58

  double tpo = tpot;
59
60
61
62
63
64
  double qvs = qv * (sal - 35.) + qc;
  double dvs = dv * (sal - 35.) + dc;
  double t = (tpo + qvs) / dvs;
  double fne = -qvs + t * (dvs + t * (qnq + t * qn3)) - tpo;
  double fst = dvs + t * (2. * qnq + 3. * qn3 * t);
  t = t - fne / fst;
65
66

  return t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67
68
}

69
70
/* compute potential temperature from insitu temperature */
/* Ref: Gill, p. 602, Section A3.5:Potential Temperature */
71
72
static double
adipot(double t, double s, double p)
73
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
75
  double a_a1 = 3.6504E-4, a_a2 = 8.3198E-5, a_a3 = 5.4065E-7, a_a4 = 4.0274E-9, a_b1 = 1.7439E-5, a_b2 = 2.9778E-7,
         a_c1 = 8.9309E-7, a_c2 = 3.1628E-8, a_c3 = 2.1987E-10, a_d = 4.1057E-9, a_e1 = 1.6056E-10, a_e2 = 5.0484E-12;
76

77
  double s_rel = s - 35.0;
78

79
80
81
82
83
  double aa = (a_a1 + t * (a_a2 - t * (a_a3 - a_a4 * t)));
  double bb = s_rel * (a_b1 - a_b2 * t);
  double cc = (a_c1 + t * (-a_c2 + a_c3 * t));
  double cc1 = a_d * s_rel;
  double dd = (-a_e1 + a_e2 * t);
84

85
  double tpot = t - p * (aa + bb + p * (cc - cc1 + p * dd));
86

87
  return tpot;
88
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89

90
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
calc_adisit(long gridsize, long nlevel, double *pressure, Field tho, Field sao, Field tis)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92
{
93
94
95
  // pressure units: hPa
  // tho units:      Celsius
  // sao units:      psu
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96

97
  for (long levelID = 0; levelID < nlevel; ++levelID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
98
    {
99
      long offset = gridsize * levelID;
100
101
102
      double *thoptr = tho.ptr + offset;
      double *saoptr = sao.ptr + offset;
      double *tisptr = tis.ptr + offset;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103

104
105
      for (long i = 0; i < gridsize; ++i)
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
106
          if (DBL_IS_EQUAL(thoptr[i], tho.missval) || DBL_IS_EQUAL(saoptr[i], sao.missval))
107
108
109
110
111
112
113
114
            {
              tisptr[i] = tis.missval;
            }
          else
            {
              tisptr[i] = adisit_1(thoptr[i], saoptr[i], pressure[levelID]);
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
115
116
117
    }
}

118
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
calc_adipot(long gridsize, long nlevel, double *pressure, Field t, Field s, Field tpot)
120
{
121
122
  // pressure units: hPa
  // t units:        Celsius
Uwe Schulzweida's avatar
Uwe Schulzweida committed
123
  // s units:        psu
124

125
  for (long levelID = 0; levelID < nlevel; ++levelID)
126
    {
127
      long offset = gridsize * levelID;
128
129
130
      double *tptr = t.ptr + offset;
      double *sptr = s.ptr + offset;
      double *tpotptr = tpot.ptr + offset;
131

132
133
      for (long i = 0; i < gridsize; ++i)
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
134
          if (DBL_IS_EQUAL(tptr[i], t.missval) || DBL_IS_EQUAL(sptr[i], s.missval))
135
136
137
138
139
140
141
142
            {
              tpotptr[i] = tpot.missval;
            }
          else
            {
              tpotptr[i] = adipot(tptr[i], sptr[i], pressure[levelID]);
            }
        }
143
144
145
    }
}

146
147
void *
Adisit(void *process)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
148
149
{
  int nrecs;
150
  int varID, levelID;
151
  size_t offset;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152
  int i;
153
  size_t nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
154
  int thoID = -1, saoID = -1;
155
  char varname[CDI_MAX_NAME], stdname[CDI_MAX_NAME];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
156
157
158
  double pin = -1;
  double *single;

159
  cdoInitialize(process);
160
161
  int ADISIT = cdoOperatorAdd("adisit", 1, 1, "");
  int ADIPOT = cdoOperatorAdd("adipot", 1, 1, "");
162

163
164
  UNUSED(ADIPOT);

165
  int operatorID = cdoOperatorID();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
166

167
168
  if (operatorArgc() == 1) pin = parameter2double(operatorArgv()[0]);

169
  int streamID1 = cdoStreamOpenRead(cdoStreamName(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
170

171
  int vlistID1 = cdoStreamInqVlist(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172

173
  int nvars = vlistNvars(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
174

175
  for (varID = 0; varID < nvars; varID++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
176
    {
177
      int code = vlistInqVarCode(vlistID1, varID);
178

179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
      if (code <= 0)
        {
          vlistInqVarName(vlistID1, varID, varname);
          vlistInqVarStdname(vlistID1, varID, stdname);
          strtolower(varname);

          if (strcmp(varname, "s") == 0)
            code = 5;
          else if (strcmp(varname, "t") == 0)
            code = 2;
          else if (strcmp(stdname, "sea_water_salinity") == 0)
            code = 5;

          if (operatorID == ADISIT)
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
194
              if (strcmp(stdname, "sea_water_potential_temperature") == 0) code = 2;
195
196
197
198
199
200
201
202
203
204
205
            }
          else
            {
              if (strcmp(stdname, "sea_water_temperature") == 0) code = 2;
            }
        }

      if (code == 2)
        thoID = varID;
      else if (code == 5)
        saoID = varID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206
207
    }

208
209
  if (saoID == -1) cdoAbort("Sea water salinity not found!");
  if (thoID == -1) cdoAbort("Potential or Insitu temperature not found!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
210

211
  int gridID = vlistGrid(vlistID1, 0);
212
  size_t gridsize = vlist_check_gridsize(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
213

214
215
  int zaxisID = vlistInqVarZaxis(vlistID1, saoID);
  int nlevel1 = zaxisInqSize(zaxisID);
216
  zaxisID = vlistInqVarZaxis(vlistID1, thoID);
217
  int nlevel2 = zaxisInqSize(zaxisID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218

Uwe Schulzweida's avatar
Uwe Schulzweida committed
219
  if (nlevel1 != nlevel2) cdoAbort("temperature and salinity have different number of levels!");
220
  int nlevel = nlevel1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
221

222
223
  std::vector<double> pressure(nlevel);
  cdoZaxisInqLevels(zaxisID, pressure.data());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
224

225
  if (pin >= 0)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
226
    for (i = 0; i < nlevel; ++i) pressure[i] = pin;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
227
  else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
228
    for (i = 0; i < nlevel; ++i) pressure[i] /= 10;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
229

230
  if (cdoVerbose)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231
232
    {
      cdoPrint("Level Pressure");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
      for (i = 0; i < nlevel; ++i) cdoPrint("%5d  %g", i + 1, pressure[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
234
235
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
236
  Field tho, sao, tis;
237
238
239
  fieldInit(tho);
  fieldInit(sao);
  fieldInit(tis);
240
241
242
  tho.ptr = (double *) Malloc(gridsize * nlevel * sizeof(double));
  sao.ptr = (double *) Malloc(gridsize * nlevel * sizeof(double));
  tis.ptr = (double *) Malloc(gridsize * nlevel * sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
244
245
246

  tho.nmiss = 0;
  sao.nmiss = 0;
  tis.nmiss = 0;
247

Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
249
250
251
  tho.missval = vlistInqVarMissval(vlistID1, thoID);
  sao.missval = vlistInqVarMissval(vlistID1, saoID);
  tis.missval = tho.missval;

252
  int datatype = CDI_DATATYPE_FLT32;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
253
  if (vlistInqVarDatatype(vlistID1, thoID) == CDI_DATATYPE_FLT64 && vlistInqVarDatatype(vlistID1, saoID) == CDI_DATATYPE_FLT64)
254
    datatype = CDI_DATATYPE_FLT64;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
255

256
  int vlistID2 = vlistCreate();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
257

258
  int tisID2 = vlistDefVar(vlistID2, gridID, zaxisID, TIME_VARYING);
259
  if (operatorID == ADISIT)
260
261
262
263
264
265
266
267
268
269
270
271
272
    {
      vlistDefVarParam(vlistID2, tisID2, cdiEncodeParam(20, 255, 255));
      vlistDefVarName(vlistID2, tisID2, "to");
      vlistDefVarLongname(vlistID2, tisID2, "Sea water temperature");
      vlistDefVarStdname(vlistID2, tisID2, "sea_water_temperature");
    }
  else
    {
      vlistDefVarParam(vlistID2, tisID2, cdiEncodeParam(2, 255, 255));
      vlistDefVarName(vlistID2, tisID2, "tho");
      vlistDefVarLongname(vlistID2, tisID2, "Sea water potential temperature");
      vlistDefVarStdname(vlistID2, tisID2, "sea_water_potential_temperature");
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
273
274
  vlistDefVarUnits(vlistID2, tisID2, "K");
  vlistDefVarMissval(vlistID2, tisID2, tis.missval);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
275
  vlistDefVarDatatype(vlistID2, tisID2, datatype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
276

277
  int saoID2 = vlistDefVar(vlistID2, gridID, zaxisID, TIME_VARYING);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
278
  vlistDefVarParam(vlistID2, saoID2, cdiEncodeParam(5, 255, 255));
279
  vlistDefVarName(vlistID2, saoID2, "s");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
280
281
282
283
  vlistDefVarLongname(vlistID2, saoID2, "Sea water salinity");
  vlistDefVarStdname(vlistID2, saoID2, "sea_water_salinity");
  vlistDefVarUnits(vlistID2, saoID2, "psu");
  vlistDefVarMissval(vlistID2, saoID2, sao.missval);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
284
  vlistDefVarDatatype(vlistID2, saoID2, datatype);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
285

286
287
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
288
289
  vlistDefTaxis(vlistID2, taxisID2);

290
  int streamID2 = cdoStreamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
291

292
  pstreamDefVlist(streamID2, vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
293

294
  int tsID = 0;
295
  while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
296
297
    {
      taxisCopyTimestep(taxisID2, taxisID1);
298
      pstreamDefTimestep(streamID2, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
299

300
301
302
303
304
      for (int recID = 0; recID < nrecs; ++recID)
        {
          pstreamInqRecord(streamID1, &varID, &levelID);
          offset = gridsize * levelID;

305
306
          if (varID == thoID) pstreamReadRecord(streamID1, tho.ptr + offset, &tho.nmiss);
          if (varID == saoID) pstreamReadRecord(streamID1, sao.ptr + offset, &sao.nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
307
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
308

309
      if (operatorID == ADISIT)
310
        calc_adisit(gridsize, nlevel, pressure.data(), tho, sao, tis);
311
      else
312
        calc_adipot(gridsize, nlevel, pressure.data(), tho, sao, tis);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
313

314
315
316
      for (levelID = 0; levelID < nlevel; ++levelID)
        {
          offset = gridsize * levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
317

318
319
320
321
          single = tis.ptr + offset;
          nmiss = arrayNumMV(gridsize, single, tis.missval);
          pstreamDefRecord(streamID2, tisID2, levelID);
          pstreamWriteRecord(streamID2, single, nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
322

323
324
325
326
327
          single = sao.ptr + offset;
          nmiss = arrayNumMV(gridsize, single, sao.missval);
          pstreamDefRecord(streamID2, saoID2, levelID);
          pstreamWriteRecord(streamID2, single, nmiss);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
328
329
330
331

      tsID++;
    }

332
333
  pstreamClose(streamID2);
  pstreamClose(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
334
335
336

  vlistDestroy(vlistID2);

337
338
339
  Free(tis.ptr);
  Free(tho.ptr);
  Free(sao.ptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
340
341
342

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
343
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
344
}