after_vertint.cc 17.9 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
#include <stdio.h>
#include <stdlib.h> /* exit */
#include <string.h>
#include <math.h>
5
#include <assert.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6

Uwe Schulzweida's avatar
Uwe Schulzweida committed
7
#include "compare.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
8
#include "constants.h"
9
#include "after_vertint.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
10
11
12
13
14
15

#define  SCALEHEIGHT     (-7000.)
#define  SCALESLP        (101325.0)

int Mars = 0;

16

17
void height2pressure(double *restrict phlev, const double *restrict hlev, long nphlev)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
18
19
20
21
{
  double exp_arg;
  double height;

22
  for ( long k = 0; k < nphlev; k++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
23
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
24
      height  = hlev[k];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
25
26
27
      /*
	unitsel == 1 : hlev[k] is given in meters
	unitsel == 2 : hlev[k] is given in kilometers
Uwe Schulzweida's avatar
Uwe Schulzweida committed
28
	height2pressure needs meters (MKSC-standard)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
29
30
      */

Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
      exp_arg = height / SCALEHEIGHT;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
32
33
34

      phlev[k] = SCALESLP * exp(exp_arg);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
36
37


38
void pressure2height(double *restrict hlev, const double *restrict plev, long nphlev)
39
{
40
  for ( long k = 0; k < nphlev; k++ )
41
42
43
    {
      hlev[k] = log(plev[k]/SCALESLP)*SCALEHEIGHT;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
44
}
45
46


47
void presh(double *restrict fullp, double * halfp, const double *restrict vct,
48
	   const double *restrict ps, long nhlev, long ngp)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
49
50
51
52
53
54
55
{
  if ( ps == NULL )
    {
      fprintf(stderr, "ps undefined!\n");
      exit(EXIT_FAILURE);
    }

56
57
58
  double zp, ze;
  double *halfpres = halfp;
  for ( long lh = 0; lh < nhlev; lh++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59
60
61
62
    {
      zp = vct[lh];
      ze = vct[lh+nhlev+1];

63
      for ( long i = 0; i < ngp; i++ ) halfpres[i] = zp + ze * ps[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
64
65
66
67
68

      halfpres += ngp;
    }
  memcpy(halfpres, ps, ngp*sizeof(double));

69
70
71
  if ( fullp )
    {
      halfpres = halfp;
72
      for ( long i = 0; i < ngp*nhlev; i++ )
73
74
	fullp[i] = 0.5 * (halfpres[i] + halfpres[i+ngp]);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
76
77
78

} /* presh */


79
void genind(int *nx, const double *restrict plev, const double *restrict fullp, long ngp, long nplev, long nhlev)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
81
82
{
  memset(nx, 0, ngp*nplev*sizeof(int));

Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
#if defined(_OPENMP)
84
#pragma omp parallel for default(none) shared(nx,plev,fullp,ngp,nplev,nhlev)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
#endif
86
  for ( long lp = 0; lp < nplev; lp++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87
    {
88
89
      const double pres = plev[lp];
      int *restrict nxl = nx + lp*ngp;
90
91
92
93
94
95
96
97
      for ( long lh = 0; lh < nhlev; lh++ )
        {
          const double *restrict fullpx = fullp + lh*ngp;
          for ( long i = 0; i < ngp ; i++ )
            {
              if ( pres > fullpx[i] ) nxl[i] = lh;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
98
99
100
101
102
    }

}  /* genind */


103
void genindmiss(int *nx, const double *restrict plev, int ngp, int nplev, const double *restrict ps_prog, int *restrict pnmiss)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
105
#if defined(_OPENMP)
106
#pragma omp parallel for default(none) shared(nx,plev,ngp,nplev,ps_prog,pnmiss)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
#endif
108
  for ( long lp = 0; lp < nplev; lp++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
109
110
    {
      pnmiss[lp] = 0;
111
112
113
      const double pres = plev[lp];
      int *restrict nxl  = nx + lp*ngp;
      for ( long i = 0; i < ngp; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
114
115
116
117
118
119
120
121
122
123
124
125
	{
	  if ( pres > ps_prog[i] )
	    {
	      nxl[i] = -1;
	      pnmiss[lp]++;
	    }
	}
    }

}  /* genindmiss */


126
127
void extra_P(double *restrict slp, const double *restrict halfp, const double *restrict fullp,
	     const double *restrict geop, const double *restrict temp, long ngp)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128
129
{
  double alpha, tstar, tmsl, zprt, zprtal;
130
131
  const double zlapse = 0.0065;
  const double zrg = 1.0 / PlanetGrav;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
132

133
  for ( long j = 0; j < ngp; ++j )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
134
135
136
137
    {
      if ( geop[j] < 0.0001 && geop[j] > -0.0001 ) slp[j] = halfp[j];
      else
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
	  alpha = PlanetRD * zlapse * zrg;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139
140
141
142
143
144
145
146
147
148
149
150
151
152
	  tstar = (1.0 + alpha * (halfp[j]/fullp[j] - 1.0)) * temp[j];

	  if ( tstar < 255.0 ) tstar = 0.5 * (255.0 + tstar);

	  tmsl = tstar + zlapse * zrg * geop[j];
	  if ( tmsl > 290.5 && tstar > 290.5 )
	    {
	      tstar = 0.5 * (290.5 + tstar);
	      tmsl  = tstar;
	    }

	  if ( tmsl-tstar < 0.000001 && tstar-tmsl < 0.000001 )
	    alpha = 0.0;
	  else if ( geop[j] > 0.0001 || geop[j] < -0.0001 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153
	    alpha = PlanetRD * (tmsl-tstar) / geop[j];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
154

Uwe Schulzweida's avatar
Uwe Schulzweida committed
155
	  zprt   = geop[j] / (PlanetRD * tstar);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
156
157
158
159
160
161
162
	  zprtal = zprt * alpha;
	  slp[j] = halfp[j] * exp(zprt*(1.0-zprtal*(0.5-zprtal/3.0)));
	}
    }
}  /* extrap */


Uwe Schulzweida's avatar
Uwe Schulzweida committed
163
164
static 
double extra_T(double pres, double halfp, double fullp, double geop, double temp)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
{
166
167
168
169
170
171
  double peval;
  const double zlapse = 0.0065;
  const double zrg   = 1.0 / PlanetGrav;
  double tstar = (1.0 + zlapse * PlanetRD * zrg * (halfp/fullp - 1.0)) * temp;
  const double ztsz  = tstar;
  const double z1    = tstar + zlapse * zrg * geop;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
173
174

  if ( tstar < 255.0 ) tstar = 0.5 * (255.0 + tstar);

175
  double ztmsl = tstar + zlapse * zrg * geop;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
176
177
178
179
180
181
182
183
184
185
186
187
188

  if ( ztmsl > 290.5 && tstar > 290.5 )
    {
      tstar = 0.5 * (290.5 + tstar);
      ztmsl = tstar;
    }

  if ( ztmsl > 290.5 && tstar <= 290.5 ) ztmsl=290.5;

  if ( pres <= halfp )
    peval = ((halfp-pres)*temp+ (pres-fullp)*tstar)/ (halfp-fullp);
  else
    {
189
      double ztmsl = z1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
      tstar = ztsz;
191
      const double zhts  = geop * zrg;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
192
193
194
195
196
197
198

      if ( zhts > 2000. && z1 > 298. )
	{
	  ztmsl = 298.;
	  if ( zhts < 2500. ) ztmsl = 0.002*((2500.-zhts)*z1+(zhts-2000.)*ztmsl);
	}

199
      double zalph;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200
201
202
      if ( (ztmsl-tstar) < 0.000001 )
	zalph = 0.;
      else if (geop > 0.0001 || geop < -0.0001)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
	zalph = PlanetRD*(ztmsl-tstar)/geop;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
204
      else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
	zalph = PlanetRD*zlapse*zrg;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206

207
      double zalp  = zalph*log(pres/halfp);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
208
209
210
211
212
213
214
      peval = tstar*(1.0+zalp*(1.0+zalp*(0.5+0.16666666667*zalp)));
    }

  return peval;
}  /* extra_T */


Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
216
static 
double extra_Z(double pres, double halfp, double fullp, double geop, double temp)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
217
{
218
219
220
221
222
  const double zlapse = 0.0065;
  const double ztlim = 290.5;
  const double zrg   = 1.0 / PlanetGrav;
  double alpha = PlanetRD * zlapse * zrg;
  double tstar = (1.0 + alpha * (halfp/fullp - 1.0)) * temp;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
223
224
225

  if ( tstar < 255.0 ) tstar = 0.5 * (255.0 + tstar);

226
  double tmsl = tstar + zlapse * zrg * geop;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
227

228
  if ( tmsl > ztlim && tstar > ztlim )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
229
    {
230
      tstar = 0.5 * (ztlim + tstar);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231
232
233
      tmsl  = tstar;
    }

234
  if ( tmsl > ztlim && tstar <= ztlim ) tmsl = ztlim;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
235
236
237
238

  if ( tmsl-tstar < 0.000001 && tstar-tmsl < 0.000001 )
    alpha = 0.0;
  else if ( geop > 0.0001 || geop < -0.0001 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
239
    alpha = PlanetRD * (tmsl-tstar) / geop;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
240

241
242
  const double zalp   = log(pres/halfp);
  const double zalpal = zalp * alpha;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243

244
  return (geop - PlanetRD*tstar*zalp*(1.0 + zalpal*(0.5 + zalpal/6.0)))*zrg;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
245
246
247
}  /* extra_Z */


248
249
void interp_X(const double *restrict gt, double *pt, const double *restrict hyb_press, const int *nx,
	      const double *restrict plev, long nplev, long ngp, long nhlev, double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
250
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
251
#if defined(_OPENMP)
252
#pragma omp parallel for default(none) shared(gt,pt,hyb_press,nx,plev,nplev,ngp,nhlev,missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
253
#endif
254
  for ( long lp = 0; lp < nplev; lp++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
255
    {
256
257
258
259
260
      long nl, nh;
      const double pres = plev[lp];
      const int *restrict nxl  = nx + lp*ngp;
      double *restrict ptl  = pt + lp*ngp;
      for ( long i = 0; i < ngp; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
261
262
263
264
265
266
267
268
	{
	  if ( nxl[i] == -1 )
	    ptl[i] = missval;
	  else
	    {
	      nl = nxl[i] * ngp + i;
	      nh = nl + ngp;
	      if ( nh >= ngp*nhlev )
269
		ptl[i] =  gt[nl];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
270
	      else
271
		ptl[i] =  gt[nl] + (pres-hyb_press[nl])
272
273
		       * (gt[nh] - gt[nl])
                       / (hyb_press[nh] - hyb_press[nl]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
274
275
276
277
278
279
	    }
	}
    }
}  /* interp_X */


280
281
void interp_T(const double *restrict geop, const double *restrict gt, double *pt, const double *restrict fullp,
	      const double *restrict halfp, const int *nx, const double *restrict plev, long nplev, long ngp,
282
	      long nhlev, double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
283
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
284
#if defined(_OPENMP)
285
#pragma omp parallel for default(none) shared(geop,gt,pt,fullp,halfp,nx,plev,nplev,ngp,nhlev,missval,Mars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
286
#endif
287
  for ( long lp = 0; lp < nplev; lp++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
288
    {
289
290
291
292
      long nl, nh;
      const double pres = plev[lp];
      const int *restrict nxl  = nx + lp*ngp;
      double *restrict ptl  = pt + lp*ngp;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
293
#if defined(CRAY)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
294
295
#pragma _CRI inline extra_T
#endif
296
      for ( long i = 0; i < ngp; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
297
298
299
300
301
302
303
304
305
306
307
	{
	  nl = nxl[i];
	  if ( nl < 0 )
	    ptl[i] = missval;
	  else
	    {
	      if ( nl > nhlev-2 )
		{
		  if ( Mars )
		    ptl[i] = gt[(nhlev-1)*ngp+i];
		  else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
308
#if defined(SX)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
#pragma cdir inline
#endif
		    ptl[i] = extra_T(pres, halfp[nhlev*ngp+i],
				     fullp[(nhlev-1)*ngp+i], geop[i],
				     gt[(nhlev-1)*ngp+i]);
		}
	      else
		{
		  nh = nl + 1;
		  ptl[i] =  gt[nl*ngp+i] + (pres-fullp[nl*ngp+i])
                         * (gt[nh*ngp+i] - gt[nl*ngp+i])
                         / (fullp[nh*ngp+i] - fullp[nl*ngp+i]);
		}
	    }
	}
    }
}  /* interp_T */


328
329
void interp_Z(const double *restrict geop, const double *restrict gz, double *pz, const double *restrict fullp,
	      const double *restrict halfp, const int *nx, const double *restrict gt, const double *restrict plev,
330
	      long nplev, long ngp, long nhlev, double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
331
{
332
333
334
335
336
  assert(geop != NULL);
  assert(gz != NULL);
  assert(pz != NULL);
  assert(fullp != NULL);
  assert(halfp != NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337

Uwe Schulzweida's avatar
Uwe Schulzweida committed
338
#if defined(_OPENMP)
339
#pragma omp parallel for default(none) shared(geop,gz,pz,fullp,halfp,nx,gt,plev,nplev,ngp,nhlev,missval,Mars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
340
#endif
341
  for ( long lp = 0; lp < nplev; lp++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
342
    {
343
344
345
346
      long nl, nh;
      const double pres = plev[lp];
      const int *restrict nxl  = nx + lp*ngp;
      double *restrict pzl  = pz + lp*ngp;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
347
#if defined(CRAY)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
348
349
#pragma _CRI inline extra_Z
#endif
350
      for ( long i = 0; i < ngp; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
351
352
353
354
355
356
357
358
359
360
361
362
363
	{
	  nl = nxl[i];
	  if ( nl < 0 )
	    pzl[i] = missval;
	  else
	    {
	      if ( pres > halfp[(nl+1)*ngp+i] ) nl++;

	      if ( nl > nhlev-1 )
		{
		  if ( Mars )
		    pzl[i] = gt[(nhlev-1)*ngp+i];
		  else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
364
#if defined(SX)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
#pragma cdir inline
#endif
		    pzl[i] = extra_Z(pres, halfp[nhlev*ngp+i],
				     fullp[(nhlev-1)*ngp+i], geop[i],
				     gt[(nhlev-1)*ngp+i]);
		}
	      else
		{
		  nh = nl + 1;
		  pzl[i] =  gz[nl*ngp+i] + (pres-halfp[nl*ngp+i])
		         * (gz[nh*ngp+i] - gz[nl*ngp+i])
                         / (halfp[nh*ngp+i] - halfp[nl*ngp+i]);
		}
	    }
	}
    }
}  /* interp_Z */
382
383
384


/*
385
 * 3d vertical interpolation routine (see vert_interp_lev() in src/Intlevel.cc)
386
 */
387
void vert_interp_lev3d(int gridsize, double missval, double *vardata1, double *vardata2,
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
		       int nlev2, int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
  int i, ilev;
  int idx1, idx2;
  int offset;
  double wgt1, wgt2;
  double w1, w2;
  double var1L1, var1L2, *var2;

  for ( ilev = 0; ilev < nlev2; ilev++ )
    {
      offset = ilev*gridsize;
      var2 = vardata2 + offset;

      for ( i = 0; i < gridsize; i++ )
	{
          idx1 = lev_idx1[offset+i];
          idx2 = lev_idx2[offset+i];
          wgt1 = lev_wgt1[offset+i];
          wgt2 = lev_wgt2[offset+i];

          /* upper/lower values from input field */
          var1L1 = *(vardata1+idx1);
          var1L2 = *(vardata1+idx2);

          /* if (cdoVerbose) printf("i:%d level %d: idx1=%d idx2=%d (offset+i:%d) wgt1=%g wgt2=%g var1L1:%g var1L2:%g ",
           *                         i,       ilev, idx1,   idx2,    offset+i,    wgt1,   wgt2,   var1L1,   var1L2);
           */
	  w1 = wgt1;
	  w2 = wgt2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
418
419
	  if ( DBL_IS_EQUAL(var1L1, missval) ) w1 = 0;
	  if ( DBL_IS_EQUAL(var1L2, missval) ) w2 = 0;
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446

	  if ( IS_EQUAL(w1, 0) && IS_EQUAL(w2, 0) )
	    {
	      var2[i] = missval;
	    }
	  else if ( IS_EQUAL(w1, 0) )
	    {
	      if ( w2 >= 0.5 )
		var2[i] = var1L2;
	      else
		var2[i] = missval;	      
	    }
	  else if ( IS_EQUAL(w2, 0) )
	    {
	      if ( w1 >= 0.5 )
		var2[i] = var1L1;
	      else
		var2[i] = missval;	      
	    }
	  else
	    {
	      var2[i] = var1L1*w1 + var1L2*w2;
	    }
	}
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
447
448
#if defined(CDO)
#include "util.h"
449
450
451
452
453
454
455
/*
 * Create weights for the 3d vertical coordinate
 *
 * The resulting index sets lev_idx1 and lev_idx2 contain absolute numbers,i.e.
 * wrt. the given gridsize. They can directly be used to read values from 3d
 * data fields.
 *
456
 * 3d version of vert_gen_weights() (src/Intlevel.cc)
457
 */
458
void vert_gen_weights3d(bool expol, int nlev1, int gridsize, double *lev1, int nlev2, double *lev2,
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
			int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
  int i,i1, i2;
  int    idx1 = 0, idx2 = 0;
  double val1, val2 = 0;

  for ( i = 0; i < gridsize; i++ )
    {
      for ( i2 = 0; i2 < nlev2; i2++ )
        {
          /* Because 2 levels were added to the source vertical coordinate (one on
           * top, one at the bottom), its loop starts at 1 */
          for ( i1 = 1; i1 < nlev1; i1++ )
            {
              if ( lev1[(i1-1)*gridsize+i] < lev1[i1*gridsize+i] )
                {
                  idx1 = (i1-1)*gridsize+i;
                  idx2 = i1*gridsize+i;
                }
              else
                {
                  idx1 = i1*gridsize+i;
                  idx2 = (i1-1)*gridsize+i;
                }
              val1 = lev1[idx1];
              val2 = lev1[idx2];

              if ( lev2[i2*gridsize+i] > val1 && lev2[i2*gridsize+i] <= val2 ) break;
            }

          if ( i1 == nlev1 ) 
            {
              if ( expol )
                cdoAbort("Level %g at index %d not found! Use extrapolation", lev2[i2*gridsize],i2);
              else
                cdoAbort("Level %g at index %d not found!");
            }

          if ( i1-1 == 0 ) /* destination levels ios not covert by the first two input z levels */
            {
              lev_idx1[i2*gridsize+i] = gridsize+i;
              lev_idx2[i2*gridsize+i] = gridsize+i;
              lev_wgt1[i2*gridsize+i] = 0;
              if ( expol || IS_EQUAL(lev2[i2*gridsize+i], val2) )
                lev_wgt2[i2*gridsize+i] = 1;
              else
                lev_wgt2[i2*gridsize+i] = 0;
            }
          else if ( i1 == nlev1-1 ) /* destination level is beyond the last value of the input z field */
            {
              lev_idx1[i2*gridsize+i] = (nlev1-2)*gridsize+i;
              lev_idx2[i2*gridsize+i] = (nlev1-2)*gridsize+i;
              if ( expol || IS_EQUAL(lev2[i2*gridsize+i], val2) )
                lev_wgt1[i2*gridsize+i] = 1;
              else
                lev_wgt1[i2*gridsize+i] = 0;
              lev_wgt2[i2*gridsize+i] = 0;
            }
          else /* target z values has two bounday values in input z field */
            {
              lev_idx1[i2*gridsize+i] = idx1;
              lev_idx2[i2*gridsize+i] = idx2;
              lev_wgt1[i2*gridsize+i] = (lev1[idx2]        - lev2[i2*gridsize+i]) / (lev1[idx2] - lev1[idx1]);
              lev_wgt2[i2*gridsize+i] = (lev2[i2*gridsize+i] - lev1[idx1])        / (lev1[idx2] - lev1[idx1]);

            }
  /*         if (cdoVerbose)
   *         {
   *           printf("i:%d i2:%d\ti2*gridsize+i:%d\tlev2[i2*gridsize+i]:%g\tidx1:%d\tidx2:%d\tlev1[idx1]:%g\tlev1[idx2]:%g\t",
   *                   i, i2, i2*gridsize+i,         lev2[i2*gridsize+i],    idx1,    idx2,    lev1[idx1],    lev1[idx2]);
   *           printf("\tlev_wgt1:%g\tlev_wgt2:%g\n", lev_wgt1[i2*gridsize+i], lev_wgt2[i2*gridsize+i]);
   *         }
   */
          /* backshift of the indices because of the two additional levels in input vertical coordinate */
          lev_idx1[i2*gridsize+i] -= gridsize;
          lev_idx2[i2*gridsize+i] -= gridsize;

        }
    }
}


/*
 * Create weights for the 1d vertical coordinate from a 3d vertical coordinate
 *
 * The resulting index sets lev_idx1 and lev_idx2 contain absolute numbers,i.e.
 * wrt. the given gridsize. They can directly be used to read values from 3d
 * data fields.
 *
548
 * 3d1d version of vert_gen_weights() (src/Intlevel.cc)
549
 */
550
void vert_gen_weights3d1d(bool expol, int nlev1, int gridsize, double *lev1, int nlev2, double *lev2,
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
			  int *lev_idx1, int *lev_idx2, double *lev_wgt1, double *lev_wgt2)
{
  int i,i1, i2;
  int    idx1 = 0, idx2 = 0;
  double val1, val2 = 0;

  for ( i = 0; i < gridsize; i++ )
    {
      for ( i2 = 0; i2 < nlev2; i2++ )
        {
          /* Because 2 levels were added to the source vertical coordinate (one on
           * top, one at the bottom), its loop starts at 1 */
          for ( i1 = 1; i1 < nlev1; i1++ )
            {
              if ( lev1[(i1-1)*gridsize+i] < lev1[i1*gridsize+i] )
                {
                  idx1 = (i1-1)*gridsize+i;
                  idx2 = i1*gridsize+i;
                }
              else
                {
                  idx1 = i1*gridsize+i;
                  idx2 = (i1-1)*gridsize+i;
                }
              val1 = lev1[idx1];
              val2 = lev1[idx2];

              if ( lev2[i2] > val1 && lev2[i2] <= val2 ) break;
            }

          if ( i1 == nlev1 ) 
            {
              if ( expol )
                cdoAbort("Level %g at index %d not found! Use extrapolation", lev2[i2],i2);
              else
                cdoAbort("Level %g at index %d not found!");
            }

          if ( i1-1 == 0 ) /* destination levels ios not covert by the first two input z levels */
            {
              lev_idx1[i2*gridsize+i] = gridsize+i;
              lev_idx2[i2*gridsize+i] = gridsize+i;
              lev_wgt1[i2*gridsize+i] = 0;
              if ( expol || IS_EQUAL(lev2[i2], val2) )
                lev_wgt2[i2*gridsize+i] = 1;
              else
                lev_wgt2[i2*gridsize+i] = 0;
            }
          else if ( i1 == nlev1-1 ) /* destination level is beyond the last value of the input z field */
            {
              lev_idx1[i2*gridsize+i] = (nlev1-2)*gridsize+i;
              lev_idx2[i2*gridsize+i] = (nlev1-2)*gridsize+i;
              if ( expol || IS_EQUAL(lev2[i2], val2) )
                lev_wgt1[i2*gridsize+i] = 1;
              else
                lev_wgt1[i2*gridsize+i] = 0;
              lev_wgt2[i2*gridsize+i] = 0;
            }
          else /* target z values has two bounday values in input z field */
            {
              lev_idx1[i2*gridsize+i] = idx1;
              lev_idx2[i2*gridsize+i] = idx2;
              lev_wgt1[i2*gridsize+i] = (lev1[idx2]  - lev2[i2]) / (lev1[idx2] - lev1[idx1]);
              lev_wgt2[i2*gridsize+i] = (lev2[i2] - lev1[idx1])  / (lev1[idx2] - lev1[idx1]);

            }
  /*         if (cdoVerbose)
   *         {
   *           printf("i:%d i2:%d\ti2*gridsize+i:%d\tlev2[i2]:%g\tidx1:%d\tidx2:%d\tlev1[idx1]:%g\tlev1[idx2]:%g\t",
   *                   i, i2, i2*gridsize+i,         lev2[i2],    idx1,    idx2,    lev1[idx1],    lev1[idx2]);
   *           printf("\tlev_wgt1:%g\tlev_wgt2:%g\n", lev_wgt1[i2*gridsize+i], lev_wgt2[i2*gridsize+i]);
   *         }
   */
          /* backshift of the indices because of the two additional levels in input vertical coordinate */
          lev_idx1[i2*gridsize+i] -= gridsize;
          lev_idx2[i2*gridsize+i] -= gridsize;

        }
    }
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
631
#endif