1 #include <algorithm>
2
3 #include "cdi.h"
4
5 #include "cimdOmp.h"
6
7 #define streamDefRecord cdo_def_record
8 #define streamWriteRecord cdo_write_record
9
10 #include "afterburner.h"
11 #include "constants.h"
12 #include "compare.h"
13 #include "vertical_interp.h"
14
15 int afterDebug = 0;
16
17 static char *
FieldName(int code,const char * text)18 FieldName(int code, const char *text)
19 {
20 static char name[256];
21
22 snprintf(name, sizeof(name), "[%3d].%s", code, text);
23
24 return name;
25 }
26
27 /* Free array space */
28 static void *
FreeMemory(void * ptr)29 FreeMemory(void *ptr)
30 {
31 Free(ptr);
32
33 return nullptr;
34 }
35
36 static void
FreeSpectral(struct Variable * vars)37 FreeSpectral(struct Variable *vars)
38 {
39 for (int code = MaxCodes - 1; code >= 0; --code)
40 if (vars[code].spectral) vars[code].spectral = (double *) FreeMemory(vars[code].spectral);
41 }
42
43 static void
FreeFourier(struct Variable * vars)44 FreeFourier(struct Variable *vars)
45 {
46 for (int code = 0; code < MaxCodes; code++)
47 if (vars[code].fourier) vars[code].fourier = (double *) FreeMemory(vars[code].fourier);
48 }
49
50 static void
FreeHybrid(struct Variable * vars)51 FreeHybrid(struct Variable *vars)
52 {
53 for (int code = 0; code < MaxCodes; code++)
54 if (vars[code].hybrid) vars[code].hybrid = (double *) FreeMemory(vars[code].hybrid);
55 }
56
57 static void
FreeGrid(struct Variable * vars)58 FreeGrid(struct Variable *vars)
59 {
60 for (int code = 0; code < MaxCodes; code++)
61 if (vars[code].grid) vars[code].grid = (double *) FreeMemory(vars[code].grid);
62 }
63
64 static void
FreeSamp(struct Variable * vars)65 FreeSamp(struct Variable *vars)
66 {
67 for (int code = 0; code < MaxCodes; code++)
68 if (vars[code].samp) vars[code].samp = (int *) FreeMemory(vars[code].samp);
69 }
70
71 /* alloc_dp - Allocate space for double array */
72 static double *
alloc_dp(int words,const char * array_name)73 alloc_dp(int words, const char *array_name)
74 {
75 double *result = nullptr;
76
77 if (words > 0)
78 {
79 result = (double *) Malloc(words * sizeof(double));
80 if (result == nullptr) cdo_sys_error(array_name, "No Memory!");
81 }
82
83 return result;
84 }
85
86 static void
IniQuaSum(double * dest,const double * restrict src,int len)87 IniQuaSum(double *dest, const double *restrict src, int len)
88 {
89 for (int i = 0; i < len; i++) dest[i] = src[i] * src[i];
90 }
91
92 static void
AddQuaSum(double * dest,const double * restrict src,int len)93 AddQuaSum(double *dest, const double *restrict src, int len)
94 {
95 for (int i = 0; i < len; i++) dest[i] += src[i] * src[i];
96 }
97
98 static void
VarQuaSum(double * Variance,const double * restrict Sum,int len,int n)99 VarQuaSum(double *Variance, const double *restrict Sum, int len, int n)
100 {
101 const double rn1 = 1.0 / (n - 1);
102
103 for (int i = 0; i < len; i++) Variance[i] = (Variance[i] - Sum[i] * Sum[i] * n) * rn1;
104
105 for (int i = 0; i < len; i++) Variance[i] = (Variance[i] > 0.0) ? std::sqrt(Variance[i]) : 0.0;
106 }
107
108 static void
AddVector(double * dest,const double * src,size_t len,size_t * nmiss,double missval)109 AddVector(double *dest, const double *src, size_t len, size_t *nmiss, double missval)
110 {
111 if (*nmiss)
112 {
113 array_add_array_mv(len, dest, src, missval);
114 *nmiss = array_num_mv(len, dest, missval);
115 if (*nmiss == 0) *nmiss = 1;
116 }
117 else
118 {
119 array_add_array(len, dest, src);
120 }
121 }
122
123 static void
Add2Vectors(double * dest,const double * restrict srcA,const double * restrict srcB,size_t len)124 Add2Vectors(double *dest, const double *restrict srcA, const double *restrict srcB, size_t len)
125 {
126 for (size_t i = 0; i < len; i++) dest[i] = srcA[i] + srcB[i];
127 }
128
129 static void
Sub2Vectors(double * dest,const double * restrict srcA,const double * restrict srcB,size_t len)130 Sub2Vectors(double *dest, const double *restrict srcA, const double *restrict srcB, size_t len)
131 {
132 for (size_t i = 0; i < len; i++) dest[i] = srcA[i] - srcB[i];
133 }
134
135 static void
MultVectorScalar(double * dest,const double * restrict src,double factor,size_t len,size_t nmiss,double missval)136 MultVectorScalar(double *dest, const double *restrict src, double factor, size_t len, size_t nmiss, double missval)
137 {
138 if (nmiss)
139 {
140 for (size_t i = 0; i < len; i++) dest[i] = (IS_EQUAL(src[i], missval)) ? missval : src[i] * factor;
141 }
142 else
143 {
144 for (size_t i = 0; i < len; i++) dest[i] = src[i] * factor;
145 }
146 }
147
148 static void
DivVectorIvector(double * dest,const double * restrict src,const int * samp,size_t len,size_t * nmiss,double missval)149 DivVectorIvector(double *dest, const double *restrict src, const int *samp, size_t len, size_t *nmiss, double missval)
150 {
151 *nmiss = 0;
152
153 for (size_t i = 0; i < len; i++)
154 {
155 if (IS_EQUAL(src[i], missval) || samp[i] == 0)
156 {
157 dest[i] = missval;
158 *nmiss = *nmiss + 1;
159 }
160 else
161 dest[i] = src[i] / samp[i];
162 }
163 }
164
165 void
after_read_vct(const char * vctfile,double ** vct,int * nvct)166 after_read_vct(const char *vctfile, double **vct, int *nvct)
167 {
168 int n;
169 char line[1024];
170 double va, vb;
171
172 FILE *fp = fopen(vctfile, "r");
173 if (fp == nullptr) cdo_sys_error("Open failed on %s", vctfile);
174
175 while (fgets(line, 1023, fp)) nvct++;
176
177 *nvct *= 2;
178 *vct = (double *) Malloc(*nvct * sizeof(double));
179
180 rewind(fp);
181 for (int i = 0; i < *nvct / 2; i++)
182 {
183 fgets(line, 1023, fp);
184 sscanf(line, "%d %lg %lg", &n, &va, &vb);
185 *vct[i] = va;
186 *vct[i + *nvct / 2] = vb;
187 }
188 fprintf(stdout, " Reading VCT for %d hybrid levels from file %s\n", *nvct / 2 - 1, vctfile);
189
190 fclose(fp);
191 }
192
193 void
after_gp2sp(const AfterControl & globs,struct Variable * vars,int ccode)194 after_gp2sp(const AfterControl &globs, struct Variable *vars, int ccode)
195 {
196 struct Variable *var = &vars[ccode];
197
198 if (var->spectral == nullptr)
199 {
200 if (var->hybrid == nullptr)
201 {
202 fprintf(stderr, "%d.hybrid not found\n", ccode);
203 exit(99);
204 }
205
206 if (var->fourier == nullptr)
207 {
208 long fieldSize = globs.DimFC * var->hlev;
209 var->fourier = alloc_dp(fieldSize, "gp2sp.fourier");
210 after_GP2FC(var->hybrid, var->fourier, globs.Latitudes, globs.Longitudes, var->hlev, globs.Fouriers);
211 }
212
213 var->spectral = alloc_dp(globs.Dim3SP, "gp2sp.spectral");
214 fc2sp(var->fourier, var->spectral, globs.pold, var->hlev, globs.Latitudes, globs.Fouriers, globs.Truncation);
215 }
216 }
217
218 void
after_GP2FC(double * gp,double * fc,long nlat,long nlon,long nlev,long nfc)219 after_GP2FC(double *gp, double *fc, long nlat, long nlon, long nlev, long nfc)
220 {
221 static long ifax[10];
222 static double *trig = nullptr;
223
224 if (ifax[9] != nlon)
225 {
226 if (trig) Free(trig);
227 trig = (double *) Malloc(nlon * sizeof(double));
228 int status = fft_set(trig, ifax, nlon);
229 if (status < 0) exit(1);
230 }
231
232 gp2fc(trig, ifax, gp, fc, nlat, nlon, nlev, nfc);
233 }
234
235 void
after_FC2GP(double * fc,double * gp,long nlat,long nlon,long nlev,long nfc)236 after_FC2GP(double *fc, double *gp, long nlat, long nlon, long nlev, long nfc)
237 {
238 static long ifax[10];
239 static double *trig = nullptr;
240
241 if (ifax[9] != nlon)
242 {
243 if (trig) Free(trig);
244 trig = (double *) Malloc(nlon * sizeof(double));
245 int status = fft_set(trig, ifax, nlon);
246 if (status < 0) exit(1);
247 }
248
249 fc2gp(trig, ifax, fc, gp, nlat, nlon, nlev, nfc);
250 }
251
252 // HUMTEST
253
254 static void
sh2rh(int AnalysisData,double * sphum,double * rhum,double * t,int lev,int dimgpout,const double * level,double * fullpresshybrid)255 sh2rh(int AnalysisData, double *sphum, double *rhum, double *t, int lev, int dimgpout, const double *level, double *fullpresshybrid)
256 {
257 int lp, i;
258 int lpi, lfp;
259 double es, qsatr;
260
261 /* ***************************************************** */
262 /* Define constants for calculation in presence of water */
263 /* ***************************************************** */
264 constexpr double RGAMW = (C_RCW - C_RCPV) / C_RV;
265 constexpr double RBETW = C_RLVTT / C_RV + RGAMW * C_RTT;
266 const double RALPW = std::log(C_RESTT) + RBETW / C_RTT + RGAMW * std::log(C_RTT);
267
268 /* ***************************************************** */
269 /* Define constants for calculation in presence of ice */
270 /* ***************************************************** */
271 /*
272 const double RGAMS = (C_RCS - C_RCPV) / C_RV;
273 const double RBETS = C_RLSTT / C_RV + RGAMS * C_RTT;
274 const double RALPS = std::log(C_RESTT) + RBETS / C_RTT + RGAMS * std::log(C_RTT);
275 */
276 const double *fullp = AnalysisData ? level : fullpresshybrid;
277
278 /***************************************************/
279 /* Diagnostics of saturation water vapour pressure */
280 /* over ice makes no sense, therefore ... */
281 /* Hint of Michael Ponater 08.10.97 */
282 /***************************************************/
283 constexpr double RGAM = RGAMW;
284 constexpr double RBET = RBETW;
285 const double RALP = RALPW;
286 for (lp = 0; lp < lev; lp++)
287 {
288 for (i = 0; i < dimgpout; i++)
289 {
290 lpi = lp * dimgpout + i;
291 lfp = (1 - AnalysisData) * lpi + AnalysisData * lp;
292 /*
293 if (t[lpi] < C_RTT) {
294 RGAM = RGAMS; RBET = RBETS; RALP = RALPS;
295 } else {
296 RGAM = RGAMW; RBET = RBETW; RALP = RALPW;
297 }
298 */
299 es = (std::exp(RALP - RBET / t[lpi] - RGAM * std::log(t[lpi]))) / fullp[lfp];
300 // qsat = es / (1. + C_RETV * (1. - es));
301 qsatr = (1. + C_RETV * (1. - es)) / es;
302 rhum[lpi] = sphum[lpi] * 100. * qsatr;
303 }
304 }
305 }
306
307 static void
rh2sh(double * sphum,double * rhum,double * t,int lev,int dimgpout,const double * level)308 rh2sh(double *sphum, double *rhum, double *t, int lev, int dimgpout, const double *level)
309 {
310 int lp, i;
311 int lpi;
312 double es, qsat;
313
314 /* ***************************************************** */
315 /* Define constants for calculation in presence of water */
316 /* ***************************************************** */
317 constexpr double RGAMW = (C_RCW - C_RCPV) / C_RV;
318 constexpr double RBETW = C_RLVTT / C_RV + RGAMW * C_RTT;
319 const double RALPW = std::log(C_RESTT) + RBETW / C_RTT + RGAMW * std::log(C_RTT);
320
321 /* ***************************************************** */
322 /* Define constants for calculation in presence of ice */
323 /* ***************************************************** */
324 // const double RGAMS = (C_RCS - C_RCPV) / C_RV;
325 // const double RBETS = C_RLSTT / C_RV + RGAMS * C_RTT;
326 // const double RALPS = std::log(C_RESTT) + RBETS / C_RTT + RGAMS * std::log(C_RTT);
327
328 /***************************************************/
329 /* Diagnostics of saturation water vapour pressure */
330 /* over ice makes no sense, therefore ... */
331 /* Hint of Michael Ponater 08.10.97 */
332 /***************************************************/
333
334 constexpr double RGAM = RGAMW;
335 constexpr double RBET = RBETW;
336 const double RALP = RALPW;
337 for (lp = 0; lp < lev; lp++)
338 {
339 for (i = 0; i < dimgpout; i++)
340 {
341 lpi = lp * dimgpout + i;
342 /* if (t[lpi] < C_RTT) { */
343 /* RGAM = RGAMS; RBET = RBETS; RALP = RALPS; */
344 /* } else { */
345 /* RGAM = RGAMW; RBET = RBETW; RALP = RALPW; */
346 /* } */
347 es = (std::exp(RALP - RBET / t[lpi] - RGAM * std::log(t[lpi]))) / level[lp];
348 qsat = es / (1. + C_RETV * (1. - es));
349 sphum[lpi] = rhum[lpi] * qsat / 100.;
350 }
351 }
352 }
353
354 void
after_FCrh2FCsh(const AfterControl & globs,struct Variable * vars)355 after_FCrh2FCsh(const AfterControl &globs, struct Variable *vars)
356 {
357 const long fieldSize = globs.DimGP * globs.NumLevelRequest;
358
359 if (vars[RHUMIDITY].grid == nullptr) vars[RHUMIDITY].grid = alloc_dp(fieldSize, "vars[RHUMIDITY].grid");
360 if (vars[TEMPERATURE].grid == nullptr) vars[TEMPERATURE].grid = alloc_dp(fieldSize, "vars[TEMPERATURE].grid");
361 if (vars[HUMIDITY].grid == nullptr) vars[HUMIDITY].grid = alloc_dp(fieldSize, "vars[HUMIDITY].grid");
362
363 after_FC2GP(vars[RHUMIDITY].fourier, vars[RHUMIDITY].grid, globs.Latitudes, globs.Longitudes, vars[RHUMIDITY].plev,
364 globs.Fouriers);
365 after_FC2GP(vars[TEMPERATURE].fourier, vars[TEMPERATURE].grid, globs.Latitudes, globs.Longitudes, vars[TEMPERATURE].plev,
366 globs.Fouriers);
367
368 rh2sh(vars[HUMIDITY].grid, vars[RHUMIDITY].grid, vars[TEMPERATURE].grid, globs.NumLevelRequest, globs.DimGP, globs.LevelRequest);
369
370 after_GP2FC(vars[HUMIDITY].grid, vars[HUMIDITY].fourier, globs.Latitudes, globs.Longitudes, vars[HUMIDITY].plev, globs.Fouriers);
371
372 vars[HUMIDITY].grid = (double *) FreeMemory(vars[HUMIDITY].grid);
373 vars[RHUMIDITY].grid = (double *) FreeMemory(vars[RHUMIDITY].grid);
374 vars[TEMPERATURE].grid = (double *) FreeMemory(vars[TEMPERATURE].grid);
375 }
376
377 void
after_SPuv2SPdv(const AfterControl & globs,struct Variable * vars)378 after_SPuv2SPdv(const AfterControl &globs, struct Variable *vars)
379 {
380 double *Div, *DivOut, *Vor, *VorOut;
381
382 Div = DivOut = vars[DIVERGENCE].spectral;
383 Vor = VorOut = vars[VORTICITY].spectral;
384 const long fieldSize = globs.DimFC * globs.NumLevelRequest;
385
386 if (vars[U_WIND].fourier == nullptr) vars[U_WIND].fourier = alloc_dp(fieldSize, "vars[U_WIND].fourier");
387 if (vars[V_WIND].fourier == nullptr) vars[V_WIND].fourier = alloc_dp(fieldSize, "vars[V_WIND].fourier");
388
389 sp2fc(vars[U_WIND].spectral, vars[U_WIND].fourier, globs.poli, globs.NumLevelRequest, globs.Latitudes, globs.Fouriers,
390 globs.Truncation);
391 sp2fc(vars[V_WIND].spectral, vars[V_WIND].fourier, globs.poli, globs.NumLevelRequest, globs.Latitudes, globs.Fouriers,
392 globs.Truncation);
393 uv2dv(vars[U_WIND].fourier, vars[V_WIND].fourier, Div, Vor, globs.pol2, globs.pol3, globs.NumLevelRequest, globs.Latitudes,
394 globs.Truncation);
395
396 vars[U_WIND].fourier = (double *) FreeMemory(vars[U_WIND].fourier);
397 vars[V_WIND].fourier = (double *) FreeMemory(vars[V_WIND].fourier);
398
399 for (int i = 0; i < globs.NumLevelRequest; ++i)
400 {
401 sp2sp(Div, globs.Truncation, DivOut, globs.Truncation);
402 sp2sp(Vor, globs.Truncation, VorOut, globs.Truncation);
403 Div += globs.DimSP;
404 Vor += globs.DimSP;
405 DivOut += globs.DimSP;
406 VorOut += globs.DimSP;
407 }
408 }
409
410 void
after_FCsh2FCrh(const AfterControl & globs,struct Variable * vars)411 after_FCsh2FCrh(const AfterControl &globs, struct Variable *vars)
412 {
413 const long fieldSize = globs.DimGP * globs.NumLevelRequest;
414
415 if (vars[RHUMIDITY].grid == nullptr) vars[RHUMIDITY].grid = alloc_dp(fieldSize, "vars[RHUMIDITY].grid");
416 if (vars[TEMPERATURE].grid == nullptr) vars[TEMPERATURE].grid = alloc_dp(fieldSize, "vars[TEMPERATURE].grid");
417 if (vars[HUMIDITY].grid == nullptr) vars[HUMIDITY].grid = alloc_dp(fieldSize, "vars[HUMIDITY].grid");
418
419 after_FC2GP(vars[HUMIDITY].fourier, vars[HUMIDITY].grid, globs.Latitudes, globs.Longitudes, vars[HUMIDITY].plev, globs.Fouriers);
420 after_FC2GP(vars[TEMPERATURE].fourier, vars[TEMPERATURE].grid, globs.Latitudes, globs.Longitudes, vars[TEMPERATURE].plev,
421 globs.Fouriers);
422
423 sh2rh(globs.AnalysisData, vars[HUMIDITY].grid, vars[RHUMIDITY].grid, vars[TEMPERATURE].grid, globs.NumLevelRequest, globs.DimGP,
424 globs.LevelRequest, vars[FULL_PRESS].hybrid);
425
426 after_GP2FC(vars[RHUMIDITY].grid, vars[RHUMIDITY].fourier, globs.Latitudes, globs.Longitudes, vars[RHUMIDITY].plev,
427 globs.Fouriers);
428
429 vars[HUMIDITY].grid = (double *) FreeMemory(vars[HUMIDITY].grid);
430 vars[RHUMIDITY].grid = (double *) FreeMemory(vars[RHUMIDITY].grid);
431 vars[TEMPERATURE].grid = (double *) FreeMemory(vars[TEMPERATURE].grid);
432 }
433 /* ENDE HUMTEST */
434
435 static void
CheckAnalyses(struct Variable * vars)436 CheckAnalyses(struct Variable *vars)
437 {
438 for (int code = 0; code < 272; code++)
439 if (vars[code].needed && code != DIVERGENCE && code != VORTICITY && code != STREAM && code != U_WIND && code != HUMIDITY
440 && code != VELOPOT && code != V_WIND && code != RHUMIDITY && code != GEOPOTHEIGHT && code != PS
441 && vars[code].spectral == nullptr && vars[code].grid == nullptr)
442 {
443 afterAbort("Code %3d not found", code);
444 }
445 }
446
447 /* Process Pressure Level data */
448 void
after_processPL(AfterControl & globs,struct Variable * vars)449 after_processPL(AfterControl &globs, struct Variable *vars)
450 {
451 int code, l;
452 long fieldSize;
453 int lindex, nlevel;
454 long offset;
455
456 globs.MeanCount++;
457 globs.TermCount++;
458
459 if (globs.MeanCount == 1)
460 {
461 if (globs.Debug) fprintf(stderr, "CheckAnalyses: %d %d\n", globs.TermCount, globs.MeanCount);
462 CheckAnalyses(vars);
463 globs.StartDate = globs.OldDate;
464 }
465 if (globs.TermCount > 120) globs.Debug = 0;
466
467 /* ============================== */
468 /* Computations in spectral space */
469 /* ============================== */
470
471 if (vars[TEMPERATURE].needed)
472 {
473 vars[TEMPERATURE].hlev = 2;
474 vars[TEMPERATURE].plev = globs.NumLevelRequest;
475 vars[TEMPERATURE].sfit = true;
476 }
477
478 if (vars[GEOPOTHEIGHT].comp)
479 {
480 vars[GEOPOTHEIGHT].hlev = 2;
481 vars[GEOPOTHEIGHT].plev = globs.NumLevelRequest;
482 vars[GEOPOTHEIGHT].sfit = true;
483 }
484
485 if (vars[GEOPOTHEIGHT].comp && vars[GEOPOTENTIAL].detected)
486 {
487 if (vars[GEOPOTHEIGHT].spectral == nullptr)
488 vars[GEOPOTHEIGHT].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "GEOPOTHEIGHT.spectral");
489 MultVectorScalar(vars[GEOPOTHEIGHT].spectral, vars[GEOPOTENTIAL].spectral, C_RG, globs.DimSP * globs.NumLevelRequest, 0, 0);
490 vars[GEOPOTENTIAL].needed = vars[GEOPOTENTIAL].selected;
491 }
492
493 if (globs.Type == 50 && vars[HUMIDITY].needed && vars[HUMIDITY].spectral == nullptr)
494 {
495 vars[HUMIDITY].plev = globs.NumLevelRequest;
496 vars[HUMIDITY].sfit = true;
497 vars[HUMIDITY].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[HUMIDITY].spectral");
498 /*
499 SPrh2SPsh();
500 */
501 vars[RHUMIDITY].needed = vars[RHUMIDITY].selected;
502 vars[TEMPERATURE].needed = vars[TEMPERATURE].selected;
503 }
504
505 if (vars[U_WIND].spectral && vars[V_WIND].spectral && (vars[DIVERGENCE].comp || vars[VORTICITY].comp))
506 {
507 vars[DIVERGENCE].hlev = vars[VORTICITY].hlev = 2;
508 vars[DIVERGENCE].plev = vars[VORTICITY].plev = globs.NumLevelRequest;
509 vars[DIVERGENCE].sfit = vars[VORTICITY].sfit = true;
510 if (vars[DIVERGENCE].spectral == nullptr)
511 vars[DIVERGENCE].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[DIVERGENCE].spectral");
512 if (vars[VORTICITY].spectral == nullptr)
513 vars[VORTICITY].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[VORTICITY].spectral");
514 after_SPuv2SPdv(globs, vars);
515 }
516
517 if (vars[U_WIND].comp || vars[V_WIND].comp)
518 {
519 vars[U_WIND].hlev = vars[V_WIND].hlev = 2;
520 vars[U_WIND].plev = vars[V_WIND].plev = globs.NumLevelRequest;
521 vars[U_WIND].sfit = vars[V_WIND].sfit = true;
522 if (vars[U_WIND].spectral == nullptr)
523 vars[U_WIND].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[U_WIND].spectral");
524 if (vars[V_WIND].spectral == nullptr)
525 vars[V_WIND].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[V_WIND].spectral");
526 dv2uv(vars[DIVERGENCE].spectral, vars[VORTICITY].spectral, vars[U_WIND].spectral, vars[V_WIND].spectral, globs.dv2uv_f1,
527 globs.dv2uv_f2, globs.Truncation, globs.DimSP, globs.NumLevelRequest);
528 }
529
530 if (vars[VELOPOT].comp)
531 {
532 vars[VELOPOT].hlev = 2;
533 vars[VELOPOT].plev = globs.NumLevelRequest;
534 vars[VELOPOT].sfit = true;
535 if (vars[VELOPOT].spectral == nullptr)
536 vars[VELOPOT].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[VELOPOT].spectral");
537 dv2ps(vars[DIVERGENCE].spectral, vars[VELOPOT].spectral, globs.NumLevelRequest, globs.Truncation);
538 }
539
540 if (vars[STREAM].comp)
541 {
542 vars[STREAM].hlev = 2;
543 vars[STREAM].plev = globs.NumLevelRequest;
544 vars[STREAM].sfit = true;
545 if (vars[STREAM].spectral == nullptr)
546 vars[STREAM].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[STREAM].spectral");
547 dv2ps(vars[VORTICITY].spectral, vars[STREAM].spectral, globs.NumLevelRequest, globs.Truncation);
548 }
549
550 /* --------------------------- */
551 /* Output of spectral fields */
552 /* --------------------------- */
553
554 if (globs.Type == 50)
555 {
556 for (code = 0; code < MaxCodes; code++)
557 if (vars[code].selected)
558 {
559 if (!vars[code].spectral) afterAbort("Code %d not available on spectral space!", code);
560
561 nlevel = zaxisInqSize(vars[code].ozaxisID);
562 for (lindex = 0; lindex < nlevel; lindex++)
563 {
564 offset = lindex * globs.DimSP;
565 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
566 streamWriteRecord(globs.ostreamID, vars[code].spectral + offset, 0);
567 }
568 }
569
570 FreeSpectral(vars);
571 return;
572 }
573
574 /* =============================== */
575 /* Transformation to fourier space */
576 /* Computations in fourier space */
577 /* =============================== */
578
579 if (globs.Type >= 60)
580 {
581 for (code = 0; code < MaxCodes; code++)
582 if (vars[code].needed && vars[code].spectral)
583 {
584 if (vars[code].fourier == nullptr)
585 {
586 fieldSize = vars[code].plev * globs.DimFC;
587 vars[code].fourier = alloc_dp(fieldSize, FieldName(code, "fourier"));
588 }
589 sp2fc(vars[code].spectral, vars[code].fourier, globs.poli, vars[code].plev, globs.Latitudes, globs.Fouriers,
590 globs.Truncation);
591 }
592 if (vars[U_WIND].needed && vars[U_WIND].fourier)
593 scaluv(vars[U_WIND].fourier, globs.rcoslat, globs.Latitudes, globs.Fouriers * globs.NumLevelRequest);
594 if (vars[V_WIND].needed && vars[V_WIND].fourier)
595 scaluv(vars[V_WIND].fourier, globs.rcoslat, globs.Latitudes, globs.Fouriers * globs.NumLevelRequest);
596
597 /* HUMTEST */
598 if (globs.Type < 70 && vars[HUMIDITY].needed && vars[HUMIDITY].fourier == nullptr)
599 {
600 vars[HUMIDITY].plev = globs.NumLevelRequest;
601 vars[HUMIDITY].sfit = true;
602 vars[HUMIDITY].fourier = alloc_dp(globs.DimFC * globs.NumLevelRequest, "vars[HUMIDITY].fourier");
603
604 after_FCrh2FCsh(globs, vars);
605
606 vars[RHUMIDITY].needed = vars[RHUMIDITY].selected;
607 vars[TEMPERATURE].needed = vars[TEMPERATURE].selected;
608 }
609
610 if (globs.Type < 70 && vars[RHUMIDITY].needed && vars[RHUMIDITY].fourier == nullptr)
611 {
612 vars[RHUMIDITY].plev = globs.NumLevelRequest;
613 vars[RHUMIDITY].sfit = true;
614 vars[RHUMIDITY].fourier = alloc_dp(globs.DimFC * globs.NumLevelRequest, "vars[RHUMIDITY].fourier");
615
616 after_FCsh2FCrh(globs, vars);
617
618 vars[HUMIDITY].needed = vars[HUMIDITY].selected;
619 vars[TEMPERATURE].needed = vars[TEMPERATURE].selected;
620 }
621 /* ENDE HUMTEST */
622 }
623
624 FreeSpectral(vars);
625
626 /* -------------------------- */
627 /* Output of fourier fields */
628 /* -------------------------- */
629
630 if (globs.Type == 60)
631 {
632 for (code = 0; code < MaxCodes; code++)
633 if (vars[code].selected)
634 {
635 if (!vars[code].fourier) afterAbort("Code %d not available on fourier space!", code);
636
637 nlevel = zaxisInqSize(vars[code].ozaxisID);
638 for (lindex = 0; lindex < nlevel; lindex++)
639 {
640 offset = lindex * globs.DimFC;
641 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
642 streamWriteRecord(globs.ostreamID, vars[code].fourier + offset, 0);
643 }
644 }
645
646 FreeFourier(vars);
647 return;
648 }
649
650 /* ----------------------- */
651 /* Output of zonal means */
652 /* ----------------------- */
653
654 if (globs.Type == 61)
655 {
656 for (code = 0; code < MaxCodes; code++)
657 if (vars[code].selected)
658 {
659 if (!vars[code].fourier) afterAbort("Code %d not available on zonal mean!", code);
660
661 nlevel = zaxisInqSize(vars[code].ozaxisID);
662 for (lindex = 0; lindex < nlevel; lindex++)
663 {
664 offset = lindex * globs.DimFC;
665 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
666 streamWriteRecord(globs.ostreamID, vars[code].fourier + offset, 0);
667 }
668 }
669
670 FreeFourier(vars);
671 return;
672 }
673
674 /* ============================ */
675 /* Transformation to gridpoints */
676 /* ============================ */
677
678 if (vars[PS].comp && vars[LNPS].grid)
679 {
680 if (vars[PS].grid == nullptr) vars[PS].grid = alloc_dp(globs.DimGP, "Ps");
681 for (l = 0; l < globs.DimGP; l++) vars[PS].grid[l] = std::exp(vars[LNPS].grid[l]);
682 }
683
684 if (globs.Type >= 70)
685 {
686 for (code = 0; code < MaxCodes; code++)
687 if (vars[code].needed && vars[code].fourier)
688 {
689 if (vars[code].grid == nullptr)
690 {
691 fieldSize = vars[code].plev * globs.DimGP;
692 vars[code].grid = alloc_dp(fieldSize, FieldName(code, "grid"));
693 }
694
695 after_FC2GP(vars[code].fourier, vars[code].grid, globs.Latitudes, globs.Longitudes, vars[code].plev, globs.Fouriers);
696 }
697 }
698
699 FreeFourier(vars);
700
701 /* HUMTEST */
702 /* -------------------------------- */
703 /* Computation in gridpoint space */
704 /* -------------------------------- */
705
706 if (vars[RHUMIDITY].needed && vars[RHUMIDITY].grid == nullptr)
707 {
708 vars[RHUMIDITY].plev = globs.NumLevelRequest;
709 vars[RHUMIDITY].sfit = true;
710 vars[RHUMIDITY].grid = alloc_dp(globs.DimGP * globs.NumLevelRequest, "vars[RHUMIDITY].grid");
711 sh2rh(globs.AnalysisData, vars[HUMIDITY].grid, vars[RHUMIDITY].grid, vars[TEMPERATURE].grid, globs.NumLevelRequest,
712 globs.DimGP, globs.LevelRequest, vars[FULL_PRESS].hybrid);
713 vars[HUMIDITY].needed = vars[HUMIDITY].selected;
714 vars[TEMPERATURE].needed = vars[TEMPERATURE].selected;
715 }
716
717 if (vars[HUMIDITY].needed && vars[HUMIDITY].grid == nullptr)
718 {
719 vars[HUMIDITY].plev = globs.NumLevelRequest;
720 vars[HUMIDITY].sfit = true;
721 vars[HUMIDITY].grid = alloc_dp(globs.DimGP * globs.NumLevelRequest, "vars[HUMIDITY].grid");
722 rh2sh(vars[HUMIDITY].grid, vars[RHUMIDITY].grid, vars[TEMPERATURE].grid, globs.NumLevelRequest, globs.DimGP,
723 globs.LevelRequest);
724 vars[RHUMIDITY].needed = vars[RHUMIDITY].selected;
725 vars[TEMPERATURE].needed = vars[TEMPERATURE].selected;
726 }
727 /* HUMTEST ENDE */
728
729 /* -------------------------- */
730 /* Computation of Means */
731 /* -------------------------- */
732
733 if (globs.Mean)
734 for (code = 0; code < MaxCodes; code++)
735 if (vars[code].needed && vars[code].grid)
736 {
737 fieldSize = globs.DimGP * vars[code].plev;
738 if (vars[code].mean == nullptr) vars[code].mean = alloc_dp(fieldSize, FieldName(code, "mean"));
739
740 if (globs.MeanCount == 1)
741 array_copy(fieldSize, vars[code].grid, vars[code].mean);
742 else
743 AddVector(vars[code].mean, vars[code].grid, fieldSize, &vars[code].nmiss, vars[code].missval);
744
745 if (globs.EndOfInterval)
746 {
747 if (vars[code].samp == nullptr)
748 MultVectorScalar(vars[code].mean, vars[code].mean, 1.0 / globs.MeanCount, fieldSize, vars[code].nmiss,
749 vars[code].missval);
750 else
751 DivVectorIvector(vars[code].mean, vars[code].mean, vars[code].samp, fieldSize, &vars[code].nmiss,
752 vars[code].missval);
753 }
754 }
755
756 /* -------------------------- */
757 /* Computation of Variances */
758 /* -------------------------- */
759
760 if (globs.Mean > 1)
761 for (code = 0; code < MaxCodes; code++)
762 if (vars[code].needed && vars[code].mean)
763 {
764 fieldSize = globs.DimGP * vars[code].plev;
765 if (vars[code].variance == nullptr) vars[code].variance = alloc_dp(fieldSize, FieldName(code, "var"));
766 if (globs.MeanCount == 1)
767 IniQuaSum(vars[code].variance, vars[code].grid, fieldSize);
768 else
769 AddQuaSum(vars[code].variance, vars[code].grid, fieldSize);
770
771 if (globs.EndOfInterval) VarQuaSum(vars[code].variance, vars[code].mean, fieldSize, globs.MeanCount);
772 }
773
774 if (globs.Mean && !globs.EndOfInterval)
775 {
776 FreeGrid(vars);
777 return;
778 }
779
780 /* ---------------------------------------------- */
781 /* Output of pressure level means and variances */
782 /* ---------------------------------------------- */
783
784 if (globs.Type == 70 && globs.Mean && globs.EndOfInterval)
785 {
786 for (code = 0; code < MaxCodes; code++)
787 if (vars[code].selected)
788 {
789 nlevel = zaxisInqSize(vars[code].ozaxisID);
790 for (lindex = 0; lindex < nlevel; lindex++)
791 {
792 offset = lindex * globs.DimGP;
793 if (globs.Mean != 2)
794 {
795 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
796 streamWriteRecord(globs.ostreamID, vars[code].mean + offset, vars[code].nmiss);
797 }
798 if (globs.Mean >= 2)
799 {
800 streamDefRecord(globs.ostreamID2, vars[code].ovarID2, lindex);
801 streamWriteRecord(globs.ostreamID2, vars[code].variance + offset, vars[code].nmiss);
802 }
803 }
804 }
805
806 FreeSamp(vars);
807 FreeGrid(vars);
808 return;
809 }
810
811 /* -------------------------------- */
812 /* Output of pressure level grids */
813 /* -------------------------------- */
814
815 if (globs.Type == 70)
816 {
817 for (code = 0; code < MaxCodes; code++)
818 if (vars[code].selected)
819 {
820 nlevel = zaxisInqSize(vars[code].ozaxisID);
821 for (lindex = 0; lindex < nlevel; lindex++)
822 {
823 offset = lindex * globs.DimGP;
824 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
825 streamWriteRecord(globs.ostreamID, vars[code].grid + offset, vars[code].nmiss);
826 }
827 }
828
829 FreeGrid(vars);
830 return;
831 }
832 }
833
834 static void
theta(double * pthetaf,double * pthetah,double * ph,double * ps,double * tf,double * ts,int levels,int dimgp,int dim3gp)835 theta(double *pthetaf, double *pthetah, double *ph, double *ps, double *tf, double *ts, int levels, int dimgp, int dim3gp)
836 {
837 double *thetah = pthetah;
838 double *thetaf = pthetaf;
839
840 const double kappa = PlanetRD / C_RCPD;
841
842 for (int h = 0; h < dimgp; h++) thetah[h] = 0.0;
843 thetah += dimgp;
844 for (int l = 0; l < levels - 1; l++)
845 {
846 for (int h = 0; h < dimgp; h++) thetah[h] = 0.5 * (tf[h] + tf[h + dimgp]) * std::pow((ps[h] / ph[h]), kappa);
847
848 ph += dimgp;
849 tf += dimgp;
850 thetah += dimgp;
851 }
852
853 array_copy(dimgp, ts, thetah);
854 thetah = pthetah;
855 for (int h = 0; h < dim3gp; h++) thetaf[h] = 0.5 * (thetah[h] + thetah[h + dimgp]);
856 }
857
858 static void
windSpeed(double * uvspeed,double * u,double * v,int dim3gp)859 windSpeed(double *uvspeed, double *u, double *v, int dim3gp)
860 {
861 for (int i = 0; i < dim3gp; i++) uvspeed[i] = std::sqrt(u[i] * u[i] + v[i] * v[i]);
862 }
863
864 static void
Omega(double * omega_in,double * divergence,double * u_wind,double * v_wind,double * halfpress,double * fullpress,double * dpsdx,double * dpsdy,double * vct,int dimgp,int nlev)865 Omega(double *omega_in, double *divergence, double *u_wind, double *v_wind, double *halfpress, double *fullpress, double *dpsdx,
866 double *dpsdy, double *vct, int dimgp, int nlev)
867 {
868 double DeltaHybrid, Cterm, Pterm;
869 double *diver, *halfp, *fullp, *uwind, *vwind;
870 double *omega = omega_in;
871
872 // Compute half level part of vertical velocity
873
874 for (int i = 0; i < dimgp; i++) omega[i] = 0.0;
875
876 for (int j = 0; j < nlev; j++)
877 {
878 omega = omega_in + j * dimgp;
879 halfp = halfpress + j * dimgp;
880 diver = divergence + j * dimgp;
881 uwind = u_wind + j * dimgp;
882 vwind = v_wind + j * dimgp;
883
884 DeltaHybrid = vct[nlev + j + 2] - vct[nlev + j + 1];
885 #ifdef HAVE_OPENMP4
886 #pragma omp parallel for simd
887 #endif
888 for (int i = 0; i < dimgp; i++)
889 {
890 omega[i + dimgp]
891 = omega[i] - diver[i] * (halfp[i + dimgp] - halfp[i]) - DeltaHybrid * (uwind[i] * dpsdx[i] + vwind[i] * dpsdy[i]);
892 }
893 }
894
895 /* interpolate to full levels */
896
897 for (int j = 0; j < nlev; j++)
898 {
899 omega = omega_in + j * dimgp;
900 #ifdef HAVE_OPENMP4
901 #pragma omp parallel for simd
902 #endif
903 for (int i = 0; i < dimgp; i++) omega[i] = 0.5 * (omega[i] + omega[i + dimgp]);
904 }
905
906 /* compute full level part of vertical velocity */
907
908 #ifdef _OPENMP
909 #pragma omp parallel for default(shared) private(omega, halfp, fullp, uwind, vwind, DeltaHybrid, Cterm, Pterm)
910 #endif
911 for (int j = 0; j < nlev; j++)
912 {
913 omega = omega_in + j * dimgp;
914 halfp = halfpress + j * dimgp;
915 fullp = fullpress + j * dimgp;
916 uwind = u_wind + j * dimgp;
917 vwind = v_wind + j * dimgp;
918
919 DeltaHybrid = vct[nlev + j + 2] - vct[nlev + j + 1];
920 if (std::fabs(DeltaHybrid) > 0)
921 {
922 Cterm = vct[j + 1] * vct[nlev + j + 1] - vct[j] * vct[nlev + j + 2];
923 #ifdef HAVE_OPENMP4
924 #pragma omp simd
925 #endif
926 for (int i = 0; i < dimgp; i++)
927 {
928 if (IS_NOT_EQUAL(Cterm, 0.0))
929 Pterm = Cterm / (halfp[i + dimgp] - halfp[i]) * std::log(halfp[i + dimgp] / halfp[i]);
930 else
931 Pterm = 0.0;
932
933 omega[i]
934 += fullp[i] * (uwind[i] * dpsdx[i] + vwind[i] * dpsdy[i]) / (halfp[i + dimgp] - halfp[i]) * (DeltaHybrid + Pterm);
935 }
936 }
937 }
938 }
939
940 void
MakeGeopotHeight(double * geop,double * gt,double * gq,double * ph,int nhor,int nlev)941 MakeGeopotHeight(double *geop, double *gt, double *gq, double *ph, int nhor, int nlev)
942 {
943 const double z2log2 = 2.0 * std::log(2.0);
944 const double vtmp = (C_RV / PlanetRD) - 1.0;
945 const double zrg = 1.0 / PlanetGrav;
946
947 if (gq) /* Humidity is present */
948 {
949 for (int j = nlev; j > 1; j--)
950 {
951 double *geopl = geop + nhor * (j - 1);
952 const double *restrict gtl = gt + nhor * (j - 1);
953 const double *restrict gql = gq + nhor * (j - 1);
954 const double *restrict phl = ph + nhor * (j - 1);
955 #ifdef HAVE_OPENMP4
956 #pragma omp parallel for simd
957 #endif
958 for (int i = 0; i < nhor; i++)
959 geopl[i] = geopl[i + nhor] + PlanetRD * gtl[i] * (1.0 + vtmp * gql[i]) * std::log(phl[i + nhor] / phl[i]);
960 }
961
962 #ifdef HAVE_OPENMP4
963 #pragma omp parallel for simd
964 #endif
965 for (int i = 0; i < nhor; i++) geop[i] = geop[i + nhor] + PlanetRD * gt[i] * (1.0 + vtmp * gq[i]) * z2log2;
966 }
967 else // No humidity
968 {
969 double *geopl = geop + nhor;
970 const double *phl = ph + nhor;
971
972 for (int j = nlev; j > 1; j--)
973 #ifdef HAVE_OPENMP4
974 #pragma omp simd
975 #endif
976 for (int i = nhor * (j - 1); i < nhor * j; i++) geop[i] = geopl[i] + PlanetRD * gt[i] * std::log(phl[i] / ph[i]);
977
978 #ifdef HAVE_OPENMP4
979 #pragma omp simd
980 #endif
981 for (int i = 0; i < nhor; i++) geop[i] = geopl[i] + PlanetRD * gt[i] * z2log2;
982 }
983
984 #ifdef HAVE_OPENMP4
985 #pragma omp parallel for simd
986 #endif
987 for (int i = 0; i < nhor * (nlev + 1); i++) geop[i] *= zrg;
988 }
989
990 constexpr double SCALESLP = 101325.0;
991
992 /* ======================================== */
993 /* LayerWater integral liquid water content */
994 /* ======================================== */
995
996 void
LayerWater(double * ww,double * ll,double pmax,double pmin,int DimGP,int HalfLevels,double * vct)997 LayerWater(double *ww, double *ll, double pmax, double pmin, int DimGP, int HalfLevels, double *vct)
998 {
999 double pph[MaxLevel];
1000
1001 int k;
1002 for (k = 0; k < HalfLevels; k++) pph[k] = vct[k] + vct[k + HalfLevels] * SCALESLP;
1003 for (k = 0; k < HalfLevels; k++)
1004 if (pph[k] > pmax) break;
1005 const auto MaxLev = k - 1;
1006 for (k = HalfLevels - 1; k >= 0; k--)
1007 if (pph[k] < pmin) break;
1008 const auto MinLev = k;
1009
1010 varray_fill(DimGP, ll, 0.0);
1011
1012 for (k = MaxLev; k <= MinLev; k++)
1013 {
1014 for (int i = 0; i < DimGP; i++) ll[i] += ww[i + k * DimGP] * (pph[k + 1] - pph[k]);
1015 }
1016 for (int i = 0; i < DimGP; i++) ll[i] /= PlanetGrav;
1017 }
1018
1019 /* ================================================= */
1020 /* LayerCloud calculates random overlap cloud cover */
1021 /* ================================================= */
1022
1023 void
LayerCloud(double * cc,double * ll,double pmax,double pmin,int DimGP,int HalfLevels,double * vct)1024 LayerCloud(double *cc, double *ll, double pmax, double pmin, int DimGP, int HalfLevels, double *vct)
1025 {
1026 double pph[MaxLevel];
1027 constexpr double ZEPSEC = 1.0e-12;
1028
1029 int k;
1030 for (k = 0; k < HalfLevels; k++) pph[k] = vct[k] + vct[k + HalfLevels] * SCALESLP;
1031 for (k = 0; k < HalfLevels; k++)
1032 if (pph[k] > pmax) break;
1033 const auto MaxLev = k - 1;
1034 for (k = HalfLevels - 1; k >= 0; k--)
1035 if (pph[k] < pmin) break;
1036 const auto MinLev = k;
1037
1038 for (int i = 0; i < DimGP; i++) ll[i] = 1. - cc[i + MaxLev * DimGP];
1039
1040 for (k = MaxLev + 1; k <= MinLev; k++)
1041 {
1042 for (int i = 0; i < DimGP; i++)
1043 ll[i]
1044 *= (1. - std::max(cc[i + (k - 1) * DimGP], cc[i + k * DimGP])) / (1. - std::min(cc[i + (k - 1) * DimGP], 1. - ZEPSEC));
1045 }
1046 for (int i = 0; i < DimGP; i++) ll[i] = 1. - ll[i];
1047 }
1048
1049 // Grid Point Computations
1050 void
after_EchamCompGP(const AfterControl & globs,struct Variable * vars)1051 after_EchamCompGP(const AfterControl &globs, struct Variable *vars)
1052 {
1053 if (vars[GEOPOTHEIGHT].comp || vars[SLP].comp || vars[THETAF].needed || vars[HALF_PRESS].needed || vars[RHUMIDITY].comp
1054 || vars[OMEGA].comp || globs.Type >= 30)
1055 {
1056 if (vars[FULL_PRESS].hybrid == nullptr) vars[FULL_PRESS].hybrid = alloc_dp(globs.Dim3GP, "vars[FULL_PRESS].hybrid");
1057
1058 vars[HALF_PRESS].hlev = globs.NumLevel + 1;
1059 vars[HALF_PRESS].plev = globs.NumLevelRequest;
1060 vars[HALF_PRESS].sfit = false;
1061
1062 if (vars[HALF_PRESS].hybrid == nullptr)
1063 vars[HALF_PRESS].hybrid = alloc_dp(globs.Dim3GP + globs.DimGP, "vars[HALF_PRESS].hybrid");
1064
1065 vct_to_hybrid_pressure(vars[FULL_PRESS].hybrid, vars[HALF_PRESS].hybrid, globs.vct, vars[PS_PROG].hybrid, globs.NumLevel,
1066 globs.DimGP);
1067 }
1068
1069 if (globs.unitsel > 2) vars[FULL_PRESS].hybrid = (double *) FreeMemory(vars[FULL_PRESS].hybrid);
1070
1071 if (vars[THETAF].needed)
1072 {
1073 vars[THETAF].hlev = globs.NumLevel;
1074 vars[THETAF].plev = globs.NumLevelRequest;
1075 vars[THETAF].sfit = true;
1076 if (vars[THETAF].hybrid == nullptr) vars[THETAF].hybrid = alloc_dp(globs.Dim3GP, "vars[THETAF].hybrid");
1077 if (vars[THETAH].hybrid == nullptr) vars[THETAH].hybrid = alloc_dp(globs.Dim3GP, "vars[THETAH].hybrid");
1078 theta(vars[THETAF].hybrid, vars[THETAH].hybrid, vars[HALF_PRESS].hybrid, vars[PS_PROG].hybrid, vars[TEMPERATURE].hybrid,
1079 vars[TS].hybrid, globs.NumLevel, globs.DimGP, globs.Dim3GP);
1080 }
1081
1082 if (vars[GEOPOTHEIGHT].comp)
1083 {
1084 vars[GEOPOTHEIGHT].hlev = globs.NumLevel + 1;
1085 vars[GEOPOTHEIGHT].plev = globs.NumLevelRequest;
1086 vars[GEOPOTHEIGHT].sfit = true;
1087 vars[GEOPOTHEIGHT].hybrid = alloc_dp(globs.Dim3GP + globs.DimGP, "vars[GEOPOTHEIGHT].hybrid");
1088
1089 array_copy(globs.DimGP, globs.Orography, vars[GEOPOTHEIGHT].hybrid + globs.Dim3GP);
1090 MakeGeopotHeight(vars[GEOPOTHEIGHT].hybrid, vars[TEMPERATURE].hybrid, vars[HUMIDITY].hybrid, vars[HALF_PRESS].hybrid,
1091 globs.DimGP, globs.NumLevel);
1092
1093 vars[HUMIDITY].needed = vars[HUMIDITY].selected;
1094 }
1095 else if (vars[GEOPOTHEIGHT].hybrid && vars[GEOPOTHEIGHT].hlev == globs.NumLevel)
1096 {
1097 vars[GEOPOTHEIGHT].hlev = globs.NumLevel + 1;
1098 vars[GEOPOTHEIGHT].sfit = true;
1099 vars[GEOPOTHEIGHT].hybrid = (double *) Realloc(vars[GEOPOTHEIGHT].hybrid, (globs.Dim3GP + globs.DimGP) * sizeof(double));
1100 array_copy(globs.DimGP, globs.Orography, vars[GEOPOTHEIGHT].hybrid + globs.Dim3GP);
1101 for (int i = 0; i < globs.DimGP; i++) vars[GEOPOTHEIGHT].hybrid[globs.Dim3GP + i] /= PlanetGrav;
1102 }
1103
1104 if (vars[DPSDX].needed || vars[DPSDY].needed)
1105 for (int l = 0; l < globs.DimGP; l++)
1106 {
1107 vars[DPSDX].hybrid[l] *= vars[PS_PROG].hybrid[l];
1108 vars[DPSDY].hybrid[l] *= vars[PS_PROG].hybrid[l];
1109 }
1110
1111 if (vars[OMEGA].comp)
1112 {
1113 vars[OMEGA].hlev = globs.NumLevel + 1;
1114 vars[OMEGA].plev = globs.NumLevelRequest;
1115 vars[OMEGA].sfit = true;
1116 vars[OMEGA].hybrid = alloc_dp(globs.Dim3GP + globs.DimGP, "OMEGA.hybrid");
1117
1118 Omega(vars[OMEGA].hybrid, vars[DIVERGENCE].hybrid, vars[U_WIND].hybrid, vars[V_WIND].hybrid, vars[HALF_PRESS].hybrid,
1119 vars[FULL_PRESS].hybrid, vars[DPSDX].hybrid, vars[DPSDY].hybrid, globs.vct, globs.DimGP, globs.NumLevel);
1120
1121 vars[DPSDX].needed = vars[DPSDX].selected;
1122 vars[DPSDY].needed = vars[DPSDY].selected;
1123 }
1124
1125 if (vars[WINDSPEED].comp)
1126 {
1127 vars[WINDSPEED].hlev = globs.NumLevel;
1128 vars[WINDSPEED].plev = globs.NumLevelRequest;
1129 vars[WINDSPEED].sfit = true;
1130 vars[WINDSPEED].hybrid = alloc_dp(globs.Dim3GP, "vars[WINDSPEED].hybrid");
1131
1132 windSpeed(vars[WINDSPEED].hybrid, vars[U_WIND].hybrid, vars[V_WIND].hybrid, globs.Dim3GP);
1133 }
1134
1135 if (vars[RHUMIDITY].comp)
1136 {
1137 vars[RHUMIDITY].hlev = globs.NumLevel;
1138 vars[RHUMIDITY].plev = globs.NumLevelRequest;
1139 vars[RHUMIDITY].sfit = false;
1140 vars[RHUMIDITY].hybrid = alloc_dp(globs.Dim3GP, "vars[RHUMIDITY].hybrid");
1141
1142 sh2rh(globs.AnalysisData, vars[HUMIDITY].hybrid, vars[RHUMIDITY].hybrid, vars[TEMPERATURE].hybrid, globs.NumLevel,
1143 globs.DimGP, globs.LevelRequest, vars[FULL_PRESS].hybrid);
1144
1145 vars[TEMPERATURE].needed = vars[TEMPERATURE].selected;
1146 vars[HUMIDITY].needed = vars[HUMIDITY].selected;
1147 }
1148
1149 if (vars[PS].comp)
1150 {
1151 vars[PS].hlev = 1;
1152 vars[PS].plev = 1;
1153 vars[PS].sfit = true; /* ??? */
1154 vars[PS].hybrid = alloc_dp(globs.DimGP, "vars[PS].hybrid");
1155 for (int l = 0; l < globs.DimGP; l++) vars[PS].hybrid[l] = std::exp(vars[LNPS].hybrid[l]);
1156 }
1157
1158 if (vars[SLP].comp)
1159 {
1160 vars[SLP].hlev = 1;
1161 vars[SLP].plev = 1;
1162 vars[SLP].sfit = true;
1163 vars[SLP].hybrid = alloc_dp(globs.DimGP, "vars[SLP].hybrid");
1164
1165 extrapolate_P(vars[SLP].hybrid, vars[HALF_PRESS].hybrid + globs.Dim3GP, vars[FULL_PRESS].hybrid + globs.Dim3GP - globs.DimGP,
1166 globs.Orography, vars[TEMPERATURE].hybrid + globs.Dim3GP - globs.DimGP, globs.DimGP);
1167 vars[TEMPERATURE].needed = vars[TEMPERATURE].selected || vars[GEOPOTHEIGHT].selected;
1168 }
1169
1170 if (vars[PRECIP].comp)
1171 {
1172 vars[PRECIP].hlev = vars[PRECIP].plev = 1;
1173 vars[PRECIP].sfit = false;
1174 vars[PRECIP].hybrid = alloc_dp(globs.DimGP, "PRECIP.hybrid");
1175 Add2Vectors(vars[PRECIP].hybrid, vars[142].hybrid, vars[143].hybrid, globs.DimGP);
1176 }
1177
1178 if (vars[NET_TOP].comp)
1179 {
1180 vars[NET_TOP].hlev = vars[NET_TOP].plev = 1;
1181 vars[NET_TOP].sfit = false;
1182 vars[NET_TOP].hybrid = alloc_dp(globs.DimGP, "NET_TOP.hybrid");
1183 Add2Vectors(vars[NET_TOP].hybrid, vars[178].hybrid, vars[179].hybrid, globs.DimGP);
1184 }
1185
1186 if (vars[NET_BOT].comp)
1187 {
1188 vars[NET_BOT].hlev = vars[NET_BOT].plev = 1;
1189 vars[NET_BOT].sfit = false;
1190 vars[NET_BOT].hybrid = alloc_dp(globs.DimGP, "NET_BOT.hybrid");
1191 Add2Vectors(vars[NET_BOT].hybrid, vars[176].hybrid, vars[177].hybrid, globs.DimGP);
1192 }
1193
1194 if (vars[NET_HEAT].comp)
1195 {
1196 vars[NET_HEAT].hlev = vars[NET_HEAT].plev = 1;
1197 vars[NET_HEAT].sfit = false;
1198 vars[NET_HEAT].hybrid = alloc_dp(globs.DimGP, "NET_HEAT.hybrid");
1199 /*
1200 if ( Source == S_ECHAM5 )
1201 {
1202 MultVectorScalar(vars[NET_HEAT].hybrid, vars[218].hybrid, (-3.345e5),
1203 globs.DimGP, vars[218].nmiss, vars[218].missval);
1204 Add2Vectors(vars[NET_HEAT].hybrid, vars[NET_HEAT].hybrid,
1205 vars[176].hybrid, globs.DimGP); Add2Vectors(vars[NET_HEAT].hybrid,
1206 vars[NET_HEAT].hybrid, vars[177].hybrid, globs.DimGP);
1207 Add2Vectors(vars[NET_HEAT].hybrid, vars[NET_HEAT].hybrid,
1208 vars[146].hybrid, globs.DimGP); Add2Vectors(vars[NET_HEAT].hybrid,
1209 vars[NET_HEAT].hybrid, vars[147].hybrid, globs.DimGP);
1210 Add2Vectors(vars[NET_HEAT].hybrid, vars[NET_HEAT].hybrid,
1211 vars[206].hybrid, globs.DimGP); Sub2Vectors(vars[NET_HEAT].hybrid,
1212 vars[NET_HEAT].hybrid, vars[208].hybrid, globs.DimGP);
1213 Sub2Vectors(vars[NET_HEAT].hybrid, vars[NET_HEAT].hybrid,
1214 vars[209].hybrid, globs.DimGP);
1215 }
1216 else
1217 */
1218 {
1219 MultVectorScalar(vars[NET_HEAT].hybrid, vars[218].hybrid, C_TIMES_RHOH2O, globs.DimGP, vars[218].nmiss, vars[218].missval);
1220 Add2Vectors(vars[NET_HEAT].hybrid, vars[NET_HEAT].hybrid, vars[176].hybrid, globs.DimGP);
1221 Add2Vectors(vars[NET_HEAT].hybrid, vars[NET_HEAT].hybrid, vars[177].hybrid, globs.DimGP);
1222 Add2Vectors(vars[NET_HEAT].hybrid, vars[NET_HEAT].hybrid, vars[146].hybrid, globs.DimGP);
1223 Add2Vectors(vars[NET_HEAT].hybrid, vars[NET_HEAT].hybrid, vars[147].hybrid, globs.DimGP);
1224 Sub2Vectors(vars[NET_HEAT].hybrid, vars[NET_HEAT].hybrid, vars[220].hybrid, globs.DimGP);
1225 }
1226 }
1227
1228 if (vars[NET_WATER].comp)
1229 {
1230 vars[NET_WATER].hlev = vars[NET_WATER].plev = 1;
1231 vars[NET_WATER].sfit = false;
1232 vars[NET_WATER].hybrid = alloc_dp(globs.DimGP, "NET_WATER.hybrid");
1233 Sub2Vectors(vars[NET_WATER].hybrid, vars[182].hybrid, vars[160].hybrid, globs.DimGP);
1234 Add2Vectors(vars[NET_WATER].hybrid, vars[NET_WATER].hybrid, vars[142].hybrid, globs.DimGP);
1235 Add2Vectors(vars[NET_WATER].hybrid, vars[NET_WATER].hybrid, vars[143].hybrid, globs.DimGP);
1236 }
1237
1238 if (vars[LOW_WATER].comp)
1239 {
1240 vars[LOW_WATER].hlev = vars[LOW_WATER].plev = 1;
1241 vars[LOW_WATER].sfit = false;
1242 vars[LOW_WATER].hybrid = alloc_dp(globs.DimGP, "vars[LOW_WATER].hybrid");
1243 LayerWater(vars[222].hybrid, vars[LOW_WATER].hybrid, 75000., 101300., globs.DimGP, globs.HalfLevels, globs.vct);
1244 }
1245
1246 if (vars[MID_WATER].comp)
1247 {
1248 vars[MID_WATER].hlev = vars[MID_WATER].plev = 1;
1249 vars[MID_WATER].sfit = false;
1250 vars[MID_WATER].hybrid = alloc_dp(globs.DimGP, "vars[MID_WATER].hybrid");
1251 LayerWater(vars[222].hybrid, vars[MID_WATER].hybrid, 46000., 73000., globs.DimGP, globs.HalfLevels, globs.vct);
1252 }
1253
1254 if (vars[HIH_WATER].comp)
1255 {
1256 vars[HIH_WATER].hlev = vars[HIH_WATER].plev = 1;
1257 vars[HIH_WATER].sfit = false;
1258 vars[HIH_WATER].hybrid = alloc_dp(globs.DimGP, "vars[HIH_WATER].hybrid");
1259 LayerWater(vars[222].hybrid, vars[HIH_WATER].hybrid, 5000., 44000., globs.DimGP, globs.HalfLevels, globs.vct);
1260 }
1261
1262 if (vars[ALL_WATER].comp)
1263 {
1264 vars[ALL_WATER].hlev = vars[ALL_WATER].plev = 1;
1265 vars[ALL_WATER].sfit = false;
1266 vars[ALL_WATER].hybrid = alloc_dp(globs.DimGP, "vars[ALL_WATER].hybrid");
1267 LayerWater(vars[222].hybrid, vars[ALL_WATER].hybrid, 5000., 101300., globs.DimGP, globs.HalfLevels, globs.vct);
1268 }
1269
1270 if (vars[LOW_CLOUD].comp)
1271 {
1272 vars[LOW_CLOUD].hlev = vars[LOW_CLOUD].plev = 1;
1273 vars[LOW_CLOUD].sfit = false;
1274 vars[LOW_CLOUD].hybrid = alloc_dp(globs.DimGP, "vars[LOW_CLOUD].hybrid");
1275 LayerCloud(vars[223].hybrid, vars[LOW_CLOUD].hybrid, 75000., 101300., globs.DimGP, globs.HalfLevels, globs.vct);
1276 }
1277
1278 if (vars[MID_CLOUD].comp)
1279 {
1280 vars[MID_CLOUD].hlev = vars[MID_CLOUD].plev = 1;
1281 vars[MID_CLOUD].sfit = false;
1282 vars[MID_CLOUD].hybrid = alloc_dp(globs.DimGP, "vars[MID_CLOUD].hybrid");
1283 LayerCloud(vars[223].hybrid, vars[MID_CLOUD].hybrid, 46000., 73000., globs.DimGP, globs.HalfLevels, globs.vct);
1284 }
1285
1286 if (vars[HIH_CLOUD].comp)
1287 {
1288 vars[HIH_CLOUD].hlev = vars[HIH_CLOUD].plev = 1;
1289 vars[HIH_CLOUD].sfit = false;
1290 vars[HIH_CLOUD].hybrid = alloc_dp(globs.DimGP, "vars[HIH_CLOUD].hybrid");
1291 LayerCloud(vars[223].hybrid, vars[HIH_CLOUD].hybrid, 5000., 44000., globs.DimGP, globs.HalfLevels, globs.vct);
1292 }
1293
1294 if (vars[SW_CLF].comp)
1295 {
1296 vars[SW_CLF].hlev = vars[SW_CLF].plev = 1;
1297 vars[SW_CLF].sfit = false;
1298 vars[SW_CLF].hybrid = alloc_dp(globs.DimGP, "SW_CLF.hybrid");
1299 Sub2Vectors(vars[SW_CLF].hybrid, vars[178].hybrid, vars[224].hybrid, globs.DimGP);
1300 }
1301
1302 if (vars[SW_BOT_CLF].comp)
1303 {
1304 vars[SW_BOT_CLF].hlev = vars[SW_BOT_CLF].plev = 1;
1305 vars[SW_BOT_CLF].sfit = false;
1306 vars[SW_BOT_CLF].hybrid = alloc_dp(globs.DimGP, "vars[SW_BOT_CLF].hybrid");
1307 Sub2Vectors(vars[SW_BOT_CLF].hybrid, vars[176].hybrid, vars[185].hybrid, globs.DimGP);
1308 }
1309
1310 if (vars[SW_TOP_CLF].comp)
1311 {
1312 vars[SW_TOP_CLF].hlev = vars[SW_TOP_CLF].plev = 1;
1313 vars[SW_TOP_CLF].sfit = false;
1314 vars[SW_TOP_CLF].hybrid = alloc_dp(globs.DimGP, "vars[SW_TOP_CLF].hybrid");
1315 Sub2Vectors(vars[SW_TOP_CLF].hybrid, vars[178].hybrid, vars[187].hybrid, globs.DimGP);
1316 }
1317
1318 if (vars[LW_CLF].comp)
1319 {
1320 vars[LW_CLF].hlev = vars[LW_CLF].plev = 1;
1321 vars[LW_CLF].sfit = false;
1322 vars[LW_CLF].hybrid = alloc_dp(globs.DimGP, "LW_CLF.hybrid");
1323 Sub2Vectors(vars[LW_CLF].hybrid, vars[179].hybrid, vars[225].hybrid, globs.DimGP);
1324 }
1325
1326 if (vars[LW_BOT_CLF].comp)
1327 {
1328 vars[LW_BOT_CLF].hlev = vars[LW_BOT_CLF].plev = 1;
1329 vars[LW_BOT_CLF].sfit = false;
1330 vars[LW_BOT_CLF].hybrid = alloc_dp(globs.DimGP, "vars[LW_BOT_CLF].hybrid");
1331 Sub2Vectors(vars[LW_BOT_CLF].hybrid, vars[177].hybrid, vars[186].hybrid, globs.DimGP);
1332 }
1333
1334 if (vars[LW_TOP_CLF].comp)
1335 {
1336 vars[LW_TOP_CLF].hlev = vars[LW_TOP_CLF].plev = 1;
1337 vars[LW_TOP_CLF].sfit = false;
1338 vars[LW_TOP_CLF].hybrid = alloc_dp(globs.DimGP, "vars[LW_TOP_CLF].hybrid");
1339 Sub2Vectors(vars[LW_TOP_CLF].hybrid, vars[179].hybrid, vars[188].hybrid, globs.DimGP);
1340 }
1341
1342 if (vars[NET_CLF].comp)
1343 {
1344 vars[NET_CLF].hlev = vars[NET_CLF].plev = 1;
1345 vars[NET_CLF].sfit = false;
1346 vars[NET_CLF].hybrid = alloc_dp(globs.DimGP, "NET_CLF.hybrid");
1347 Add2Vectors(vars[NET_CLF].hybrid, vars[178].hybrid, vars[179].hybrid, globs.DimGP);
1348 Sub2Vectors(vars[NET_CLF].hybrid, vars[NET_CLF].hybrid, vars[224].hybrid, globs.DimGP);
1349 Sub2Vectors(vars[NET_CLF].hybrid, vars[NET_CLF].hybrid, vars[225].hybrid, globs.DimGP);
1350 }
1351
1352 if (vars[SW_ATM].comp)
1353 {
1354 vars[SW_ATM].hlev = vars[SW_ATM].plev = 1;
1355 vars[SW_ATM].sfit = false;
1356 vars[SW_ATM].hybrid = alloc_dp(globs.DimGP, "vars[SW_ATM].hybrid");
1357 Sub2Vectors(vars[SW_ATM].hybrid, vars[178].hybrid, vars[176].hybrid, globs.DimGP);
1358 }
1359
1360 if (vars[LW_ATM].comp)
1361 {
1362 vars[LW_ATM].hlev = vars[LW_ATM].plev = 1;
1363 vars[LW_ATM].sfit = false;
1364 vars[LW_ATM].hybrid = alloc_dp(globs.DimGP, "vars[LW_ATM].hybrid");
1365 Sub2Vectors(vars[LW_ATM].hybrid, vars[179].hybrid, vars[177].hybrid, globs.DimGP);
1366 }
1367
1368 if (vars[NET_ATM].comp)
1369 {
1370 vars[NET_ATM].hlev = vars[NET_ATM].plev = 1;
1371 vars[NET_ATM].sfit = false;
1372 vars[NET_ATM].hybrid = alloc_dp(globs.DimGP, "vars[NET_ATM].hybrid");
1373 Add2Vectors(vars[NET_ATM].hybrid, vars[178].hybrid, vars[179].hybrid, globs.DimGP);
1374 Sub2Vectors(vars[NET_ATM].hybrid, vars[NET_ATM].hybrid, vars[176].hybrid, globs.DimGP);
1375 Sub2Vectors(vars[NET_ATM].hybrid, vars[NET_ATM].hybrid, vars[177].hybrid, globs.DimGP);
1376 }
1377
1378 if (vars[SURF_RUNOFF].comp)
1379 {
1380 vars[SURF_RUNOFF].hlev = vars[SURF_RUNOFF].plev = 1;
1381 vars[SURF_RUNOFF].sfit = false;
1382 vars[SURF_RUNOFF].hybrid = alloc_dp(globs.DimGP, "vars[SURF_RUNOFF].hybrid");
1383 Sub2Vectors(vars[SURF_RUNOFF].hybrid, vars[182].hybrid, vars[221].hybrid, globs.DimGP);
1384 Add2Vectors(vars[SURF_RUNOFF].hybrid, vars[SURF_RUNOFF].hybrid, vars[142].hybrid, globs.DimGP);
1385 Add2Vectors(vars[SURF_RUNOFF].hybrid, vars[SURF_RUNOFF].hybrid, vars[143].hybrid, globs.DimGP);
1386 }
1387
1388 if (vars[FRESH_WATER].comp)
1389 {
1390 vars[FRESH_WATER].hlev = vars[FRESH_WATER].plev = 1;
1391 vars[FRESH_WATER].sfit = false;
1392 vars[FRESH_WATER].hybrid = alloc_dp(globs.DimGP, "vars[FRESH_WATER].hybrid");
1393 Add2Vectors(vars[FRESH_WATER].hybrid, vars[142].hybrid, vars[143].hybrid, globs.DimGP);
1394 Add2Vectors(vars[FRESH_WATER].hybrid, vars[FRESH_WATER].hybrid, vars[182].hybrid, globs.DimGP);
1395 }
1396 }
1397
1398 static void
Derivate(double field[],double derilam[],int levels,int Waves,int Latitudes,double DerivationFactor[])1399 Derivate(double field[], double derilam[], int levels, int Waves, int Latitudes, double DerivationFactor[])
1400 {
1401 int l, n, lev;
1402 int i;
1403
1404 i = 0;
1405 for (lev = 0; lev < levels; lev++)
1406 for (n = 0; n < Waves; n++)
1407 {
1408 for (l = 0; l < Latitudes; l++)
1409 {
1410 derilam[i] = -n * field[i + Latitudes] * DerivationFactor[l];
1411 i++;
1412 }
1413 for (l = 0; l < Latitudes; l++)
1414 {
1415 derilam[i] = n * field[i - Latitudes] * DerivationFactor[l];
1416 i++;
1417 }
1418 }
1419 }
1420
1421 /* Process Model Level data */
1422 void
after_processML(AfterControl & globs,struct Variable * vars)1423 after_processML(AfterControl &globs, struct Variable *vars)
1424 {
1425 int code, l, i;
1426 long fieldSize = 0;
1427 int lindex, nlevel;
1428 long offset;
1429 int leveltype;
1430 size_t nmiss;
1431 double *pressureLevel = nullptr;
1432
1433 globs.MeanCount++;
1434 globs.TermCount++;
1435
1436 if (globs.MeanCount == 1)
1437 {
1438 if (globs.Debug) cdo_print("TermCount = %d MeanCount = %d", globs.TermCount, globs.MeanCount);
1439 globs.StartDate = globs.OldDate;
1440 }
1441
1442 if (globs.TermCount > 120) globs.Debug = 0;
1443
1444 /* ============================== */
1445 /* Computations in spectral space */
1446 /* ============================== */
1447
1448 if (vars[U_WIND].comp || vars[V_WIND].comp)
1449 {
1450 vars[U_WIND].hlev = vars[V_WIND].hlev = vars[DIVERGENCE].hlev;
1451 vars[U_WIND].plev = vars[V_WIND].plev = vars[DIVERGENCE].plev;
1452 vars[U_WIND].sfit = vars[V_WIND].sfit = true;
1453 vars[U_WIND].spectral = alloc_dp(globs.Dim3SP, "vars[U_WIND].spectral");
1454 vars[V_WIND].spectral = alloc_dp(globs.Dim3SP, "vars[V_WIND].spectral");
1455
1456 if (vars[DIVERGENCE].spectral == nullptr) after_gp2sp(globs, vars, DIVERGENCE);
1457 if (vars[VORTICITY].spectral == nullptr) after_gp2sp(globs, vars, VORTICITY);
1458
1459 dv2uv(vars[DIVERGENCE].spectral, vars[VORTICITY].spectral, vars[U_WIND].spectral, vars[V_WIND].spectral, globs.dv2uv_f1,
1460 globs.dv2uv_f2, globs.Truncation, globs.DimSP, vars[DIVERGENCE].hlev);
1461 }
1462
1463 if (vars[VELOPOT].comp && globs.Type < 30)
1464 {
1465 vars[VELOPOT].hlev = vars[DIVERGENCE].hlev;
1466 vars[VELOPOT].plev = vars[DIVERGENCE].plev;
1467 vars[VELOPOT].spectral = alloc_dp(globs.Dim3SP, "vars[VELOPOT].spectral");
1468
1469 if (vars[DIVERGENCE].spectral == nullptr) after_gp2sp(globs, vars, DIVERGENCE);
1470
1471 dv2ps(vars[DIVERGENCE].spectral, vars[VELOPOT].spectral, vars[DIVERGENCE].hlev, globs.Truncation);
1472 }
1473
1474 if (vars[STREAM].comp && globs.Type < 30)
1475 {
1476 vars[STREAM].hlev = vars[VORTICITY].hlev;
1477 vars[STREAM].plev = vars[VORTICITY].plev;
1478 vars[STREAM].spectral = alloc_dp(globs.Dim3SP, "vars[STREAM].spectral");
1479
1480 if (vars[VORTICITY].spectral == nullptr) after_gp2sp(globs, vars, VORTICITY);
1481
1482 dv2ps(vars[VORTICITY].spectral, vars[STREAM].spectral, vars[VORTICITY].hlev, globs.Truncation);
1483 }
1484
1485 if (vars[VORTICITY].spectral && !vars[VORTICITY].needed)
1486 vars[VORTICITY].spectral = (double *) FreeMemory(vars[VORTICITY].spectral);
1487
1488 if (vars[DIVERGENCE].spectral && !vars[DIVERGENCE].needed)
1489 vars[DIVERGENCE].spectral = (double *) FreeMemory(vars[DIVERGENCE].spectral);
1490
1491 /* --------------------------- */
1492 /* Output of spectral fields */
1493 /* --------------------------- */
1494
1495 if (globs.Type == 0)
1496 {
1497 for (code = 0; code < MaxCodes; code++)
1498 if (vars[code].selected)
1499 {
1500 if (!vars[code].spectral) afterAbort("Code %d not available on spectral space!", code);
1501
1502 nlevel = zaxisInqSize(vars[code].ozaxisID);
1503 for (lindex = 0; lindex < nlevel; lindex++)
1504 {
1505 offset = lindex * globs.DimSP;
1506 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
1507 streamWriteRecord(globs.ostreamID, vars[code].spectral + offset, 0);
1508 }
1509 }
1510
1511 FreeSpectral(vars);
1512 return;
1513 }
1514
1515 /* ------------------------------- */
1516 /* Computations in fourier space */
1517 /* ------------------------------- */
1518
1519 if (globs.Type >= 10)
1520 {
1521 for (code = 0; code < MaxCodes; code++)
1522 if (vars[code].needed && vars[code].spectral)
1523 {
1524 if (vars[code].fourier == nullptr)
1525 {
1526 fieldSize = vars[code].hlev * globs.DimFC;
1527 vars[code].fourier = alloc_dp(fieldSize, FieldName(code, "fourier"));
1528 sp2fc(vars[code].spectral, vars[code].fourier, globs.poli, vars[code].hlev, globs.Latitudes, globs.Fouriers,
1529 globs.Truncation);
1530 }
1531 if (code != LNPS) vars[code].spectral = (double *) FreeMemory(vars[code].spectral);
1532 }
1533
1534 if (vars[U_WIND].needed && vars[U_WIND].fourier)
1535 scaluv(vars[U_WIND].fourier, globs.rcoslat, globs.Latitudes, globs.Fouriers * globs.NumLevel);
1536 if (vars[V_WIND].needed && vars[V_WIND].fourier)
1537 scaluv(vars[V_WIND].fourier, globs.rcoslat, globs.Latitudes, globs.Fouriers * globs.NumLevel);
1538
1539 if (vars[DPSDX].needed)
1540 {
1541 vars[DPSDX].hlev = 1;
1542 vars[DPSDX].plev = 1;
1543 vars[DPSDX].sfit = false;
1544 vars[DPSDX].fourier = alloc_dp(globs.DimFC, "vars[DPSDX].fourier");
1545 if (vars[LNPS].fourier == nullptr) after_gp2sp(globs, vars, LNPS);
1546 Derivate(vars[LNPS].fourier, vars[DPSDX].fourier, 1, globs.Waves, globs.Latitudes, globs.DerivationFactor);
1547 }
1548 if (vars[DPSDY].needed)
1549 {
1550 vars[DPSDY].hlev = 1;
1551 vars[DPSDY].plev = 1;
1552 vars[DPSDY].sfit = false;
1553 vars[DPSDY].fourier = alloc_dp(globs.DimFC, "vars[DPSDY].fourier");
1554 if (vars[LNPS].spectral == nullptr) after_gp2sp(globs, vars, LNPS);
1555 sp2fc(vars[LNPS].spectral, vars[DPSDY].fourier, globs.pdev, vars[DPSDY].hlev, globs.Latitudes, globs.Fouriers,
1556 globs.Truncation);
1557 }
1558 }
1559
1560 FreeSpectral(vars);
1561
1562 /* -------------------------- */
1563 /* Output of fourier fields */
1564 /* -------------------------- */
1565
1566 if (globs.Type == 10)
1567 {
1568 for (code = 0; code < MaxCodes; code++)
1569 if (vars[code].selected)
1570 {
1571 if (!vars[code].fourier) afterAbort("Code %d not available on fourier space!", code);
1572
1573 nlevel = zaxisInqSize(vars[code].ozaxisID);
1574 for (lindex = 0; lindex < nlevel; lindex++)
1575 {
1576 offset = lindex * globs.DimFC;
1577 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
1578 streamWriteRecord(globs.ostreamID, vars[code].fourier + offset, 0);
1579 }
1580 }
1581
1582 FreeFourier(vars);
1583 return;
1584 }
1585
1586 /* ----------------------- */
1587 /* Output of zonal means */
1588 /* ----------------------- */
1589
1590 if (globs.Type == 11)
1591 {
1592 for (code = 0; code < MaxCodes; code++)
1593 if (vars[code].selected)
1594 {
1595 if (!vars[code].fourier) afterAbort("Code %d not available on zonal mean!", code);
1596
1597 nlevel = zaxisInqSize(vars[code].ozaxisID);
1598 for (lindex = 0; lindex < nlevel; lindex++)
1599 {
1600 offset = lindex * globs.DimFC;
1601 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
1602 streamWriteRecord(globs.ostreamID, vars[code].fourier + offset, 0);
1603 }
1604 }
1605
1606 FreeFourier(vars);
1607 return;
1608 }
1609
1610 /* ------------------------------ */
1611 /* Transformation to gridpoints */
1612 /* ------------------------------ */
1613
1614 if (globs.Type >= 20)
1615 {
1616 for (code = 0; code < MaxCodes; code++)
1617 if (vars[code].needed && vars[code].fourier)
1618 {
1619 if (vars[code].hybrid == nullptr)
1620 {
1621 fieldSize = globs.DimGP * vars[code].hlev;
1622 vars[code].hybrid = alloc_dp(fieldSize, FieldName(code, "hybrid"));
1623 after_FC2GP(vars[code].fourier, vars[code].hybrid, globs.Latitudes, globs.Longitudes, vars[code].hlev,
1624 globs.Fouriers);
1625 }
1626 vars[code].fourier = (double *) FreeMemory(vars[code].fourier);
1627 }
1628
1629 if (vars[PS_PROG].comp && vars[PS_PROG].hybrid == nullptr)
1630 {
1631 vars[PS_PROG].hybrid = alloc_dp(globs.DimGP, "PS_PROG");
1632 if (vars[LNPS].hybrid)
1633 {
1634 for (l = 0; l < globs.DimGP; l++) vars[PS_PROG].hybrid[l] = std::exp(vars[LNPS].hybrid[l]);
1635 }
1636 else if (vars[PS].hybrid)
1637 {
1638 cdo_warning("log surface pressure (code 152) not found - using surface pressure (code 134)!");
1639 array_copy(globs.DimGP, vars[PS].hybrid, vars[PS_PROG].hybrid);
1640 }
1641 else
1642 {
1643 afterAbort("surface pressure not found!");
1644 }
1645 }
1646 vars[LNPS].needed = vars[LNPS].selected;
1647
1648 if (globs.Orography == nullptr)
1649 {
1650 globs.Orography = alloc_dp(globs.DimGP, "Orography");
1651 if (vars[GEOPOTENTIAL].hybrid)
1652 array_copy(globs.DimGP, vars[GEOPOTENTIAL].hybrid, globs.Orography);
1653 else
1654 {
1655 if (vars[GEOPOTENTIAL].selected || globs.Type >= 30)
1656 {
1657 cdo_warning("Orography not found - using zero orography!");
1658 varray_fill(globs.DimGP, globs.Orography, 0.0);
1659 }
1660 }
1661 }
1662 vars[GEOPOTENTIAL].needed = vars[GEOPOTENTIAL].selected;
1663
1664 after_EchamCompGP(globs, vars);
1665 }
1666
1667 FreeFourier(vars);
1668
1669 if (globs.Type == 20)
1670 {
1671 /* ----------------------- */
1672 /* Means on hybrid grids */
1673 /* ----------------------- */
1674
1675 if (globs.Mean)
1676 {
1677 for (code = 0; code < MaxCodes; code++)
1678 {
1679 if (vars[code].selected && vars[code].hybrid)
1680 {
1681 fieldSize = globs.DimGP * vars[code].hlev;
1682
1683 if (vars[code].mean == nullptr) vars[code].mean = alloc_dp(fieldSize, FieldName(code, "mean"));
1684
1685 if (globs.Mean > 1 && vars[code].variance == nullptr)
1686 vars[code].variance = alloc_dp(fieldSize, FieldName(code, "variance"));
1687
1688 if (globs.MeanCount == 1)
1689 {
1690 array_copy(fieldSize, vars[code].hybrid, vars[code].mean);
1691 if (globs.Mean > 1) IniQuaSum(vars[code].variance, vars[code].hybrid, fieldSize);
1692 }
1693 else
1694 {
1695 AddVector(vars[code].mean, vars[code].hybrid, fieldSize, &vars[code].nmiss, vars[code].missval);
1696 if (globs.Mean > 1) AddQuaSum(vars[code].variance, vars[code].hybrid, fieldSize);
1697 }
1698
1699 if (globs.EndOfInterval)
1700 {
1701 if (vars[code].samp == nullptr)
1702 MultVectorScalar(vars[code].hybrid, vars[code].mean, 1.0 / globs.MeanCount, fieldSize, vars[code].nmiss,
1703 vars[code].missval);
1704 else
1705 DivVectorIvector(vars[code].hybrid, vars[code].mean, vars[code].samp, fieldSize, &vars[code].nmiss,
1706 vars[code].missval);
1707
1708 if (globs.Mean > 1) VarQuaSum(vars[code].variance, vars[code].hybrid, fieldSize, globs.MeanCount);
1709 }
1710 }
1711 }
1712 }
1713
1714 /* ---------------------------- */
1715 /* Output of hybrid level grids */
1716 /* ---------------------------- */
1717
1718 if (globs.Mean == 0 || globs.EndOfInterval)
1719 {
1720 for (code = 0; code < MaxCodes; code++)
1721 if (vars[code].selected)
1722 {
1723 if (vars[code].hybrid == nullptr) afterAbort("Internal problem. Code %d not allocated!", code);
1724
1725 nlevel = zaxisInqSize(vars[code].ozaxisID);
1726 for (lindex = 0; lindex < nlevel; lindex++)
1727 {
1728 offset = lindex * globs.DimGP;
1729 if (globs.Mean != 2)
1730 {
1731 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
1732 streamWriteRecord(globs.ostreamID, vars[code].hybrid + offset, vars[code].nmiss);
1733 }
1734 if (globs.Mean >= 2)
1735 {
1736 streamDefRecord(globs.ostreamID2, vars[code].ovarID2, lindex);
1737 streamWriteRecord(globs.ostreamID2, vars[code].variance + offset, vars[code].nmiss);
1738 }
1739 }
1740 }
1741
1742 FreeSamp(vars);
1743 }
1744
1745 FreeHybrid(vars);
1746 return;
1747 }
1748
1749 /* -------------------------------------- */
1750 /* Vertical interpolation / extrapolation */
1751 /* -------------------------------------- */
1752
1753 if (globs.Type >= 30)
1754 {
1755 if (globs.vertIndex == nullptr) globs.vertIndex = (int *) Malloc(globs.NumLevelRequest * globs.DimGP * sizeof(int));
1756
1757 if (globs.unitsel)
1758 {
1759 if (globs.p_of_height == nullptr) globs.p_of_height = alloc_dp(globs.NumLevelRequest, "p_of_height");
1760 height_to_pressure(globs.p_of_height, globs.LevelRequest, globs.NumLevelRequest);
1761 pressureLevel = globs.p_of_height;
1762 }
1763 else
1764 {
1765 pressureLevel = globs.LevelRequest;
1766 }
1767
1768 gen_vert_index(globs.vertIndex, pressureLevel, vars[FULL_PRESS].hybrid, globs.DimGP, globs.NumLevelRequest, globs.NumLevel);
1769
1770 nmiss = 0;
1771 if (!globs.extrapolate)
1772 {
1773 if (globs.pnmiss == nullptr) globs.pnmiss = (size_t *) Malloc(globs.NumLevelRequest * sizeof(size_t));
1774 gen_vert_index_mv(globs.vertIndex, pressureLevel, globs.DimGP, globs.NumLevelRequest, vars[PS_PROG].hybrid, globs.pnmiss);
1775 for (i = 0; i < globs.NumLevelRequest; i++) nmiss += globs.pnmiss[i];
1776 }
1777
1778 for (code = 0; code < MaxCodes; code++)
1779 if (vars[code].needed && vars[code].hybrid)
1780 {
1781 leveltype = zaxisInqType(vars[code].izaxisID);
1782 if (vars[code].hlev == 1 || leveltype != ZAXIS_HYBRID || (vars[code].hlev < globs.NumLevel))
1783 {
1784 if (vars[code].grid) FreeMemory(vars[code].grid);
1785 vars[code].grid = vars[code].hybrid;
1786 vars[code].hybrid = nullptr;
1787 }
1788 else
1789 {
1790 if (vars[code].grid == nullptr)
1791 {
1792 fieldSize = globs.DimGP * globs.NumLevelRequest;
1793 vars[code].grid = alloc_dp(fieldSize, FieldName(code, "grid"));
1794 }
1795
1796 if (code == TEMPERATURE)
1797 {
1798 vertical_interp_T(globs.Orography, vars[TEMPERATURE].hybrid, vars[TEMPERATURE].grid, vars[FULL_PRESS].hybrid,
1799 vars[HALF_PRESS].hybrid, globs.vertIndex, pressureLevel, globs.NumLevelRequest, globs.DimGP,
1800 globs.NumLevel, vars[code].missval);
1801 }
1802 else if (code == GEOPOTHEIGHT)
1803 {
1804 if (vars[TEMPERATURE].hybrid == nullptr) afterAbort("Code 130 not found!");
1805 vertical_interp_Z(globs.Orography, vars[GEOPOTHEIGHT].hybrid, vars[GEOPOTHEIGHT].grid, vars[FULL_PRESS].hybrid,
1806 vars[HALF_PRESS].hybrid, globs.vertIndex, vars[TEMPERATURE].hybrid, pressureLevel,
1807 globs.NumLevelRequest, globs.DimGP, globs.NumLevel, vars[code].missval);
1808 }
1809 else
1810 {
1811 int numlevel = vars[code].hlev;
1812 if (code == OMEGA) numlevel = globs.NumLevel;
1813 double *hyb_press = vars[FULL_PRESS].hybrid;
1814 if (numlevel == (globs.NumLevel + 1)) hyb_press = vars[HALF_PRESS].hybrid;
1815
1816 vertical_interp_X(vars[code].hybrid, vars[code].grid, hyb_press, globs.vertIndex, pressureLevel,
1817 globs.NumLevelRequest, globs.DimGP, numlevel, vars[code].missval);
1818 }
1819
1820 if (!globs.extrapolate) vars[code].nmiss = nmiss;
1821
1822 if (code != TEMPERATURE) vars[code].hybrid = (double *) FreeMemory(vars[code].hybrid);
1823 }
1824 }
1825 }
1826
1827 vars[TEMPERATURE].needed = vars[TEMPERATURE].selected;
1828 FreeHybrid(vars);
1829 if (vars[HALF_PRESS].hybrid) vars[HALF_PRESS].hybrid = (double *) FreeMemory(vars[HALF_PRESS].hybrid);
1830
1831 /* -------------------------------- */
1832 /* Output of pressure level grids */
1833 /* -------------------------------- */
1834
1835 if (globs.Type == 30 && globs.Mean == 0)
1836 {
1837 for (code = 0; code < MaxCodes; code++)
1838 if (vars[code].selected && vars[code].grid)
1839 {
1840 nlevel = zaxisInqSize(vars[code].ozaxisID);
1841 for (lindex = 0; lindex < nlevel; lindex++)
1842 {
1843 offset = lindex * globs.DimGP;
1844 if (globs.Mean != 2)
1845 {
1846 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
1847 streamWriteRecord(globs.ostreamID, vars[code].grid + offset, vars[code].nmiss);
1848 }
1849 }
1850 }
1851
1852 FreeGrid(vars);
1853 return;
1854 }
1855
1856 /* ---------------------- */
1857 /* Computation of Means */
1858 /* ---------------------- */
1859
1860 if (globs.Type >= 30 && globs.Mean)
1861 for (code = 0; code < MaxCodes; code++)
1862 if (vars[code].needed && vars[code].grid)
1863 {
1864 fieldSize = globs.DimGP * vars[code].plev;
1865
1866 if (vars[code].mean == nullptr) vars[code].mean = alloc_dp(fieldSize, FieldName(code, "mean"));
1867
1868 if (globs.MeanCount == 1)
1869 array_copy(fieldSize, vars[code].grid, vars[code].mean);
1870 else
1871 AddVector(vars[code].mean, vars[code].grid, fieldSize, &vars[code].nmiss, vars[code].missval);
1872
1873 if (globs.EndOfInterval)
1874 {
1875 if (vars[code].samp == nullptr)
1876 MultVectorScalar(vars[code].mean, vars[code].mean, 1.0 / globs.MeanCount, fieldSize, vars[code].nmiss,
1877 vars[code].missval);
1878 else
1879 DivVectorIvector(vars[code].mean, vars[code].mean, vars[code].samp, fieldSize, &vars[code].nmiss,
1880 vars[code].missval);
1881 }
1882 }
1883
1884 /* -------------------------- */
1885 /* Computation of Variances */
1886 /* -------------------------- */
1887
1888 if (globs.Type >= 30 && globs.Mean > 1)
1889 for (code = 0; code < MaxCodes; code++)
1890 if (vars[code].needed && vars[code].mean)
1891 {
1892 fieldSize = globs.DimGP * vars[code].plev;
1893
1894 if (vars[code].variance == nullptr) vars[code].variance = alloc_dp(fieldSize, FieldName(code, "var"));
1895
1896 if (globs.MeanCount == 1)
1897 IniQuaSum(vars[code].variance, vars[code].grid, fieldSize);
1898 else
1899 AddQuaSum(vars[code].variance, vars[code].grid, fieldSize);
1900
1901 if (globs.EndOfInterval) VarQuaSum(vars[code].variance, vars[code].mean, fieldSize, globs.MeanCount);
1902 }
1903
1904 if (globs.Mean && !globs.EndOfInterval)
1905 {
1906 FreeGrid(vars);
1907 return;
1908 }
1909
1910 /* --------------------------------------------- */
1911 /* Output of pressure level means and variances */
1912 /* --------------------------------------------- */
1913
1914 if (globs.Type == 30 && globs.Mean)
1915 {
1916 for (code = 0; code < MaxCodes; code++)
1917 if (vars[code].selected && vars[code].mean)
1918 {
1919 nlevel = zaxisInqSize(vars[code].ozaxisID);
1920 for (lindex = 0; lindex < nlevel; lindex++)
1921 {
1922 offset = lindex * globs.DimGP;
1923 if (globs.Mean != 2)
1924 {
1925 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
1926 streamWriteRecord(globs.ostreamID, vars[code].mean + offset, vars[code].nmiss);
1927 }
1928 if (globs.Mean >= 2)
1929 {
1930 streamDefRecord(globs.ostreamID2, vars[code].ovarID2, lindex);
1931 streamWriteRecord(globs.ostreamID2, vars[code].variance + offset, vars[code].nmiss);
1932 }
1933 }
1934 }
1935
1936 FreeSamp(vars);
1937 FreeGrid(vars);
1938 return;
1939 }
1940
1941 /* ------------------ */
1942 /* Free mean fields */
1943 /* ------------------ */
1944
1945 if (globs.Type >= 40 && globs.Mean)
1946 for (code = 0; code < MaxCodes; code++)
1947 if (vars[code].mean)
1948 {
1949 if (vars[code].variance) vars[code].variance = (double *) FreeMemory(vars[code].variance);
1950 if (vars[code].grid) vars[code].grid = (double *) FreeMemory(vars[code].grid);
1951 vars[code].grid = vars[code].mean;
1952 vars[code].mean = nullptr;
1953 }
1954
1955 /* --------------------------------- */
1956 /* Transformation to fourier space */
1957 /* --------------------------------- */
1958
1959 if (globs.Type >= 40)
1960 {
1961 for (code = 0; code < MaxCodes; code++)
1962 if (vars[code].needed && vars[code].grid && (vars[code].sfit || globs.Type < 70))
1963 {
1964 if (vars[code].nmiss) afterAbort("Missing values for code %d unsupported with TYPE > 30!", code);
1965
1966 if (vars[code].fourier == nullptr)
1967 {
1968 fieldSize = globs.DimFC * vars[code].plev;
1969 vars[code].fourier = alloc_dp(fieldSize, FieldName(code, "fourier"));
1970 }
1971
1972 after_GP2FC(vars[code].grid, vars[code].fourier, globs.Latitudes, globs.Longitudes, vars[code].plev, globs.Fouriers);
1973
1974 if (vars[code].grid && (vars[code].sfit || globs.Type < 70)) vars[code].grid = (double *) FreeMemory(vars[code].grid);
1975 }
1976 }
1977
1978 for (code = 0; code < MaxCodes; code++)
1979 if (vars[code].grid && (vars[code].sfit || globs.Type < 70)) vars[code].grid = (double *) FreeMemory(vars[code].grid);
1980
1981 /* -------------------------- */
1982 /* Output of fourier fields */
1983 /* -------------------------- */
1984
1985 if (globs.Type == 40)
1986 {
1987 for (code = 0; code < MaxCodes; code++)
1988 if (vars[code].selected)
1989 {
1990 if (!vars[code].fourier) afterAbort("Code %d not available on fourier space!", code);
1991
1992 nlevel = zaxisInqSize(vars[code].ozaxisID);
1993 for (lindex = 0; lindex < nlevel; lindex++)
1994 {
1995 offset = lindex * globs.DimFC;
1996 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
1997 streamWriteRecord(globs.ostreamID, vars[code].fourier + offset, vars[code].nmiss);
1998 }
1999 }
2000
2001 FreeFourier(vars);
2002 return;
2003 }
2004
2005 /* --------------------- */
2006 /* Output of zonal means */
2007 /* --------------------- */
2008
2009 if (globs.Type == 41)
2010 {
2011 for (code = 0; code < MaxCodes; code++)
2012 if (vars[code].selected)
2013 {
2014 if (!vars[code].fourier) afterAbort("Code %d not available on zonal mean!", code);
2015
2016 nlevel = zaxisInqSize(vars[code].ozaxisID);
2017 for (lindex = 0; lindex < nlevel; lindex++)
2018 {
2019 offset = lindex * globs.DimFC;
2020 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
2021 streamWriteRecord(globs.ostreamID, vars[code].fourier + offset, vars[code].nmiss);
2022 }
2023 }
2024
2025 FreeFourier(vars);
2026 return;
2027 }
2028
2029 /* ---------------------------------- */
2030 /* Transformation to spectral space */
2031 /* ---------------------------------- */
2032
2033 if (globs.Type >= 50)
2034 {
2035 if (vars[U_WIND].needed && vars[U_WIND].fourier)
2036 scaluv(vars[U_WIND].fourier, globs.coslat, globs.Latitudes, globs.Fouriers * globs.NumLevelRequest);
2037 if (vars[V_WIND].needed && vars[V_WIND].fourier)
2038 scaluv(vars[V_WIND].fourier, globs.coslat, globs.Latitudes, globs.Fouriers * globs.NumLevelRequest);
2039
2040 for (code = 0; code < MaxCodes; code++)
2041 if (vars[code].needed && vars[code].fourier)
2042 {
2043 if (vars[code].spectral == nullptr)
2044 {
2045 fieldSize = vars[code].plev * globs.DimSP;
2046 vars[code].spectral = alloc_dp(fieldSize, FieldName(code, "spectral"));
2047 }
2048
2049 fc2sp(vars[code].fourier, vars[code].spectral, globs.pold, vars[code].plev, globs.Latitudes, globs.Fouriers,
2050 globs.Truncation);
2051 }
2052
2053 if (vars[DIVERGENCE].needed || vars[VORTICITY].needed || vars[VELOPOT].needed || vars[STREAM].needed)
2054 {
2055 if (vars[DIVERGENCE].spectral == nullptr)
2056 vars[DIVERGENCE].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[DIVERGENCE].spectral");
2057 if (vars[VORTICITY].spectral == nullptr)
2058 vars[VORTICITY].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[VORTICITY].spectral");
2059 if ((vars[U_WIND].fourier == 0 || vars[V_WIND].fourier == 0) && globs.NumLevelRequest)
2060 afterAbort("uwind or vwind missing!");
2061 uv2dv(vars[U_WIND].fourier, vars[V_WIND].fourier, vars[DIVERGENCE].spectral, vars[VORTICITY].spectral, globs.pol2,
2062 globs.pol3, globs.NumLevelRequest, globs.Latitudes, globs.Truncation);
2063 }
2064
2065 if (vars[VELOPOT].needed)
2066 {
2067 vars[VELOPOT].hlev = vars[DIVERGENCE].hlev;
2068 vars[VELOPOT].plev = vars[DIVERGENCE].plev;
2069 vars[VELOPOT].sfit = true;
2070 if (vars[VELOPOT].spectral == nullptr)
2071 vars[VELOPOT].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[VELOPOT].spectral");
2072 dv2ps(vars[DIVERGENCE].spectral, vars[VELOPOT].spectral, globs.NumLevelRequest, globs.Truncation);
2073 }
2074
2075 if (vars[STREAM].needed)
2076 {
2077 vars[STREAM].hlev = vars[VORTICITY].hlev;
2078 vars[STREAM].plev = vars[VORTICITY].plev;
2079 vars[STREAM].sfit = true;
2080 if (vars[STREAM].spectral == nullptr)
2081 vars[STREAM].spectral = alloc_dp(globs.DimSP * globs.NumLevelRequest, "vars[STREAM].spectral");
2082 dv2ps(vars[VORTICITY].spectral, vars[STREAM].spectral, globs.NumLevelRequest, globs.Truncation);
2083 }
2084 }
2085
2086 for (code = 0; code < MaxCodes; code++)
2087 if (vars[code].fourier && (vars[code].sfit || globs.Type < 61)) vars[code].fourier = (double *) FreeMemory(vars[code].fourier);
2088
2089 /* --------------------------- */
2090 /* Output of spectral fields */
2091 /* --------------------------- */
2092
2093 if (globs.Type == 50)
2094 {
2095 for (code = 0; code < MaxCodes; code++)
2096 if (vars[code].selected && vars[code].spectral)
2097 {
2098 nlevel = zaxisInqSize(vars[code].ozaxisID);
2099 for (lindex = 0; lindex < nlevel; lindex++)
2100 {
2101 offset = lindex * globs.DimSP;
2102 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
2103 streamWriteRecord(globs.ostreamID, vars[code].spectral + offset, 0);
2104 }
2105 }
2106
2107 FreeSpectral(vars);
2108 return;
2109 }
2110
2111 /* -------------------------------*/
2112 /* Computations in fourier space */
2113 /* -------------------------------*/
2114
2115 if (globs.Type >= 60)
2116 {
2117 for (code = 0; code < MaxCodes; code++)
2118 if (vars[code].needed && vars[code].spectral)
2119 {
2120 if (vars[code].fourier == nullptr)
2121 {
2122 fieldSize = vars[code].plev * globs.DimFC;
2123 vars[code].fourier = alloc_dp(fieldSize, FieldName(code, "fourier"));
2124 }
2125 sp2fc(vars[code].spectral, vars[code].fourier, globs.poli, vars[code].plev, globs.Latitudes, globs.Fouriers,
2126 globs.Truncation);
2127 }
2128 if (vars[U_WIND].needed && vars[U_WIND].fourier)
2129 scaluv(vars[U_WIND].fourier, globs.rcoslat, globs.Latitudes, globs.Fouriers * globs.NumLevelRequest);
2130 if (vars[V_WIND].needed && vars[V_WIND].fourier)
2131 scaluv(vars[V_WIND].fourier, globs.rcoslat, globs.Latitudes, globs.Fouriers * globs.NumLevelRequest);
2132 }
2133
2134 FreeSpectral(vars);
2135
2136 /* -------------------------- */
2137 /* Output of fourier fields */
2138 /* -------------------------- */
2139
2140 if (globs.Type == 60)
2141 {
2142 for (code = 0; code < MaxCodes; code++)
2143 if (vars[code].selected)
2144 {
2145 if (!vars[code].fourier) afterAbort("Code %d not available on fourier space!", code);
2146
2147 nlevel = zaxisInqSize(vars[code].ozaxisID);
2148 for (lindex = 0; lindex < nlevel; lindex++)
2149 {
2150 offset = lindex * globs.DimFC;
2151 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
2152 streamWriteRecord(globs.ostreamID, vars[code].fourier + offset, 0);
2153 }
2154 }
2155
2156 FreeFourier(vars);
2157 return;
2158 }
2159
2160 /* ----------------------- */
2161 /* Output of zonal means */
2162 /* ----------------------- */
2163
2164 if (globs.Type == 61)
2165 {
2166 for (code = 0; code < MaxCodes; code++)
2167 if (vars[code].selected)
2168 {
2169 if (!vars[code].fourier) afterAbort("Code %d not available on zonal mean!", code);
2170
2171 nlevel = zaxisInqSize(vars[code].ozaxisID);
2172 for (lindex = 0; lindex < nlevel; lindex++)
2173 {
2174 offset = lindex * globs.DimFC;
2175 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
2176 streamWriteRecord(globs.ostreamID, vars[code].fourier + offset, vars[code].nmiss);
2177 }
2178 }
2179
2180 FreeFourier(vars);
2181 return;
2182 }
2183
2184 /* ------------------------------ */
2185 /* Transformation to gridpoints */
2186 /* ------------------------------ */
2187
2188 if (globs.Type >= 70)
2189 {
2190 for (code = 0; code < MaxCodes; code++)
2191 if (vars[code].needed && vars[code].fourier)
2192 {
2193 fieldSize = vars[code].plev * globs.DimGP;
2194 if (vars[code].grid == nullptr) vars[code].grid = alloc_dp(fieldSize, FieldName(code, "grid"));
2195
2196 after_FC2GP(vars[code].fourier, vars[code].grid, globs.Latitudes, globs.Longitudes, vars[code].plev, globs.Fouriers);
2197 }
2198 }
2199
2200 FreeFourier(vars);
2201
2202 /* -------------------------------- */
2203 /* Output of pressure level grids */
2204 /* -------------------------------- */
2205
2206 if (globs.Type == 70)
2207 {
2208 for (code = 0; code < MaxCodes; code++)
2209 if (vars[code].selected)
2210 {
2211 nlevel = zaxisInqSize(vars[code].ozaxisID);
2212 for (lindex = 0; lindex < nlevel; lindex++)
2213 {
2214 offset = lindex * globs.DimGP;
2215 streamDefRecord(globs.ostreamID, vars[code].ovarID, lindex);
2216 streamWriteRecord(globs.ostreamID, vars[code].grid + offset, vars[code].nmiss);
2217 }
2218 }
2219
2220 FreeSamp(vars);
2221 FreeGrid(vars);
2222 return;
2223 }
2224 }
2225
2226 void
after_AnalysisAddRecord(const AfterControl * globs,struct Variable * vars,int code,int gridID,int zaxisID,int levelID,size_t nmiss)2227 after_AnalysisAddRecord(const AfterControl *globs, struct Variable *vars, int code, int gridID, int zaxisID, int levelID,
2228 size_t nmiss)
2229 {
2230 size_t fieldSize;
2231 int truncation;
2232
2233 auto gridtype = gridInqType(gridID);
2234 auto leveltype = zaxisInqType(zaxisID);
2235 auto nlevel = zaxisInqSize(zaxisID);
2236 auto gridSize = gridInqSize(gridID);
2237 auto dataSize = gridSize * nlevel;
2238 auto dataOffset = gridSize * levelID;
2239
2240 vars[code].nmiss0 += nmiss;
2241
2242 if (gridtype == GRID_SPECTRAL)
2243 {
2244 vars[code].sfit = true;
2245 vars[code].hlev = globs->NumLevelRequest;
2246 vars[code].plev = globs->NumLevelRequest;
2247 if (nlevel > 1 && leveltype == ZAXIS_PRESSURE)
2248 {
2249 if (code != U_WIND && code != V_WIND)
2250 {
2251 fieldSize = globs->Dim3SP;
2252 if (vars[code].spectral0 == nullptr) vars[code].spectral0 = alloc_dp(fieldSize, FieldName(code, "spectral"));
2253 truncation = gridInqTrunc(gridID);
2254 sp2sp(globs->Field, truncation, vars[code].spectral0 + levelID * globs->DimSP, globs->Truncation);
2255 }
2256 else
2257 {
2258 fieldSize = globs->Dim3SP;
2259 if (vars[code].spectral0 == nullptr) vars[code].spectral0 = alloc_dp(fieldSize, FieldName(code, "spectral"));
2260 array_copy(globs->DimSP, globs->Field, vars[code].spectral0 + levelID * globs->DimSP);
2261 }
2262 }
2263 else
2264 afterAbort("Only pressure level data supported for spectral data!");
2265 }
2266 else
2267 {
2268 if (nlevel > 1 && leveltype == ZAXIS_PRESSURE)
2269 {
2270 vars[code].sfit = true;
2271 fieldSize = globs->Dim3GP;
2272 vars[code].hlev = globs->NumLevelRequest;
2273 vars[code].plev = globs->NumLevelRequest;
2274 if (vars[code].grid0 == nullptr) vars[code].grid0 = alloc_dp(fieldSize, FieldName(code, "grid0"));
2275 array_copy(globs->DimGP, globs->Field, vars[code].grid0 + levelID * globs->DimGP);
2276 }
2277 else
2278 {
2279 vars[code].sfit = false;
2280 fieldSize = globs->DimGP;
2281 vars[code].hlev = 1;
2282 vars[code].plev = 1;
2283 if (vars[code].grid0 == nullptr) vars[code].grid0 = alloc_dp(fieldSize, FieldName(code, "grid0"));
2284 array_copy(globs->DimGP, globs->Field, vars[code].grid0);
2285 }
2286
2287 if (globs->Mean > 0 && (nmiss || vars[code].samp))
2288 {
2289 if (vars[code].samp == nullptr)
2290 {
2291 vars[code].samp = (int *) Malloc(dataSize * sizeof(int));
2292 for (size_t i = 0; i < dataSize; i++) vars[code].samp[i] = globs->MeanCount0;
2293 }
2294
2295 for (size_t i = 0; i < gridSize; i++)
2296 if (IS_NOT_EQUAL(globs->Field[i], vars[code].missval)) vars[code].samp[i + dataOffset]++;
2297 }
2298 }
2299 }
2300
2301 double *
after_get_dataptr(struct Variable * vars,int code,int gridID,int zaxisID,int levelID)2302 after_get_dataptr(struct Variable *vars, int code, int gridID, int zaxisID, int levelID)
2303 {
2304 auto gridtype = gridInqType(gridID);
2305 auto nlevel = zaxisInqSize(zaxisID);
2306 auto gridSize = gridInqSize(gridID);
2307 auto dataSize = gridSize * nlevel;
2308 auto dataOffset = gridSize * levelID;
2309
2310 double *dataptr = nullptr;
2311
2312 if (gridtype == GRID_SPECTRAL)
2313 {
2314 if (vars[code].spectral0 == nullptr) vars[code].spectral0 = alloc_dp(dataSize, FieldName(code, "spectral0"));
2315
2316 dataptr = vars[code].spectral0 + dataOffset;
2317 }
2318 else
2319 {
2320 if (vars[code].hybrid0 == nullptr) vars[code].hybrid0 = alloc_dp(dataSize, FieldName(code, "hybrid0"));
2321
2322 dataptr = vars[code].hybrid0 + dataOffset;
2323 }
2324
2325 return dataptr;
2326 }
2327
2328 void
after_EchamAddRecord(const AfterControl * globs,struct Variable * vars,int code,int gridID,int zaxisID,int levelID,size_t nmiss)2329 after_EchamAddRecord(const AfterControl *globs, struct Variable *vars, int code, int gridID, int zaxisID, int levelID, size_t nmiss)
2330 {
2331 auto gridtype = gridInqType(gridID);
2332 auto leveltype = zaxisInqType(zaxisID);
2333 auto nlevel = zaxisInqSize(zaxisID);
2334 auto gridSize = gridInqSize(gridID);
2335 auto dataSize = gridSize * nlevel;
2336 auto dataOffset = gridSize * levelID;
2337
2338 vars[code].nmiss0 += nmiss;
2339
2340 if (gridtype == GRID_SPECTRAL)
2341 {
2342 vars[code].sfit = true;
2343 vars[code].hlev = nlevel;
2344 vars[code].plev = 1;
2345 if (nlevel > 1 && leveltype == ZAXIS_HYBRID) vars[code].plev = globs->NumLevelRequest;
2346
2347 if (gridInqTrunc(gridID) != globs->Truncation)
2348 afterAbort("Resolution error. Required %d - Found %d", globs->Truncation, gridInqTrunc(gridID));
2349 }
2350 else
2351 {
2352 vars[code].hlev = nlevel;
2353 vars[code].plev = nlevel;
2354 vars[code].sfit = false;
2355 if (nlevel > 1 && leveltype == ZAXIS_HYBRID && (nlevel == globs->NumLevel || nlevel == globs->NumLevel + 1))
2356 {
2357 vars[code].plev = globs->NumLevelRequest;
2358 vars[code].sfit = true;
2359 }
2360
2361 if (globs->Mean > 0 && (nmiss || vars[code].samp))
2362 {
2363 if (vars[code].samp == nullptr)
2364 {
2365 vars[code].samp = (int *) Malloc(dataSize * sizeof(int));
2366 for (size_t i = 0; i < dataSize; i++) vars[code].samp[i] = globs->MeanCount0;
2367 }
2368
2369 double *dataptr = vars[code].hybrid0 + dataOffset;
2370 for (size_t i = 0; i < gridSize; i++)
2371 if (IS_NOT_EQUAL(dataptr[i], vars[code].missval)) vars[code].samp[i + dataOffset]++;
2372 }
2373 }
2374 }
2375
2376 static void
MakeDependencies(struct Variable * vars,int varcode,int depcode)2377 MakeDependencies(struct Variable *vars, int varcode, int depcode)
2378 {
2379 if (vars[varcode].needed && !vars[varcode].detected)
2380 {
2381 vars[depcode].needed = true;
2382 vars[varcode].comp = true;
2383
2384 if (afterDebug) fprintf(stderr, "Needed code %d to compute code %d\n", depcode, varcode);
2385
2386 if (vars[depcode].ivarID == -1)
2387 {
2388 if (depcode == U_WIND)
2389 {
2390 MakeDependencies(vars, U_WIND, DIVERGENCE);
2391 MakeDependencies(vars, U_WIND, VORTICITY);
2392 }
2393 if (depcode == V_WIND)
2394 {
2395 MakeDependencies(vars, V_WIND, DIVERGENCE);
2396 MakeDependencies(vars, V_WIND, VORTICITY);
2397 }
2398 }
2399
2400 if (vars[varcode].ivarID == -1)
2401 {
2402 if (vars[depcode].ivarID == -1)
2403 {
2404 afterAbort("code %d undefined, needed to compute code %d", depcode, varcode);
2405 }
2406 else
2407 {
2408 vars[varcode].ivarID = vars[depcode].ivarID;
2409 vars[varcode].igridID = vars[depcode].igridID;
2410 vars[varcode].ogridID = vars[depcode].ogridID;
2411 vars[varcode].izaxisID = vars[depcode].izaxisID;
2412 vars[varcode].ozaxisID = vars[depcode].ozaxisID;
2413 }
2414 }
2415 }
2416 }
2417
2418 static void
CheckDependencies(struct Variable * vars,int analysisdata)2419 CheckDependencies(struct Variable *vars, int analysisdata)
2420 {
2421 MakeDependencies(vars, VELOPOT, U_WIND);
2422 MakeDependencies(vars, VELOPOT, V_WIND);
2423 MakeDependencies(vars, VELOPOT, VORTICITY);
2424 MakeDependencies(vars, VELOPOT, DIVERGENCE);
2425
2426 MakeDependencies(vars, STREAM, U_WIND);
2427 MakeDependencies(vars, STREAM, V_WIND);
2428 MakeDependencies(vars, STREAM, VORTICITY);
2429 MakeDependencies(vars, STREAM, DIVERGENCE);
2430
2431 MakeDependencies(vars, VORTICITY, U_WIND);
2432 MakeDependencies(vars, VORTICITY, V_WIND);
2433
2434 MakeDependencies(vars, DIVERGENCE, U_WIND);
2435 MakeDependencies(vars, DIVERGENCE, V_WIND);
2436
2437 MakeDependencies(vars, U_WIND, VORTICITY);
2438 MakeDependencies(vars, U_WIND, DIVERGENCE);
2439 MakeDependencies(vars, U_WIND, V_WIND);
2440
2441 MakeDependencies(vars, V_WIND, VORTICITY);
2442 MakeDependencies(vars, V_WIND, DIVERGENCE);
2443 MakeDependencies(vars, V_WIND, U_WIND);
2444
2445 MakeDependencies(vars, WINDSPEED, U_WIND);
2446 MakeDependencies(vars, WINDSPEED, V_WIND);
2447
2448 if (analysisdata)
2449 {
2450 MakeDependencies(vars, RHUMIDITY, HUMIDITY);
2451 MakeDependencies(vars, RHUMIDITY, TEMPERATURE);
2452 MakeDependencies(vars, HUMIDITY, RHUMIDITY);
2453 MakeDependencies(vars, HUMIDITY, TEMPERATURE);
2454 MakeDependencies(vars, GEOPOTHEIGHT, GEOPOTENTIAL);
2455 }
2456 else
2457 {
2458 MakeDependencies(vars, THETAF, TEMPERATURE);
2459 MakeDependencies(vars, SLP, TEMPERATURE);
2460 }
2461
2462 MakeDependencies(vars, SW_ATM, 176);
2463 MakeDependencies(vars, SW_ATM, 178);
2464 MakeDependencies(vars, LW_ATM, 177);
2465 MakeDependencies(vars, LW_ATM, 179);
2466 MakeDependencies(vars, NET_ATM, 176);
2467 MakeDependencies(vars, NET_ATM, 177);
2468 MakeDependencies(vars, NET_ATM, 178);
2469 MakeDependencies(vars, NET_ATM, 179);
2470
2471 MakeDependencies(vars, SW_BOT_CLF, 176);
2472 MakeDependencies(vars, SW_BOT_CLF, 185);
2473 MakeDependencies(vars, LW_BOT_CLF, 177);
2474 MakeDependencies(vars, LW_BOT_CLF, 186);
2475 MakeDependencies(vars, SW_TOP_CLF, 178);
2476 MakeDependencies(vars, SW_TOP_CLF, 187);
2477 MakeDependencies(vars, LW_TOP_CLF, 179);
2478 MakeDependencies(vars, LW_TOP_CLF, 188);
2479 MakeDependencies(vars, NET_TOP_CLF, 178);
2480 MakeDependencies(vars, NET_TOP_CLF, 179);
2481 MakeDependencies(vars, NET_TOP_CLF, 187);
2482 MakeDependencies(vars, NET_TOP_CLF, 188);
2483
2484 MakeDependencies(vars, ALL_WATER, 222);
2485 MakeDependencies(vars, LOW_WATER, 222);
2486 MakeDependencies(vars, MID_WATER, 222);
2487 MakeDependencies(vars, HIH_WATER, 222);
2488
2489 MakeDependencies(vars, LOW_CLOUD, 223);
2490 MakeDependencies(vars, MID_CLOUD, 223);
2491 MakeDependencies(vars, HIH_CLOUD, 223);
2492
2493 if (vars[LOW_CLOUD].comp || vars[MID_CLOUD].comp || vars[HIH_CLOUD].comp)
2494 {
2495 static int zaxisID = -999;
2496 if (zaxisID == -999) zaxisID = zaxisCreate(ZAXIS_SURFACE, 1);
2497
2498 vars[LOW_CLOUD].izaxisID = zaxisID;
2499 vars[LOW_CLOUD].ozaxisID = zaxisID;
2500 vars[MID_CLOUD].izaxisID = zaxisID;
2501 vars[MID_CLOUD].ozaxisID = zaxisID;
2502 vars[HIH_CLOUD].izaxisID = zaxisID;
2503 vars[HIH_CLOUD].ozaxisID = zaxisID;
2504 }
2505 }
2506
2507 void
after_AnalysisDependencies(struct Variable * vars,int ncodes)2508 after_AnalysisDependencies(struct Variable *vars, int ncodes)
2509 {
2510 int code;
2511
2512 for (code = 0; code < ncodes; code++) vars[code].needed = vars[code].selected;
2513
2514 MakeDependencies(vars, PS, LNPS);
2515
2516 CheckDependencies(vars, 1);
2517 }
2518
2519 void
after_EchamDependencies(struct Variable * vars,int ncodes,int type,int source)2520 after_EchamDependencies(struct Variable *vars, int ncodes, int type, int source)
2521 {
2522 int code;
2523
2524 for (code = 0; code < ncodes; code++) vars[code].needed = vars[code].selected;
2525
2526 for (code = 0; code < ncodes; code++)
2527 if (vars[code].detected == false) vars[code].ivarID = -1;
2528
2529 if (type >= 50)
2530 {
2531 vars[U_WIND].needed |= vars[DIVERGENCE].needed;
2532 vars[V_WIND].needed |= vars[DIVERGENCE].needed;
2533 vars[U_WIND].needed |= vars[VORTICITY].needed;
2534 vars[V_WIND].needed |= vars[VORTICITY].needed;
2535 vars[U_WIND].needed |= vars[VELOPOT].needed;
2536 vars[V_WIND].needed |= vars[VELOPOT].needed;
2537 vars[U_WIND].needed |= vars[STREAM].needed;
2538 vars[V_WIND].needed |= vars[STREAM].needed;
2539 }
2540
2541 if (type >= 30) vars[LNPS].needed = true;
2542
2543 if (type >= 20)
2544 {
2545 MakeDependencies(vars, THETAF, LNPS);
2546 MakeDependencies(vars, SLP, LNPS);
2547 MakeDependencies(vars, SLP, GEOPOTENTIAL);
2548 /*
2549 MakeDependencies(vars, SLP, HALF_PRESS);
2550 MakeDependencies(vars, SLP, FULL_PRESS);
2551 */
2552 MakeDependencies(vars, RHUMIDITY, TEMPERATURE);
2553 MakeDependencies(vars, RHUMIDITY, HUMIDITY);
2554 MakeDependencies(vars, RHUMIDITY, LNPS);
2555 MakeDependencies(vars, GEOPOTHEIGHT, TEMPERATURE);
2556 MakeDependencies(vars, GEOPOTHEIGHT, HUMIDITY);
2557 MakeDependencies(vars, GEOPOTHEIGHT, LNPS);
2558 MakeDependencies(vars, OMEGA, DIVERGENCE);
2559 MakeDependencies(vars, OMEGA, U_WIND);
2560 MakeDependencies(vars, OMEGA, V_WIND);
2561 MakeDependencies(vars, OMEGA, LNPS);
2562 MakeDependencies(vars, OMEGA, DPSDX);
2563 MakeDependencies(vars, OMEGA, DPSDY);
2564 }
2565
2566 MakeDependencies(vars, DPSDX, LNPS);
2567 MakeDependencies(vars, HALF_PRESS, LNPS);
2568 MakeDependencies(vars, PS, LNPS);
2569
2570 MakeDependencies(vars, SURF_RUNOFF, 142);
2571 MakeDependencies(vars, SURF_RUNOFF, 143);
2572 MakeDependencies(vars, SURF_RUNOFF, 182);
2573 MakeDependencies(vars, SURF_RUNOFF, 221); /* snow depth change */
2574
2575 MakeDependencies(vars, THETAF, TS);
2576
2577 MakeDependencies(vars, FRESH_WATER, 142);
2578 MakeDependencies(vars, FRESH_WATER, 143);
2579 MakeDependencies(vars, FRESH_WATER, 182);
2580
2581 MakeDependencies(vars, PRECIP, 142);
2582 MakeDependencies(vars, PRECIP, 143);
2583
2584 if (source != S_ECHAM5)
2585 {
2586 MakeDependencies(vars, NET_WATER, 142);
2587 MakeDependencies(vars, NET_WATER, 143);
2588 MakeDependencies(vars, NET_WATER, 160);
2589 MakeDependencies(vars, NET_WATER, 182);
2590
2591 MakeDependencies(vars, NET_TOP, 178);
2592 MakeDependencies(vars, NET_TOP, 179);
2593
2594 MakeDependencies(vars, NET_BOT, 176);
2595 MakeDependencies(vars, NET_BOT, 177);
2596
2597 MakeDependencies(vars, NET_HEAT, 146);
2598 MakeDependencies(vars, NET_HEAT, 147);
2599 MakeDependencies(vars, NET_HEAT, 176);
2600 MakeDependencies(vars, NET_HEAT, 177);
2601 MakeDependencies(vars, NET_HEAT, 218);
2602 /*
2603 if ( source == S_ECHAM5 )
2604 {
2605 MakeDependencies(vars, NET_HEAT, 206);
2606 MakeDependencies(vars, NET_HEAT, 208);
2607 MakeDependencies(vars, NET_HEAT, 209);
2608 }
2609 else
2610 */
2611 {
2612 MakeDependencies(vars, NET_HEAT, 220);
2613 }
2614 }
2615
2616 MakeDependencies(vars, SW_CLF, 178);
2617 MakeDependencies(vars, SW_CLF, 224);
2618
2619 MakeDependencies(vars, LW_CLF, 179);
2620 MakeDependencies(vars, LW_CLF, 225);
2621
2622 MakeDependencies(vars, NET_CLF, 178);
2623 MakeDependencies(vars, NET_CLF, 179);
2624 MakeDependencies(vars, NET_CLF, 224);
2625 MakeDependencies(vars, NET_CLF, 225);
2626
2627 if (vars[DPSDX].needed || vars[DPSDY].needed || vars[GEOPOTHEIGHT].comp || vars[SLP].comp || vars[THETAF].needed
2628 || vars[HALF_PRESS].needed || vars[RHUMIDITY].comp || vars[OMEGA].comp || type >= 30)
2629 vars[PS_PROG].comp = true;
2630
2631 CheckDependencies(vars, 0);
2632 }
2633
2634
2635 void
after_legini_setup(AfterControl & globs,struct Variable * vars)2636 after_legini_setup(AfterControl &globs, struct Variable *vars)
2637 {
2638 const long ntr = globs.Truncation;
2639 const long nlat = globs.Latitudes;
2640 const long dimsp = (ntr + 1) * (ntr + 2);
2641 const long pdim = (dimsp / 2) * nlat;
2642
2643 globs.poli = (double *) Malloc(pdim * sizeof(double));
2644
2645 if (!globs.AnalysisData)
2646 {
2647 if (globs.Type >= 20) globs.pold = (double *) Malloc(pdim * sizeof(double));
2648 if (vars[DPSDY].needed) globs.pdev = (double *) Malloc(pdim * sizeof(double));
2649 }
2650
2651 if ((vars[DIVERGENCE].needed || vars[VORTICITY].needed || vars[VELOPOT].needed || vars[STREAM].needed) && globs.Type > 20)
2652 {
2653 globs.pol2 = (double *) Malloc(pdim * sizeof(double));
2654 globs.pol3 = (double *) Malloc(pdim * sizeof(double));
2655 }
2656
2657 if (globs.AnalysisData && (globs.Type == 70) && globs.Gaussian && !globs.Spectral)
2658 {
2659 if (globs.poli)
2660 {
2661 Free(globs.poli);
2662 globs.poli = nullptr;
2663 }
2664 if (globs.pol2)
2665 {
2666 Free(globs.pol2);
2667 globs.pol2 = nullptr;
2668 }
2669 if (globs.pol3)
2670 {
2671 Free(globs.pol3);
2672 globs.pol3 = nullptr;
2673 }
2674 return;
2675 }
2676
2677 after_legini_full(ntr, nlat, globs.poli, globs.pold, globs.pdev, globs.pol2, globs.pol3, globs.coslat);
2678
2679 for (long jgl = 0; jgl < nlat; ++jgl) globs.rcoslat[jgl] = 1.0 / globs.coslat[jgl];
2680
2681 for (long jgl = 0; jgl < nlat; ++jgl) globs.DerivationFactor[jgl] = globs.rcoslat[jgl] / PlanetRadius;
2682 }
2683