Ticket 686: WRWP
[baltrad-wrwp.git] / lib / wrwp.c
1 /* --------------------------------------------------------------------
2 Copyright (C) 2011 Swedish Meteorological and Hydrological Institute, SMHI
3
4 This is free software: you can redistribute it and/or modify
5 it under the terms of the GNU Lesser General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
8
9 This software is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 GNU Lesser General Public License for more details.
13
14 You should have received a copy of the GNU Lesser General Public License
15 along with HLHDF.  If not, see <http://www.gnu.org/licenses/>.
16 ------------------------------------------------------------------------*/
17
18 /** Function for deriving weather radar wind and reflectivity profiles
19  *
20  * @file
21  * @author Gunther Haase, SMHI
22  * @date 2013-02-06
23  *
24  * @co-author Ulf E. Nordh, SMHI
25  * @date 2017-02-23, started overhaul of the code to achieve better
26  * resemblance with N2 and requirements from customer E-profile
27  * 
28  */
29
30 #include "wrwp.h"
31 #include "vertical_profile.h"
32 #include "rave_debug.h"
33 #include "rave_alloc.h"
34 #include <string.h>
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include "rave_attribute.h"
38 /**
39  * Represents one wrwp generator
40  */
41 struct _Wrwp_t {
42   RAVE_OBJECT_HEAD /** Always on top */
43   int dz; /**< Height interval for deriving a profile [m] */
44   int hmax; /**< Maximum height of the profile [m] */
45   int dmin; /**< Minimum distance for deriving a profile [m] */
46   int dmax; /**< Maximum distance for deriving a profile [m]*/
47   double emin; /**< Minimum elevation angle [deg] */
48   double vmin; /**< Radial velocity threshold [m/s] */
49   double nodata_VP; /**< Nodata value for vertical profile */
50   double gain_VP; /**< Gain for VP fields */
51   double offset_VP; /**< Offset for VP fields */
52   double undetect_VP; /**<Undetect for VP fields */
53 };
54
55 /*@{ Private functions */
56 /**
57  * Constructor
58  */
59 static int Wrwp_constructor(RaveCoreObject* obj)
60 {
61   Wrwp_t* wrwp = (Wrwp_t*)obj;
62   wrwp->hmax = HMAX;
63   wrwp->dz = DZ;
64   wrwp->dmin = DMIN;
65   wrwp->dmax = DMAX;
66   wrwp->emin = EMIN;
67   wrwp->vmin = VMIN;
68   wrwp->nodata_VP = NODATA_VP;
69   wrwp->gain_VP = GAIN_VP;
70   wrwp->offset_VP = OFFSET_VP;
71   wrwp->undetect_VP = UNDETECT_VP;
72   return 1;
73 }
74
75 /**
76  * Destructor
77  */
78 static void Wrwp_destructor(RaveCoreObject* obj)
79 {
80 }
81
82 static int WrwpInternal_findAndAddAttribute(VerticalProfile_t* vp, PolarVolume_t* pvol, const char* name, double minSelAng)
83 {
84   int nscans = PolarVolume_getNumberOfScans(pvol);
85   int i = 0;
86   int found = 0;
87   double elangle = 0.0;
88   
89   for (i = 0; i < nscans && found == 0; i++) {
90     PolarScan_t* scan = PolarVolume_getScan(pvol, i);
91     if (scan != NULL && PolarScan_hasAttribute(scan, name)) {
92         elangle = PolarScan_getElangle(scan);
93         
94         /* Filter with respect to the selected min elangle
95            that is given in the web GUI route for WRWP generation */
96         if ((elangle * RAD2DEG) >= minSelAng) {
97           RaveAttribute_t* attr = PolarScan_getAttribute(scan, name);
98           VerticalProfile_addAttribute(vp, attr);
99           found = 1;
100           RAVE_OBJECT_RELEASE(attr);
101         }
102     }
103     RAVE_OBJECT_RELEASE(scan);
104   }
105   return found;
106 }
107
108 static int WrwpInternal_addIntAttribute(VerticalProfile_t* vp, const char* name, int value)
109 {
110   RaveAttribute_t* attr = RaveAttributeHelp_createLong(name, value);
111   int result = 0;
112   if (attr != NULL) {
113     result = VerticalProfile_addAttribute(vp, attr);
114   }
115   RAVE_OBJECT_RELEASE(attr);
116   return result;
117 }
118
119 /* Function that adds various quantities under a field's what in order to
120    better resemble the function existing in vertical profiles from N2 */
121 static int WrwpInternal_addDoubleAttr2Field(RaveField_t* field, const char* name, double quantity)
122 {
123   int result = 0;
124   RaveAttribute_t* attr = RaveAttributeHelp_createDouble(name, quantity);
125   result = RaveField_addAttribute(field, attr);
126   if (attr == NULL || !result) {
127     RAVE_ERROR0("Failed to add what/quantity attribute to field");
128     goto done;
129   }
130   RAVE_OBJECT_RELEASE(attr);
131 done:
132   return result;
133 }
134
135 /*@} End of Private functions */
136
137 /*@{ Interface functions */
138 void Wrwp_setDZ(Wrwp_t* self, int dz)
139 {
140   RAVE_ASSERT((self != NULL), "self == NULL");
141   self->dz = dz;
142 }
143
144 int Wrwp_getDZ(Wrwp_t* self)
145 {
146   RAVE_ASSERT((self != NULL), "self == NULL");
147   return self->dz;
148 }
149
150 void Wrwp_setNODATA_VP(Wrwp_t* self, int nodata_VP)
151 {
152   RAVE_ASSERT((self != NULL), "self == NULL");
153   self->nodata_VP = nodata_VP;
154 }
155
156 int Wrwp_getNODATA_VP(Wrwp_t* self)
157 {
158   RAVE_ASSERT((self != NULL), "self == NULL");
159   return self->nodata_VP;
160 }
161
162 void Wrwp_setUNDETECT_VP(Wrwp_t* self, int undetect_VP)
163 {
164   RAVE_ASSERT((self != NULL), "self == NULL");
165   self->undetect_VP = undetect_VP;
166 }
167
168 int Wrwp_getUNDETECT_VP(Wrwp_t* self)
169 {
170   RAVE_ASSERT((self != NULL), "self == NULL");
171   return self->undetect_VP;
172 }
173
174 void Wrwp_setGAIN_VP(Wrwp_t* self, int gain_VP)
175 {
176   RAVE_ASSERT((self != NULL), "self == NULL");
177   self->gain_VP = gain_VP;
178 }
179
180 int Wrwp_getGAIN_VP(Wrwp_t* self)
181 {
182   RAVE_ASSERT((self != NULL), "self == NULL");
183   return self->gain_VP;
184 }
185
186 void Wrwp_setOFFSET_VP(Wrwp_t* self, int offset_VP)
187 {
188   RAVE_ASSERT((self != NULL), "self == NULL");
189   self->offset_VP = offset_VP;
190 }
191
192 int Wrwp_getOFFSET_VP(Wrwp_t* self)
193 {
194   RAVE_ASSERT((self != NULL), "self == NULL");
195   return self->offset_VP;
196 }
197
198 void Wrwp_setHMAX(Wrwp_t* self, int hmax)
199 {
200   RAVE_ASSERT((self != NULL), "self == NULL");
201   self->hmax = hmax;
202 }
203
204 int Wrwp_getHMAX(Wrwp_t* self)
205 {
206   RAVE_ASSERT((self != NULL), "self == NULL");
207   return self->hmax;
208 }
209
210 void Wrwp_setDMIN(Wrwp_t* self, int dmin)
211 {
212   RAVE_ASSERT((self != NULL), "self == NULL");
213   self->dmin = dmin;
214 }
215
216 int Wrwp_getDMIN(Wrwp_t* self)
217 {
218   RAVE_ASSERT((self != NULL), "self == NULL");
219   return self->dmin;
220 }
221
222 void Wrwp_setDMAX(Wrwp_t* self, int dmax)
223 {
224   RAVE_ASSERT((self != NULL), "self == NULL");
225   self->dmax = dmax;
226 }
227
228 int Wrwp_getDMAX(Wrwp_t* self)
229 {
230   RAVE_ASSERT((self != NULL), "self == NULL");
231   return self->dmax;
232 }
233
234 void Wrwp_setEMIN(Wrwp_t* self, double emin)
235 {
236   RAVE_ASSERT((self != NULL), "self == NULL");
237   self->emin = emin;
238 }
239
240 double Wrwp_getEMIN(Wrwp_t* self)
241 {
242   RAVE_ASSERT((self != NULL), "self == NULL");
243   return self->emin;
244 }
245
246 void Wrwp_setVMIN(Wrwp_t* self, double vmin)
247 {
248   RAVE_ASSERT((self != NULL), "self == NULL");
249   self->vmin = vmin;
250 }
251
252 double Wrwp_getVMIN(Wrwp_t* self)
253 {
254   RAVE_ASSERT((self != NULL), "self == NULL");
255   return self->vmin;
256 }
257
258 VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj) {
259   VerticalProfile_t* result = NULL;
260   PolarScan_t* scan = NULL;
261   PolarScanParam_t* vrad = NULL;
262   PolarScanParam_t* dbz = NULL;
263   PolarNavigator_t* polnav = NULL;
264   int nrhs = NRHS, lda = LDA, ldb = LDB;
265   int nscans = 0, nbins = 0, nrays, nv, nz, i, iz, is, ib, ir;
266   int findEarliestTime;
267   int findLatestTime;
268   int firstInit;
269
270   double rscale, elangle, gain, offset, nodata, undetect, val;
271   double d, h, elangleForThisScan;
272   double alpha, beta, gamma, vvel, vdir, vstd, zsum, zmean, zstd;
273   double centerOfLayer, u_wnd_comp, v_wnd_comp, vdir_rad;
274   int ysize = 0, yindex = 0;
275   
276   const char* starttime = NULL;
277   const char* endtime = NULL;
278   const char* startdate = NULL;
279   const char* enddate = NULL;
280   const char* product = "VP";
281   
282   const char* startTimeOfThisScan = NULL;
283   const char* endTimeOfThisScan = NULL;
284   const char* startDateOfThisScan = NULL;
285   const char* endDateOfThisScan = NULL;
286     
287   const char* starttimeFirst = NULL;
288   const char* endtimeFirst = NULL;
289   const char* startdateFirst = NULL;
290   const char* enddateFirst = NULL;
291       
292   const char* starttimeNext = NULL;
293   const char* endtimeNext = NULL;
294   const char* startdateNext = NULL;
295   const char* enddateNext = NULL;
296   
297   char* startDateTimeStrFirst = NULL;
298   char* endDateTimeStrFirst = NULL;
299   
300   char* startDateTimeStrNext = NULL;
301   char* endDateTimeStrNext = NULL;
302    
303   char starttimeArrFirst[10];
304   char endtimeArrFirst[10]; 
305   char startdateArrFirst[10];
306   char enddateArrFirst[10];  
307  
308   char starttimeArrNext[10];
309   char endtimeArrNext[10];
310   char startdateArrNext[10];
311   char enddateArrNext[10];
312
313   char* malfuncString = NULL;
314    
315   /* The four field below are the ones fulfilling the requirements from the user E-profiles */
316   RaveField_t *nv_field = NULL, *HGHT_field = NULL;
317   RaveField_t *UWND_field = NULL, *VWND_field = NULL;
318   
319   /* The following fields are removed from the VP but kept in the code
320      for convenience reasons, add to initialization when needed:     
321      *RaveField_t *ff_field = NULL, *ff_dev_field = NULL, *dd_field = NULL;
322       RaveField_t *dbzh_field = NULL, *dbzh_dev_field = NULL, *nz_field = NULL; */
323
324   RAVE_ASSERT((self != NULL), "self == NULL");
325
326   nv_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
327   HGHT_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
328   UWND_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
329   VWND_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
330     
331   /* Field defs kept in the code for convenience reasons, add when needed:
332   ff_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
333   ff_dev_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
334   dd_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
335   dbzh_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
336   dbzh_dev_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
337   nz_field = RAVE_OBJECT_NEW(&RaveField_TYPE); */
338
339   if (nv_field == NULL || HGHT_field == NULL || UWND_field == NULL ||
340       VWND_field == NULL) {
341     /* Removed, add when needed:ff_field == NULL || ff_dev_field == NULL || 
342        dd_field == NULL || dbzh_field == NULL || dbzh_dev_field == NULL || 
343        nz_field == NULL || */
344     RAVE_ERROR0("Failed to allocate memory for the resulting vp fields");
345     goto done;
346   }
347   ysize = self->hmax / self->dz;
348   if (!RaveField_createData(nv_field, 1, ysize, RaveDataType_INT) || 
349       !RaveField_createData(HGHT_field, 1, ysize, RaveDataType_DOUBLE) ||
350       !RaveField_createData(UWND_field, 1, ysize, RaveDataType_DOUBLE) ||
351       !RaveField_createData(VWND_field, 1, ysize, RaveDataType_DOUBLE)) {
352       /* Add when needed: !RaveField_createData(ff_field, 1, ysize, RaveDataType_DOUBLE) ||
353       !RaveField_createData(ff_dev_field, 1, ysize, RaveDataType_DOUBLE) ||
354       !RaveField_createData(dd_field, 1, ysize, RaveDataType_DOUBLE) ||
355       !RaveField_createData(dbzh_field, 1, ysize, RaveDataType_DOUBLE) ||
356       !RaveField_createData(dbzh_dev_field, 1, ysize, RaveDataType_DOUBLE) ||
357       !RaveField_createData(nz_field, 1, ysize, RaveDataType_INT) || */
358     RAVE_ERROR0("Failed to allocate arrays for the resulting vp fields");
359     goto done;
360   }
361
362         polnav = RAVE_OBJECT_NEW(&PolarNavigator_TYPE);
363         PolarNavigator_setLat0(polnav, PolarVolume_getLatitude(inobj));
364         PolarNavigator_setLon0(polnav, PolarVolume_getLongitude(inobj));
365         PolarNavigator_setAlt0(polnav, PolarVolume_getHeight(inobj));
366             
367         nscans = PolarVolume_getNumberOfScans (inobj);
368     
369
370         // We use yindex for filling in the arrays even though we loop to hmax...
371         yindex = 0;
372
373
374         // loop over atmospheric layers
375   for (iz = 0; iz < self->hmax; iz += self->dz) {
376     /* allocate memory and initialize with zeros */
377     double *A = RAVE_CALLOC((size_t)(NOR*NOC), sizeof (double));
378     double *b = RAVE_CALLOC((size_t)(NOR), sizeof (double));
379     double *v = RAVE_CALLOC((size_t)(NOR), sizeof (double));
380     double *z = RAVE_CALLOC((size_t)(NOR), sizeof (double));
381     double *az = RAVE_CALLOC((size_t)(NOR), sizeof (double));
382
383     vdir = -9999.0;
384     vvel = -9999.0;
385     vstd = 0.0;
386     zsum = 0.0;
387     zmean = -9999.0;
388     zstd = 0.0;
389     nv = 0;
390     nz = 0;
391     
392     /* Define the center height of each verical layer, this will later
393        become the HGHT array */
394     centerOfLayer = iz + (self->dz / 2.0);
395     firstInit = 0;
396     
397     // loop over scans
398     for (is = 0; is < nscans; is++) {      
399       scan = PolarVolume_getScan(inobj, is);
400       nbins = PolarScan_getNbins(scan);
401       nrays = PolarScan_getNrays(scan);
402       rscale = PolarScan_getRscale(scan);
403       elangleForThisScan = PolarScan_getElangle(scan);
404       elangle = elangleForThisScan;
405       startTimeOfThisScan = PolarScan_getStartTime(scan);
406       endTimeOfThisScan = PolarScan_getEndTime(scan);
407       startDateOfThisScan = PolarScan_getStartDate(scan);
408       endDateOfThisScan = PolarScan_getEndDate(scan);
409       
410       RaveAttribute_t* attr = PolarScan_getAttribute(scan, "how/malfunc");
411       if (attr != NULL) {
412         RaveAttribute_getString(attr, &malfuncString); /* Set the malfuncString if attr is not NULL */
413       }
414       if (malfuncString == NULL || strcmp(malfuncString, "False") == 0) { /* Assuming malfuncString = NULL means no malfunc */
415                   
416         if (elangleForThisScan * RAD2DEG >= self->emin && firstInit == 0) {
417           /* Initialize using the first scan and define 2 strings for the combined datetime */
418           starttimeFirst = startTimeOfThisScan;
419           endtimeFirst = endTimeOfThisScan;
420           startdateFirst = startDateOfThisScan;
421           enddateFirst = endDateOfThisScan;
422           firstInit = 1;
423         }
424         
425         if (elangleForThisScan * RAD2DEG >= self->emin && firstInit == 1 && is > 0) {
426           starttimeNext = startTimeOfThisScan;
427           endtimeNext = endTimeOfThisScan;
428           startdateNext = startDateOfThisScan;
429           enddateNext = endDateOfThisScan;
430
431           strcpy(starttimeArrFirst, starttimeFirst);
432           strcpy(endtimeArrFirst, endtimeFirst);
433           strcpy(startdateArrFirst, startdateFirst);
434           strcpy(enddateArrFirst, enddateFirst);
435
436           startDateTimeStrFirst = strcat(startdateArrFirst, starttimeArrFirst);
437           endDateTimeStrFirst = strcat(enddateArrFirst, endtimeArrFirst);
438
439           strcpy(starttimeArrNext, starttimeNext);
440           strcpy(endtimeArrNext, endtimeNext);
441           strcpy(startdateArrNext, startdateNext);
442           strcpy(enddateArrNext, enddateNext);
443
444           startDateTimeStrNext = strcat(startdateArrNext, starttimeArrNext);
445           endDateTimeStrNext = strcat(enddateArrNext, endtimeArrNext);
446
447           /* Find the earliest and latest datetime's and replace the initial ones
448              if they represent an earlier starttime/startdate and a later endtime/enddate */
449           findEarliestTime = strcmp(startDateTimeStrFirst, startDateTimeStrNext);
450           findLatestTime = strcmp(endDateTimeStrFirst, endDateTimeStrNext);
451
452           /* The starttime and startdate */
453           if (findEarliestTime > 0) {
454             starttimeFirst = starttimeNext;
455             startdateFirst = startdateNext;
456           }
457
458           /* The endtime and enddate */
459           if (findLatestTime < 0) {
460             endtimeFirst = endtimeNext;
461             enddateFirst = enddateNext;
462           }             
463         }
464                      
465         // radial wind scans
466         if (PolarScan_hasParameter(scan, "VRAD") || PolarScan_hasParameter(scan, "VRADH")) {
467           if (PolarScan_hasParameter(scan, "VRAD")) {
468             vrad = PolarScan_getParameter(scan, "VRAD");
469           } else {
470             vrad = PolarScan_getParameter(scan, "VRADH");
471           } 
472           gain = PolarScanParam_getGain(vrad);
473           offset = PolarScanParam_getOffset(vrad);
474           nodata = PolarScanParam_getNodata(vrad);
475           undetect = PolarScanParam_getUndetect(vrad);
476
477           for (ir = 0; ir < nrays; ir++) {
478             for (ib = 0; ib < nbins; ib++) {
479               PolarNavigator_reToDh(polnav, (ib+0.5)*rscale, elangle, &d, &h);
480               PolarScanParam_getValue(vrad, ib, ir, &val);
481
482               if ((h >= iz) &&
483                   (h < iz + self->dz) &&
484                   (d >= self->dmin) &&
485                   (d <= self->dmax) &&
486                   (elangle * RAD2DEG >= self->emin) &&
487                   (val != nodata) &&
488                   (val != undetect) &&
489                   (abs(offset + gain * val) >= self->vmin)) {
490                 *(v+nv) = offset+gain*val;
491                 *(az+nv) = 360./nrays*ir*DEG2RAD;
492                 *(A+nv*NOC) = sin(*(az+nv));
493                 *(A+nv*NOC+1) = cos(*(az+nv));
494                 *(A+nv*NOC+2) = 1;
495                 *(b+nv) = *(v+nv);
496                 nv = nv+1;
497               }
498             }
499           }
500           RAVE_OBJECT_RELEASE(vrad);
501         }
502
503         // reflectivity scans
504         if (PolarScan_hasParameter(scan, "DBZH")) {
505           dbz = PolarScan_getParameter(scan, "DBZH");
506           gain = PolarScanParam_getGain(dbz);
507           offset = PolarScanParam_getOffset(dbz);
508           nodata = PolarScanParam_getNodata(dbz);
509           undetect = PolarScanParam_getUndetect(dbz);
510
511           for (ir = 0; ir < nrays; ir++) {
512             for (ib = 0; ib < nbins; ib++) {
513               PolarNavigator_reToDh(polnav, (ib+0.5)*rscale, elangle, &d, &h);
514               PolarScanParam_getValue (dbz, ib, ir, &val);
515               if ((h >= iz) &&
516                   (h < iz+self->dz) &&
517                   (d >= self->dmin) &&
518                   (d <= self->dmax) &&
519                   (elangle*RAD2DEG >= self->emin) &&
520                   (val != nodata) &&
521                   (val != undetect)) {
522                 *(z+nz) = dBZ2Z(offset+gain*val);
523                 zsum = zsum + *(z+nz);
524                 nz = nz+1;
525               }
526             }
527           }
528           RAVE_OBJECT_RELEASE(dbz);
529         }
530         RAVE_OBJECT_RELEASE(scan);
531       }
532     }
533
534     // radial wind calculations
535     if (nv>3) {
536       //***************************************************************
537       // fitting: y = gamma+alpha*sin(x+beta)                         *
538       // alpha -> amplitude                                           *
539       // beta -> phase shift                                          *
540       // gamma -> consider an y-shift due to the terminal velocity of *
541       //          falling rain drops                                  *
542       //***************************************************************
543
544       /* QR decomposition */
545       /*info = */LAPACKE_dgels(LAPACK_ROW_MAJOR, 'N', NOR, NOC, nrhs, A, lda, b, ldb);
546
547       /* parameter of the wind model */
548       alpha = sqrt(pow(*(b),2) + pow(*(b+1),2));
549       beta = atan2(*(b+1), *b);
550       gamma = *(b+2);
551
552       /* wind velocity */
553       vvel = alpha;
554
555       /* wind direction */
556       vdir = 0;
557       if (alpha < 0) {
558         vdir = (PI/2-beta) * RAD2DEG;
559       } else if (alpha > 0) {
560         vdir = (3*PI/2-beta) * RAD2DEG;
561       }
562       if (vdir < 0) {
563         vdir = vdir + 360;
564       } else if (vdir > 360) {
565         vdir = vdir - 360;
566       }
567
568       /* standard deviation of the wind velocity */
569       for (i=0; i<nv; i++) {
570         vstd = vstd + pow (*(v+i) - (gamma+alpha*sin (*(az+i)+beta)),2);
571       }
572       vstd = sqrt (vstd/(nv-1));
573     }
574
575     // reflectivity calculations
576     if (nz > 1) {
577       /* standard deviation of reflectivity */
578       for (i = 0; i < nz; i++) {
579         zstd = zstd + pow(*(z+i) - (zsum/nz),2);
580       }
581       zmean = Z2dBZ(zsum/nz);
582       zstd = sqrt(zstd/(nz-1));
583       zstd = Z2dBZ(zstd);
584     }
585
586     /* Calculate the x-component (East) and y-component (North) of the wind 
587        velocity using the wind direction and the magnitude of the wind velocity */
588     vdir_rad = vdir * DEG2RAD;
589     u_wnd_comp = vvel * cos(vdir_rad);
590     v_wnd_comp = vvel * sin(vdir_rad);
591     
592     if ((nv < NMIN) || (nz < NMIN)) {
593       /* Add when needed: RaveField_setValue(ff_field, 0, yindex, self->nodata_VP);
594       RaveField_setValue(ff_dev_field, 0, yindex, self->nodata_VP);
595       RaveField_setValue(dd_field, 0, yindex, self->nodata_VP);
596       RaveField_setValue(dbzh_field, 0, yindex, self->nodata_VP);
597       RaveField_setValue(dbzh_dev_field, 0, yindex, self->nodata_VP);
598       RaveField_setValue(nz_field, 0, yindex, 0); */
599       RaveField_setValue(nv_field, 0, yindex, 0);
600       RaveField_setValue(HGHT_field, 0, yindex,centerOfLayer);
601       RaveField_setValue(UWND_field, 0, yindex,self->nodata_VP);
602       RaveField_setValue(VWND_field, 0, yindex,self->nodata_VP);
603     } else {
604       /* Add wgen needed: RaveField_setValue(ff_field, 0, yindex, vvel);
605       RaveField_setValue(ff_dev_field, 0, yindex, vstd);
606       RaveField_setValue(dd_field, 0, yindex, vdir);
607       RaveField_setValue(dbzh_field, 0, yindex, zmean);
608       RaveField_setValue(dbzh_dev_field, 0, yindex, zstd);
609       RaveField_setValue(nz_field, 0, yindex, nz); */
610       RaveField_setValue(nv_field, 0, yindex, nv);
611       RaveField_setValue(HGHT_field, 0, yindex,centerOfLayer);
612       RaveField_setValue(UWND_field, 0, yindex,u_wnd_comp);
613       RaveField_setValue(VWND_field, 0, yindex,v_wnd_comp);
614     }
615
616     RAVE_FREE(A);
617     RAVE_FREE(b);
618     RAVE_FREE(v);
619     RAVE_FREE(z);
620     RAVE_FREE(az);
621
622     yindex++; /* Next interval*/
623   }
624
625   /* Set values used for /dataset1/what attributes */  
626   starttime = starttimeFirst;
627   endtime = endtimeFirst;
628   startdate = startdateFirst;
629   enddate = enddateFirst;
630
631   
632   result = RAVE_OBJECT_NEW(&VerticalProfile_TYPE);
633   if (result != NULL) {
634     if (!VerticalProfile_setUWND(result, UWND_field) ||
635         !VerticalProfile_setVWND(result, VWND_field) ||
636         !VerticalProfile_setNV(result, nv_field) ||
637         !VerticalProfile_setHGHT(result, HGHT_field)) {
638         /* Add when needed: !VerticalProfile_setFF(result, ff_field) ||
639         !VerticalProfile_setFFDev(result, ff_dev_field) ||
640         !VerticalProfile_setDD(result, dd_field) ||
641         !VerticalProfile_setDBZ(result, dbzh_field) ||
642         !VerticalProfile_setDBZDev(result, dbzh_dev_field)) */
643       RAVE_ERROR0("Failed to set vertical profile fields");
644       RAVE_OBJECT_RELEASE(result);
645     }
646   }
647
648   VerticalProfile_setLongitude(result, PolarVolume_getLongitude(inobj));
649   VerticalProfile_setLatitude(result, PolarVolume_getLatitude(inobj));
650   VerticalProfile_setHeight(result, PolarVolume_getHeight(inobj));
651   VerticalProfile_setSource(result, PolarVolume_getSource(inobj));
652   VerticalProfile_setInterval(result, self->dz);
653   VerticalProfile_setLevels(result, ysize);
654   VerticalProfile_setMinheight(result, 0);
655   VerticalProfile_setMaxheight(result, self->hmax);
656   VerticalProfile_setDate(result, PolarVolume_getDate(inobj));
657   VerticalProfile_setTime(result, PolarVolume_getTime(inobj));
658   
659   /* Set the times and product, starttime is the starttime for the lowest elev
660      endtime is the endtime for the highest elev. */
661   VerticalProfile_setStartTime(result, starttime);
662   VerticalProfile_setEndTime(result, endtime);
663   VerticalProfile_setStartDate(result, startdate);
664   VerticalProfile_setEndDate(result, enddate);
665   VerticalProfile_setProduct(result, product);
666    
667   /* Need to filter the attribute data with respect to selected min elangle
668      Below attributes satisfy reguirements from E-profile, add other when needed */
669   
670   WrwpInternal_findAndAddAttribute(result, inobj, "how/highprf", self->emin);
671   WrwpInternal_findAndAddAttribute(result, inobj, "how/lowprf", self->emin);
672   /*WrwpInternal_findAndAddAttribute(result, inobj, "how/pulsewidth", self->emin); */
673   WrwpInternal_findAndAddAttribute(result, inobj, "how/wavelength", self->emin);
674  /* WrwpInternal_findAndAddAttribute(result, inobj, "how/RXbandwidth", self->emin);
675   WrwpInternal_findAndAddAttribute(result, inobj, "how/RXlossH", self->emin);
676   WrwpInternal_findAndAddAttribute(result, inobj, "how/TXlossH", self->emin);
677   WrwpInternal_findAndAddAttribute(result, inobj, "how/antgainH", self->emin);
678   WrwpInternal_findAndAddAttribute(result, inobj, "how/azmethod", self->emin);
679   WrwpInternal_findAndAddAttribute(result, inobj, "how/binmethod", self->emin);
680   WrwpInternal_findAndAddAttribute(result, inobj, "how/malfunc", self->emin);
681   WrwpInternal_findAndAddAttribute(result, inobj, "how/nomTXpower", self->emin);
682   WrwpInternal_findAndAddAttribute(result, inobj, "how/radar_msg", self->emin);
683   WrwpInternal_findAndAddAttribute(result, inobj, "how/radconstH", self->emin);
684   WrwpInternal_findAndAddAttribute(result, inobj, "how/radomeloss", self->emin);
685   WrwpInternal_findAndAddAttribute(result, inobj, "how/rpm", self->emin);
686   WrwpInternal_findAndAddAttribute(result, inobj, "how/software", self->emin);
687   WrwpInternal_findAndAddAttribute(result, inobj, "how/system", self->emin);*/
688
689   WrwpInternal_addIntAttribute(result, "how/minrange", Wrwp_getDMIN(self));
690   WrwpInternal_addIntAttribute(result, "how/maxrange", Wrwp_getDMAX(self));
691   
692   
693   /* Put in selected attributes in selected fields, currently a nodata, gain and
694      offset are defined for UWND and VWND and are put under what/nodata. */
695   WrwpInternal_addDoubleAttr2Field(UWND_field, "what/nodata", self->nodata_VP);
696   WrwpInternal_addDoubleAttr2Field(UWND_field, "what/gain", self->gain_VP);
697   WrwpInternal_addDoubleAttr2Field(UWND_field, "what/offset", self->offset_VP);
698   WrwpInternal_addDoubleAttr2Field(UWND_field, "what/undetect", self->undetect_VP);
699   
700   VerticalProfile_addField(result, UWND_field);
701   
702   WrwpInternal_addDoubleAttr2Field(VWND_field, "what/nodata", self->nodata_VP);
703   WrwpInternal_addDoubleAttr2Field(VWND_field, "what/gain", self->gain_VP);
704   WrwpInternal_addDoubleAttr2Field(VWND_field, "what/offset", self->offset_VP);
705   WrwpInternal_addDoubleAttr2Field(VWND_field, "what/undetect", self->undetect_VP);
706   
707   VerticalProfile_addField(result, VWND_field);
708   
709   WrwpInternal_addDoubleAttr2Field(HGHT_field, "what/nodata", self->nodata_VP);
710   WrwpInternal_addDoubleAttr2Field(HGHT_field, "what/gain", self->gain_VP);
711   WrwpInternal_addDoubleAttr2Field(HGHT_field, "what/offset", self->offset_VP);
712   WrwpInternal_addDoubleAttr2Field(HGHT_field, "what/undetect", self->undetect_VP);
713   
714   VerticalProfile_addField(result, HGHT_field);
715   
716   WrwpInternal_addDoubleAttr2Field(nv_field, "what/nodata", self->nodata_VP);
717   WrwpInternal_addDoubleAttr2Field(nv_field, "what/gain", self->gain_VP);
718   WrwpInternal_addDoubleAttr2Field(nv_field, "what/offset", self->offset_VP);
719   WrwpInternal_addDoubleAttr2Field(nv_field, "what/undetect", self->undetect_VP);
720   
721   VerticalProfile_addField(result, nv_field);
722
723 done:
724   RAVE_OBJECT_RELEASE(polnav);
725   /* 
726   RAVE_OBJECT_RELEASE(ff_field);
727   RAVE_OBJECT_RELEASE(ff_dev_field);
728   RAVE_OBJECT_RELEASE(dd_field);
729   RAVE_OBJECT_RELEASE(dbzh_field);
730   RAVE_OBJECT_RELEASE(dbzh_dev_field);
731   RAVE_OBJECT_RELEASE(nz_field); */
732   RAVE_OBJECT_RELEASE(nv_field);
733   RAVE_OBJECT_RELEASE(HGHT_field);
734   RAVE_OBJECT_RELEASE(UWND_field);
735   RAVE_OBJECT_RELEASE(UWND_field);
736
737   return result;
738 }
739
740 /*@} End of Interface functions */
741
742 RaveCoreObjectType Wrwp_TYPE = {
743     "Wrwp",
744     sizeof(Wrwp_t),
745     Wrwp_constructor,
746     Wrwp_destructor
747 };
748