Cloudlayer.cc 9.12 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.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2017 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
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.
*/


Ralf Mueller's avatar
Ralf Mueller committed
19
#include <cdi.h>
20
21
22
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
23
#include "after_vertint.h"
24
25
26
27
28
29
30
31
32
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
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


#define  SCALESLP        (101325.0)


/* ================================================= */
/* LayerCloud calculates random overlap cloud cover */
/* ================================================= */

static
void layer_cloud(const double *cc, double *ll, long MaxLev, long MinLev, long dimgp)
{
  long i, k;
  double maxval, minval;
  double ZEPSEC;

  ZEPSEC = 1.-1.0e-12;

  for ( i = 0; i < dimgp; i++ ) ll[i] = 1. - cc[i+MaxLev*dimgp];

  //  printf("maxlev %d minlev %d\n", MaxLev, MinLev);

  for ( k = MaxLev + 1; k <= MinLev; k++ )
    {
      for ( i = 0; i < dimgp; i++ )
	{
	  maxval = MAX(cc[i+(k-1)*dimgp], cc[i+k*dimgp]);
	  minval = MIN(cc[i+(k-1)*dimgp], ZEPSEC);
	  ll[i] *= (1. - maxval) / (1. - minval);
	}
    }

  for ( i = 0; i < dimgp; i++ ) ll[i] = 1. - ll[i];
}

static
void vct2plev(const double *vct, double *plevs, long nlevs)
{
  long k;

  for ( k = 0; k < nlevs; k++ )
    plevs[k] = vct[k] + vct[k+nlevs] * SCALESLP;
  /*
  for ( k = 0; k < nlevs; k++ )
    printf("plevs %ld %g\n", k, plevs[k]);
  
  for ( k = 1; k < nlevs; k++ )
    printf("plevs %ld %g\n", k-1, (plevs[k]+plevs[k-1])/2);
  */
}

static
void hl_index(int *kmax, int *kmin, double pmax, double pmin, long nhlevs, double *pph)
{
  long k;
  long MaxLev, MinLev;
   
  for ( k = 0; k < nhlevs; k++ )
    if ( pph[k] > pmax ) break;
  
  MaxLev = k - 1;
  
  for ( k  = nhlevs - 1; k >= 0; k-- )
    if ( pph[k] < pmin ) break;
   
  MinLev = k;

  *kmax = MaxLev;
  *kmin = MinLev;
}

static
void pl_index(int *kmax, int *kmin, double pmax, double pmin, long nlevs, double *plevs)
{
  long k;
  long MaxLev = -1, MinLev = -1;
   
  for ( k = 0; k < nlevs; k++ )
    if ( plevs[k] >= pmax )
      {
	MaxLev = k;
	break;
      }
  
  for ( k  = nlevs - 1; k >= 0; k-- )
    if ( plevs[k] < pmin )
      {
	MinLev = k;
	break;
      }
   
  *kmax = MaxLev;
  *kmin = MinLev;
}


#define NVARS  3

