not-so-dirty implementation (actually tolerance) of unstructured grids

volgrid6d/grid_class.F90

@@ -773,6 +773,14 @@ SUBROUTINE griddim_import_gribapi(this, gaid)
 #endif
 CALL grib_get(gaid,'GRIBEditionNumber',EditionNumber)
 
+IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
+  CALL grib_get(gaid,'numberOfDataPoints', this%dim%nx)
+  this%dim%ny = 1
+  this%grid%grid%component_flag = 0
+  CALL griddim_import_ellipsoid(this, gaid)
+  RETURN
+ENDIF
+
 ! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
 CALL grib_get(gaid, 'Ni', this%dim%nx)
 CALL grib_get(gaid, 'Nj', this%dim%ny)
@@ -1150,6 +1158,14 @@ SUBROUTINE griddim_export_gribapi(this, gaid)
  "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
 #endif
 
+IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
+  CALL grib_set(gaid,'numberOfDataPoints', this%dim%nx)
+! reenable when possible
+!  CALL griddim_export_ellipsoid(this, gaid)
+  RETURN
+ENDIF
+
+
 ! Edition dependent setup
 IF (EditionNumber == 1) THEN
   ratio = 1.d3
@@ -1471,19 +1487,20 @@ SUBROUTINE griddim_export_ellipsoid(this, gaid)
 TYPE(griddim_def),INTENT(in) :: this
 INTEGER,INTENT(in) :: gaid
 
-INTEGER :: ellips_type
+INTEGER :: ellips_type, ierr
 DOUBLE PRECISION :: r1, r2, f
 
 CALL get_val(this, ellips_smaj_axis=r1, ellips_flatt=f, ellips_type=ellips_type)
 
 IF (EditionNumber == 2) THEN
 
-  CALL grib_set_missing(gaid, 'scaleFactorOfRadiusOfSphericalEarth')
-  CALL grib_set_missing(gaid, 'scaledValueOfRadiusOfSphericalEarth')
-  CALL grib_set_missing(gaid, 'scaleFactorOfEarthMajorAxis')
-  CALL grib_set_missing(gaid, 'scaledValueOfEarthMajorAxis')
-  CALL grib_set_missing(gaid, 'scaleFactorOfEarthMinorAxis')
-  CALL grib_set_missing(gaid, 'scaledValueOfEarthMinorAxis')
+! why does it produce a message even with ierr?
+  CALL grib_set_missing(gaid, 'scaleFactorOfRadiusOfSphericalEarth', ierr)
+  CALL grib_set_missing(gaid, 'scaledValueOfRadiusOfSphericalEarth', ierr)
+  CALL grib_set_missing(gaid, 'scaleFactorOfEarthMajorAxis', ierr)
+  CALL grib_set_missing(gaid, 'scaledValueOfEarthMajorAxis', ierr)
+  CALL grib_set_missing(gaid, 'scaleFactorOfEarthMinorAxis', ierr)
+  CALL grib_set_missing(gaid, 'scaledValueOfEarthMinorAxis', ierr)
 
   SELECT CASE(ellips_type)
   CASE(ellips_grs80) ! iag-grs80

volgrid6d/grid_id_class.F90

@@ -817,21 +817,21 @@ SUBROUTINE grid_id_decode_data_gribapi(gaid, field)
  iScansNegatively, jScansPositively, jPointsAreConsecutive
 INTEGER :: numberOfValues,numberOfPoints
 REAL :: vector(SIZE(field))
-INTEGER :: x1, x2, xs, y1, y2, ys, ord(2)
+INTEGER :: x1, x2, xs, y1, y2, ys, ord(2), ierr
 
 
 call grib_get(gaid,'GRIBEditionNumber',EditionNumber)
 
 if (EditionNumber == 2) then
 