void *Cloudlayer(void *argument)
{
124
  int gridID, zaxisID;
125
  int nlevel, nlevs, nrecs, code;
126
  int varID, levelID;
127
  bool zrev = false;
128
  int i;
129
  int offset;
130
  size_t nmiss;
131
  int aclcacID = -1;
132
  int nvars2 = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
133
  int aclcac_code_found = 0;
134
  int kmin[NVARS], kmax[NVARS];
135
  char varname[CDI_MAX_NAME];
136
137
  double sfclevel = 0;
  double *cloud[NVARS];
138
  double pmin = 0, pmax = 0;
139
140
141

  cdoInitialize(argument);

142
143
144
145
  if ( operatorArgc() > 0 )
    {
      operatorCheckArgc(2);
      nvars2 = 1;
146
147
      pmin = parameter2double(operatorArgv()[0]);
      pmax = parameter2double(operatorArgv()[1]);
148
149
150
151
152
153
    }
  else
    {
      nvars2 = NVARS;
    }

154
  int streamID1 = pstreamOpenRead(cdoStreamName(0));
155

156
  int vlistID1 = pstreamInqVlist(streamID1);
157

158
  int gridsize = vlist_check_gridsize(vlistID1);
159

160
  int aclcac_code = 223;
161

162
  int nvars = vlistNvars(vlistID1);
163
164
165
166
167
168
169
170
171
172
173
174
  for ( varID = 0; varID < nvars; ++varID )
    {
      zaxisID  = vlistInqVarZaxis(vlistID1, varID);
      code = vlistInqVarCode(vlistID1, varID);

      if ( code <= 0 )
	{
	  vlistInqVarName(vlistID1, varID, varname);
	  strtolower(varname);
	  if ( strcmp(varname, "aclcac") == 0 ) code = 223;
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
175
176
177
178
179
180
181
182
183
      if  ( code == aclcac_code )
	{
	  aclcac_code_found = 1;
	  if ( zaxisInqType(zaxisID) == ZAXIS_PRESSURE || zaxisInqType(zaxisID) == ZAXIS_HYBRID )
	    {
	      aclcacID  = varID;
	      break;
	    }
	}
184
185
186
    }

  if ( aclcacID == -1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
187
188
189
190
191
192
    {
      if ( aclcac_code_found )
	cdoAbort("Cloud cover (parameter 223) not found on pressure or hybrid levels!");
      else
	cdoAbort("Cloud cover (parameter 223) not found!");
    }
193

194
  double missval = vlistInqVarMissval(vlistID1, aclcacID);
195
196
197
198
  gridID  = vlistInqVarGrid(vlistID1, aclcacID);
  zaxisID = vlistInqVarZaxis(vlistID1, aclcacID);

  nlevel = zaxisInqSize(zaxisID);
199
  int nhlev  = nlevel+1;
200

201
  double *aclcac = (double*) Malloc(gridsize*nlevel*sizeof(double));
202
  for ( varID = 0; varID < nvars2; ++varID )
203
    cloud[varID] = (double*) Malloc(gridsize*sizeof(double));
204
205
206

  if ( zaxisInqType(zaxisID) == ZAXIS_PRESSURE )
    {
207
      double *plevs = (double*) Malloc(nlevel*sizeof(double));
208
209
210
211
      zaxisInqLevels(zaxisID, plevs);
      if ( plevs[0] > plevs[nlevel-1] )
	{
	  double ptmp;
212
	  zrev = true;
213
214
215
216
217
218
219
220
221
222
223
224
225
	  for ( levelID = 0; levelID < nlevel/2; ++levelID )
	    {
	      ptmp = plevs[levelID];
	      plevs[levelID] = plevs[nlevel-1-levelID];
	      plevs[nlevel-1-levelID] = ptmp;
	    }
	}
      /*
      for ( levelID = 0; levelID < nlevel; ++levelID )
	{
	  printf("level %d %g\n", levelID, plevs[levelID]);
	}
      */
226
227
228
229
230
231
232
233
234
235
      if ( nvars2 == 1 )
	{
	  pl_index(&kmax[0], &kmin[0], pmin, pmax, nlevel, plevs);
	}
      else
	{
	  pl_index(&kmax[2], &kmin[2],  5000., 44000., nlevel, plevs);
	  pl_index(&kmax[1], &kmin[1], 46000., 73000., nlevel, plevs);
	  pl_index(&kmax[0], &kmin[0], 75000.,101300., nlevel, plevs);
	}
236

237
      Free(plevs);
238
239
240
    }
  else if ( zaxisInqType(zaxisID) == ZAXIS_HYBRID )
    {
241
      int nvct = zaxisInqVctSize(zaxisID);
242
243
      if ( nlevel == (nvct/2 - 1) )
	{
244
	  double *vct = (double*) Malloc(nvct*sizeof(double));
245
	  zaxisInqVct(zaxisID, vct);
246
247

	  nlevs = nlevel + 1;
248
	  double *plevs = (double*) Malloc(nlevs*sizeof(double));
249
	  vct2plev(vct, plevs, nlevs);
250
	  Free(vct);
251

252
253
254
255
256
257
258
259
260
261
	  if ( nvars2 == 1 )
	    {
	      hl_index(&kmax[0], &kmin[0], pmin, pmax, nhlev, plevs);
	    }
	  else
	    {
	      hl_index(&kmax[2], &kmin[2],  5000., 44000., nhlev, plevs);
	      hl_index(&kmax[1], &kmin[1], 46000., 73000., nhlev, plevs);
	      hl_index(&kmax[0], &kmin[0], 75000.,101300., nhlev, plevs);
	    }
262

263
	  Free(plevs);
264
265
266
267
268
269
270
271
 	}
      else
	cdoAbort("Unsupported vertical coordinate table format!");
   }
  else
    cdoAbort("Unsupported Z-Axis type!");


272
  int surfaceID = zaxisCreate(ZAXIS_SURFACE, 1);
273
274
  zaxisDefLevels(surfaceID, &sfclevel);

275
  int vlistID2 = vlistCreate();
276

277
278
  if ( nvars2 == 1 )
    {
279
      varID = vlistDefVar(vlistID2, gridID, surfaceID, TIME_VARYING);
280
      vlistDefVarParam(vlistID2, varID, cdiEncodeParam(33, 128, 255));
281
282
283
284
285
286
      vlistDefVarName(vlistID2, varID, "cld_lay");
      vlistDefVarLongname(vlistID2, varID, "cloud layer");
      vlistDefVarMissval(vlistID2, varID, missval);
    }
  else
    {
287
      varID = vlistDefVar(vlistID2, gridID, surfaceID, TIME_VARYING);
288
      vlistDefVarParam(vlistID2, varID, cdiEncodeParam(34, 128, 255));
289
290
291
292
      vlistDefVarName(vlistID2, varID, "low_cld");
      vlistDefVarLongname(vlistID2, varID, "low cloud");
      vlistDefVarMissval(vlistID2, varID, missval);

293
      varID = vlistDefVar(vlistID2, gridID, surfaceID, TIME_VARYING);
294
      vlistDefVarParam(vlistID2, varID, cdiEncodeParam(35, 128, 255));
295
296
297
298
      vlistDefVarName(vlistID2, varID, "mid_cld");
      vlistDefVarLongname(vlistID2, varID, "mid cloud");
      vlistDefVarMissval(vlistID2, varID, missval);

299
      varID = vlistDefVar(vlistID2, gridID, surfaceID, TIME_VARYING);
300
      vlistDefVarParam(vlistID2, varID, cdiEncodeParam(36, 128, 255));
301
302
303
304
      vlistDefVarName(vlistID2, varID, "hih_cld");
      vlistDefVarLongname(vlistID2, varID, "high cloud");
      vlistDefVarMissval(vlistID2, varID, missval);
    }
305

306
307
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
308
309
  vlistDefTaxis(vlistID2, taxisID2);

310
311
  int streamID2 = pstreamOpenWrite(cdoStreamName(1), cdoFiletype());
  pstreamDefVlist(streamID2, vlistID2);
312

313
  int tsID = 0;
314
  while ( (nrecs = pstreamInqTimestep(streamID1, tsID)) )
315
316
317
    {
      taxisCopyTimestep(taxisID2, taxisID1);

318
      pstreamDefTimestep(streamID2, tsID);
319
     
320
      for ( int recID = 0; recID < nrecs; recID++ )
321
	{
322
	  pstreamInqRecord(streamID1, &varID, &levelID);
323
324
325
326
327
328
329
330

	  if ( zrev )
	    offset = (nlevel-1-levelID)*gridsize;
	  else
	    offset = levelID*gridsize;

	  if ( varID == aclcacID )
	    {
331
	      pstreamReadRecord(streamID1, aclcac+offset, &nmiss);
Uwe Schulzweida's avatar
cleanup    
Uwe Schulzweida committed
332
	      if ( nmiss != 0 ) cdoAbort("Missing values unsupported!");
333
334
335
	    }
	}

336
      for ( varID = 0; varID < nvars2; ++varID )
337
338
339
340
	{
	  for ( i = 0; i < gridsize; i++ ) cloud[varID][i] = missval;
	}

341
      for ( varID = 0; varID < nvars2; ++varID )
342
343
344
345
346
	{
	  if ( kmax[varID] != -1 && kmin[varID] != -1 )
	    layer_cloud(aclcac, cloud[varID], kmax[varID], kmin[varID], gridsize);
	}

347
      for ( varID = 0; varID < nvars2; ++varID )
348
349
350
351
352
	{
	  nmiss = 0;
	  for ( i = 0; i < gridsize; i++ )
	    if ( DBL_IS_EQUAL(cloud[varID][i], missval) ) nmiss++;

353
354
	  pstreamDefRecord(streamID2, varID, 0);
	  pstreamWriteRecord(streamID2, cloud[varID], nmiss);
355
356
357
358
359
	}

      tsID++;
    }

360
361
  pstreamClose(streamID2);
  pstreamClose(streamID1);
362
363
364
 
  vlistDestroy(vlistID2);

365
  Free(aclcac);
366
  for ( varID = 0; varID < nvars2; ++varID )
367
    Free(cloud[varID]);
368
369
370

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
371
  return 0;
372
}