-  call grib_get(gaid,'alternativeRowScanning',alternativeRowScanning)
-  if (alternativeRowScanning /= 0)then
-    call l4f_log(L4F_ERROR, "grib_api alternativeRowScanning not supported: " &
+  CALL grib_get(gaid,'alternativeRowScanning',alternativeRowScanning,ierr)
+  IF (ierr == GRIB_SUCCESS .AND. alternativeRowScanning /= 0) THEN
+    CALL l4f_log(L4F_ERROR, "grib_api alternativeRowScanning not supported: " &
      //t2c(alternativeRowScanning))
-    call raise_error()
+    CALL raise_error()
     field(:,:) = rmiss
     RETURN
-  end if
+  ENDIF
 
 else if (EditionNumber /= 1) then
 
@@ -843,9 +843,12 @@ SUBROUTINE grid_id_decode_data_gribapi(gaid, field)
 
 end if
 
-call grib_get(gaid,'iScansNegatively',iScansNegatively)
-call grib_get(gaid,'jScansPositively',jScansPositively)
-call grib_get(gaid,'jPointsAreConsecutive',jPointsAreConsecutive)
+CALL grib_get(gaid,'iScansNegatively',iScansNegatively,ierr)
+IF (ierr /= GRIB_SUCCESS) iScansNegatively=0
+CALL grib_get(gaid,'jScansPositively',jScansPositively,ierr)
+IF (ierr /= GRIB_SUCCESS) jScansPositively=1
+CALL grib_get(gaid,'jPointsAreConsecutive',jPointsAreConsecutive,ierr)
+IF (ierr /= GRIB_SUCCESS) jPointsAreConsecutive=0
 
 call grib_get(gaid,'numberOfPoints',numberOfPoints)
 call grib_get(gaid,'numberOfValues',numberOfValues)
@@ -924,19 +927,19 @@ SUBROUTINE grid_id_encode_data_gribapi(gaid, field)
 INTEGER :: EditionNumber
 INTEGER :: alternativeRowScanning, iScansNegatively, &
  jScansPositively, jPointsAreConsecutive
-INTEGER :: x1, x2, xs, y1, y2, ys
+INTEGER :: x1, x2, xs, y1, y2, ys, ierr
 
 call grib_get(gaid,'GRIBEditionNumber',EditionNumber)
 
 if (EditionNumber == 2) then
 
-  call grib_get(gaid,'alternativeRowScanning',alternativeRowScanning)
-  if (alternativeRowScanning /= 0)then
-    call l4f_log(L4F_ERROR, "grib_api alternativeRowScanning not supported: " &
+  CALL grib_get(gaid,'alternativeRowScanning',alternativeRowScanning,ierr)
+  IF (ierr == GRIB_SUCCESS .AND. alternativeRowScanning /= 0)THEN
+    CALL l4f_log(L4F_ERROR, "grib_api alternativeRowScanning not supported: " &
      //TRIM(to_char(alternativeRowScanning)))
-    call raise_error()
+    CALL raise_error()
     RETURN
-  end if
+  ENDIF
 
 else if( EditionNumber /= 1) then
 
@@ -947,16 +950,19 @@ SUBROUTINE grid_id_encode_data_gribapi(gaid, field)
 
 end if
 
-call grib_get(gaid,'iScansNegatively',iScansNegatively)
-call grib_get(gaid,'jScansPositively',jScansPositively)
-call grib_get(gaid,'jPointsAreConsecutive',jPointsAreConsecutive)
+CALL grib_get(gaid,'iScansNegatively',iScansNegatively,ierr)
+IF (ierr /= GRIB_SUCCESS) iScansNegatively=0
+CALL grib_get(gaid,'jScansPositively',jScansPositively,ierr)
+IF (ierr /= GRIB_SUCCESS) jScansPositively=1
+CALL grib_get(gaid,'jPointsAreConsecutive',jPointsAreConsecutive,ierr)
+IF (ierr /= GRIB_SUCCESS) jPointsAreConsecutive=0
 
 ! queste sono gia` fatte in export_gridinfo, si potrebbero evitare?!
 #ifdef DEBUG
 CALL l4f_log(L4F_DEBUG, 'grib_api, Ni,Nj:'//t2c(SIZE(field,1))//','//t2c(SIZE(field,2)))
 #endif
-call grib_set(gaid,'Ni',SIZE(field,1))
-call grib_set(gaid,'Nj',SIZE(field,2))
+!call grib_set(gaid,'Ni',SIZE(field,1))
+!call grib_set(gaid,'Nj',SIZE(field,2))
 
 ! Transfer data field changing scanning mode from 64
 IF (iScansNegatively  == 0) THEN

volgrid6d/grid_rect_class.F90

@@ -273,11 +273,15 @@ SUBROUTINE grid_rect_coordinates(this, x, y)
 #endif
 
 CALL grid_rect_steps(this, nx, ny, dx, dy)
-
-x(:,:) = RESHAPE((/ ((this%xmin+(dx*DBLE(i)), i=0,nx-1), j=0,ny-1) /),&
- (/nx,ny/))
-y(:,:) = RESHAPE((/ ((this%ymin+(dy*DBLE(j)), i=0,nx-1), j=0,ny-1) /),&
- (/nx,ny/))
+IF (c_e(dx) .AND. c_e(dy)) THEN
+  x(:,:) = RESHAPE((/ ((this%xmin+(dx*DBLE(i)), i=0,nx-1), j=0,ny-1) /),&
+   (/nx,ny/))
+  y(:,:) = RESHAPE((/ ((this%ymin+(dy*DBLE(j)), i=0,nx-1), j=0,ny-1) /),&
+   (/nx,ny/))
+ELSE
+  x(:,:) = dmiss
+  y(:,:) = dmiss
+ENDIF
 
 END SUBROUTINE grid_rect_coordinates
 

volgrid6d/volgrid6d_class.F90

@@ -518,7 +518,7 @@ END SUBROUTINE volgrid_set_vol_2d
 !! This method works both with volumes having allocated and
 !! non-allocated this%voldati array, and it updates the requested
 !! slice.  In case \a this%voldati is already allocated, this is a
-!! no-operation while in the other case this method encodes the filed
+!! no-operation while in the other case this method encodes the field
 !! provided into the grid_id object on file or in memory. Since this
 !! method may be called many times by a program, it is optimized for
 !! speed and it does not make any check about the matching size of the