libsim  Versione 7.1.9
grid_class.F90
1 ! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2 ! authors:
3 ! Davide Cesari <dcesari@arpa.emr.it>
4 ! Paolo Patruno <ppatruno@arpa.emr.it>
5 
6 ! This program is free software; you can redistribute it and/or
7 ! modify it under the terms of the GNU General Public License as
8 ! published by the Free Software Foundation; either version 2 of
9 ! the License, or (at your option) any later version.
10 
11 ! This program is distributed in the hope that it will be useful,
12 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
13 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 ! GNU General Public License for more details.
15 
16 ! You should have received a copy of the GNU General Public License
17 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
18 #include "config.h"
19 
49 MODULE grid_class
50 use geo_proj_class
52 use grid_rect_class
53 use grid_id_class
54 use err_handling
57 use log4fortran
58 implicit none
59 
60 
61 character (len=255),parameter:: subcategory="grid_class"
62 
63 
69 type grid_def
70  private
71  type(geo_proj) :: proj
72  type(grid_rect) :: grid
73  integer :: category = 0
74 end type grid_def
75 
76 
82 type griddim_def
83  type(grid_def) :: grid
84  type(grid_dim) :: dim
85  integer :: category = 0
86 end type griddim_def
87 
88 
92 INTERFACE OPERATOR (==)
93  MODULE PROCEDURE grid_eq, griddim_eq
94 END INTERFACE
95 
99 INTERFACE OPERATOR (/=)
100  MODULE PROCEDURE grid_ne, griddim_ne
101 END INTERFACE
102 
104 INTERFACE init
105  MODULE PROCEDURE griddim_init
106 END INTERFACE
107 
109 INTERFACE delete
110  MODULE PROCEDURE griddim_delete
111 END INTERFACE
112 
114 INTERFACE copy
115  MODULE PROCEDURE griddim_copy
116 END INTERFACE
117 
120 INTERFACE proj
121  MODULE PROCEDURE griddim_coord_proj!, griddim_proj
122 END INTERFACE
123 
126 INTERFACE unproj
127  MODULE PROCEDURE griddim_coord_unproj, griddim_unproj
128 END INTERFACE
129 
131 INTERFACE get_val
132  MODULE PROCEDURE griddim_get_val
133 END INTERFACE
134 
136 INTERFACE set_val
137  MODULE PROCEDURE griddim_set_val
138 END INTERFACE
139 
141 INTERFACE write_unit
142  MODULE PROCEDURE griddim_write_unit
143 END INTERFACE
144 
146 INTERFACE read_unit
147  MODULE PROCEDURE griddim_read_unit
148 END INTERFACE
149 
151 INTERFACE import
152  MODULE PROCEDURE griddim_import_grid_id
153 END INTERFACE
154 
156 INTERFACE export
157  MODULE PROCEDURE griddim_export_grid_id
158 END INTERFACE
159 
161 INTERFACE display
162  MODULE PROCEDURE griddim_display
163 END INTERFACE
164 
165 #define VOL7D_POLY_TYPE TYPE(grid_def)
166 #define VOL7D_POLY_TYPES _grid
167 #include "array_utilities_pre.F90"
168 #undef VOL7D_POLY_TYPE
169 #undef VOL7D_POLY_TYPES
170 
171 #define VOL7D_POLY_TYPE TYPE(griddim_def)
172 #define VOL7D_POLY_TYPES _griddim
173 #include "array_utilities_pre.F90"
174 #undef VOL7D_POLY_TYPE
175 #undef VOL7D_POLY_TYPES
176 
177 INTERFACE wind_unrot
178  MODULE PROCEDURE griddim_wind_unrot
179 END INTERFACE
180 
181 
182 PRIVATE
183 
184 PUBLIC proj, unproj, griddim_unproj, griddim_gen_coord, &
185  griddim_zoom_coord, griddim_zoom_projcoord, &
186  griddim_setsteps, griddim_def, grid_def, grid_dim
187 PUBLIC init, delete, copy
189 PUBLIC OPERATOR(==),OPERATOR(/=)
190 PUBLIC count_distinct, pack_distinct, count_and_pack_distinct, &
191  map_distinct, map_inv_distinct,index
192 PUBLIC wind_unrot, import, export
193 PUBLIC griddim_central_lon, griddim_set_central_lon
194 CONTAINS
195 
197 SUBROUTINE griddim_init(this, nx, ny, &
198  xmin, xmax, ymin, ymax, dx, dy, component_flag, &
199  proj_type, lov, zone, xoff, yoff, &
200  longitude_south_pole, latitude_south_pole, angle_rotation, &
201  longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
202  latin1, latin2, lad, projection_center_flag, &
203  ellips_smaj_axis, ellips_flatt, ellips_type, &
204  categoryappend)
205 TYPE(griddim_def),INTENT(inout) :: this
206 INTEGER,INTENT(in),OPTIONAL :: nx
207 INTEGER,INTENT(in),OPTIONAL :: ny
208 DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin
209 DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmax
210 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymin
211 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymax
212 DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx
213 DOUBLE PRECISION,INTENT(in),OPTIONAL :: dy
216 INTEGER,INTENT(in),OPTIONAL :: component_flag
217 CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
218 DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
219 INTEGER,INTENT(in),OPTIONAL :: zone
220 DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
221 DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
222 DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
223 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
224 DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
225 DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
226 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
227 DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
228 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
229 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
230 DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
231 INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
232 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
233 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
234 INTEGER,INTENT(in),OPTIONAL :: ellips_type
235 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
236 
237 CHARACTER(len=512) :: a_name
238 
239 IF (PRESENT(categoryappend)) THEN
240  CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
241  trim(categoryappend))
242 ELSE
243  CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
244 ENDIF
245 this%category=l4f_category_get(a_name)
246 
247 ! geographical projection
248 this%grid%proj = geo_proj_new( &
249  proj_type=proj_type, lov=lov, zone=zone, xoff=xoff, yoff=yoff, &
250  longitude_south_pole=longitude_south_pole, &
251  latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
252  longitude_stretch_pole=longitude_stretch_pole, &
253  latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
254  latin1=latin1, latin2=latin2, lad=lad, projection_center_flag=projection_center_flag, &
255  ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, ellips_type=ellips_type)
256 ! grid extension
257 this%grid%grid = grid_rect_new( &
258  xmin, xmax, ymin, ymax, dx, dy, component_flag)
259 ! grid size
260 this%dim = grid_dim_new(nx, ny)
261 
262 #ifdef DEBUG
263 call l4f_category_log(this%category,l4f_debug,"init gtype: "//this%grid%proj%proj_type )
264 #endif
265 
266 END SUBROUTINE griddim_init
267 
268 
270 SUBROUTINE griddim_delete(this)
271 TYPE(griddim_def),INTENT(inout) :: this
272 
273 CALL delete(this%dim)
274 CALL delete(this%grid%proj)
275 CALL delete(this%grid%grid)
276 
277 call l4f_category_delete(this%category)
278 
279 END SUBROUTINE griddim_delete
280 
281 
283 SUBROUTINE griddim_copy(this, that, categoryappend)
284 TYPE(griddim_def),INTENT(in) :: this
285 TYPE(griddim_def),INTENT(out) :: that
286 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
287 
288 CHARACTER(len=512) :: a_name
289 
290 CALL init(that)
291 
292 CALL copy(this%grid%proj, that%grid%proj)
293 CALL copy(this%grid%grid, that%grid%grid)
294 CALL copy(this%dim, that%dim)
295 
296 ! new category
297 IF (PRESENT(categoryappend)) THEN
298  CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
299  trim(categoryappend))
300 ELSE
301  CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
302 ENDIF
303 that%category=l4f_category_get(a_name)
304 
305 END SUBROUTINE griddim_copy
306 
307 
310 ELEMENTAL SUBROUTINE griddim_coord_proj(this, lon, lat, x, y)
311 TYPE(griddim_def),INTENT(in) :: this
313 DOUBLE PRECISION,INTENT(in) :: lon, lat
315 DOUBLE PRECISION,INTENT(out) :: x, y
316 
317 CALL proj(this%grid%proj, lon, lat, x, y)
318 
319 END SUBROUTINE griddim_coord_proj
320 
321 
324 ELEMENTAL SUBROUTINE griddim_coord_unproj(this, x, y, lon, lat)
325 TYPE(griddim_def),INTENT(in) :: this
327 DOUBLE PRECISION,INTENT(in) :: x, y
329 DOUBLE PRECISION,INTENT(out) :: lon, lat
330 
331 CALL unproj(this%grid%proj, x, y, lon, lat)
332 
333 END SUBROUTINE griddim_coord_unproj
334 
335 
336 ! Computes and sets the grid parameters required to compute
337 ! coordinates of grid points in the projected system.
338 ! probably meaningless
339 !SUBROUTINE griddim_proj(this)
340 !TYPE(griddim_def),INTENT(inout) :: this !< definition of projection and grid
341 !
342 !CALL proj(this, this%dim%lon(1,1), this%dim%lat(1,1), &
343 ! this%grid%grid%xmin, this%grid%grid%ymin)
344 !
345 !CALL proj(this, this%dim%lon(this%dim%nx,this%dim%ny), &
346 ! this%dim%lat(this%dim%nx,this%dim%ny), &
347 ! this%grid%grid%xmax, this%grid%grid%ymax)
348 !
349 !END SUBROUTINE griddim_proj
350 
358 SUBROUTINE griddim_unproj(this)
359 TYPE(griddim_def),INTENT(inout) :: this
360 
361 IF (.NOT.c_e(this%dim%nx) .OR. .NOT.c_e(this%dim%ny)) RETURN
362 CALL alloc(this%dim)
363 CALL griddim_unproj_internal(this)
364 
365 END SUBROUTINE griddim_unproj
366 
367 ! internal subroutine needed for allocating automatic arrays
368 SUBROUTINE griddim_unproj_internal(this)
369 TYPE(griddim_def),INTENT(inout) ::this ! definition of projection and grid
370 
371 DOUBLE PRECISION :: x(this%dim%nx,this%dim%ny), y(this%dim%nx,this%dim%ny)
372 
373 CALL grid_rect_coordinates(this%grid%grid, x, y)
374 CALL unproj(this, x, y, this%dim%lon, this%dim%lat)
375 
376 END SUBROUTINE griddim_unproj_internal
377 
378 
380 SUBROUTINE griddim_get_val(this, nx, ny, &
381  xmin, xmax, ymin, ymax, dx, dy, component_flag, &
382  proj, proj_type, lov, zone, xoff, yoff, &
383  longitude_south_pole, latitude_south_pole, angle_rotation, &
384  longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
385  latin1, latin2, lad, projection_center_flag, &
386  ellips_smaj_axis, ellips_flatt, ellips_type)
387 TYPE(griddim_def),INTENT(in) :: this
388 INTEGER,INTENT(out),OPTIONAL :: nx
389 INTEGER,INTENT(out),OPTIONAL :: ny
391 DOUBLE PRECISION,INTENT(out),OPTIONAL :: xmin, xmax, ymin, ymax
393 DOUBLE PRECISION,INTENT(out),OPTIONAL :: dx, dy
396 INTEGER,INTENT(out),OPTIONAL :: component_flag
397 TYPE(geo_proj),INTENT(out),OPTIONAL :: proj
398 CHARACTER(len=*),INTENT(out),OPTIONAL :: proj_type
399 DOUBLE PRECISION,INTENT(out),OPTIONAL :: lov
400 INTEGER,INTENT(out),OPTIONAL :: zone
401 DOUBLE PRECISION,INTENT(out),OPTIONAL :: xoff
402 DOUBLE PRECISION,INTENT(out),OPTIONAL :: yoff
403 DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_south_pole
404 DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_south_pole
405 DOUBLE PRECISION,INTENT(out),OPTIONAL :: angle_rotation
406 DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_stretch_pole
407 DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_stretch_pole
408 DOUBLE PRECISION,INTENT(out),OPTIONAL :: stretch_factor
409 DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin1
410 DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin2
411 DOUBLE PRECISION,INTENT(out),OPTIONAL :: lad
412 INTEGER,INTENT(out),OPTIONAL :: projection_center_flag
413 DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_smaj_axis
414 DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_flatt
415 INTEGER,INTENT(out),OPTIONAL :: ellips_type
416 
417 IF (PRESENT(nx)) nx = this%dim%nx
418 IF (PRESENT(ny)) ny = this%dim%ny
419 
420 IF (PRESENT(proj)) proj = this%grid%proj
421 
422 CALL get_val(this%grid%proj, proj_type=proj_type, lov=lov, zone=zone, &
423  xoff=xoff, yoff=yoff, &
424  longitude_south_pole=longitude_south_pole, &
425  latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
426  longitude_stretch_pole=longitude_stretch_pole, &
427  latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
428  latin1=latin1, latin2=latin2, lad=lad, &
429  projection_center_flag=projection_center_flag, &
430  ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
431  ellips_type=ellips_type)
432 
433 CALL get_val(this%grid%grid, &
434  xmin, xmax, ymin, ymax, dx, dy, component_flag)
435 
436 END SUBROUTINE griddim_get_val
437 
438 
440 SUBROUTINE griddim_set_val(this, nx, ny, &
441  xmin, xmax, ymin, ymax, dx, dy, component_flag, &
442  proj_type, lov, zone, xoff, yoff, &
443  longitude_south_pole, latitude_south_pole, angle_rotation, &
444  longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
445  latin1, latin2, lad, projection_center_flag, &
446  ellips_smaj_axis, ellips_flatt, ellips_type)
447 TYPE(griddim_def),INTENT(inout) :: this
448 INTEGER,INTENT(in),OPTIONAL :: nx
449 INTEGER,INTENT(in),OPTIONAL :: ny
451 DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin, xmax, ymin, ymax
453 DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx, dy
456 INTEGER,INTENT(in),OPTIONAL :: component_flag
457 CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
458 DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
459 INTEGER,INTENT(in),OPTIONAL :: zone
460 DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
461 DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
462 DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
463 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
464 DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
465 DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
466 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
467 DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
468 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
469 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
470 DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
471 INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
472 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
473 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
474 INTEGER,INTENT(in),OPTIONAL :: ellips_type
475 
476 IF (PRESENT(nx)) this%dim%nx = nx
477 IF (PRESENT(ny)) this%dim%ny = ny
478 
479 CALL set_val(this%grid%proj, proj_type=proj_type, lov=lov, zone=zone, &
480  xoff=xoff, yoff=yoff, longitude_south_pole=longitude_south_pole, &
481  latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
482  longitude_stretch_pole=longitude_stretch_pole, &
483  latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
484  latin1=latin1, latin2=latin2, lad=lad, &
485  projection_center_flag=projection_center_flag, &
486  ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
487  ellips_type=ellips_type)
488 
489 CALL set_val(this%grid%grid, &
490  xmin, xmax, ymin, ymax, dx, dy, component_flag)
491 
492 END SUBROUTINE griddim_set_val
493 
494 
499 SUBROUTINE griddim_read_unit(this, unit)
500 TYPE(griddim_def),INTENT(out) :: this
501 INTEGER, INTENT(in) :: unit
502 
503 
504 CALL read_unit(this%dim, unit)
505 CALL read_unit(this%grid%proj, unit)
506 CALL read_unit(this%grid%grid, unit)
507 
508 END SUBROUTINE griddim_read_unit
509 
510 
515 SUBROUTINE griddim_write_unit(this, unit)
516 TYPE(griddim_def),INTENT(in) :: this
517 INTEGER, INTENT(in) :: unit
518 
519 
520 CALL write_unit(this%dim, unit)
521 CALL write_unit(this%grid%proj, unit)
522 CALL write_unit(this%grid%grid, unit)
523 
524 END SUBROUTINE griddim_write_unit
525 
526 
530 FUNCTION griddim_central_lon(this) RESULT(lon)
531 TYPE(griddim_def),INTENT(inout) :: this
532 
533 DOUBLE PRECISION :: lon
534 
535 CALL griddim_pistola_central_lon(this, lon)
536 
537 END FUNCTION griddim_central_lon
538 
539 
543 SUBROUTINE griddim_set_central_lon(this, lonref)
544 TYPE(griddim_def),INTENT(inout) :: this
545 DOUBLE PRECISION,INTENT(in) :: lonref
546 
547 DOUBLE PRECISION :: lon
548 
549 CALL griddim_pistola_central_lon(this, lon, lonref)
550 
551 END SUBROUTINE griddim_set_central_lon
552 
553 
554 ! internal subroutine for performing tasks common to the prevous two
555 SUBROUTINE griddim_pistola_central_lon(this, lon, lonref)
556 TYPE(griddim_def),INTENT(inout) :: this ! grid descriptor
557 DOUBLE PRECISION,INTENT(inout) :: lon ! central longitude
558 DOUBLE PRECISION,INTENT(in),OPTIONAL :: lonref ! reference longitude
559 
560 INTEGER :: unit
561 DOUBLE PRECISION :: lonsp, latsp, londelta, lov, lonrot
562 CHARACTER(len=80) :: ptype
563 
564 lon = dmiss
565 CALL get_val(this%grid%proj, unit=unit)
566 IF (unit == geo_proj_unit_meter) THEN ! it is a plane projection
567  CALL get_val(this%grid%proj, lov=lon)
568  IF (PRESENT(lonref)) THEN
569  CALL long_reset_to_cart_closest(lov, lonref)
570  CALL set_val(this%grid%proj, lov=lon)
571  ENDIF
572 
573 ELSE IF (unit == geo_proj_unit_degree) THEN ! it is a spheric projection
574  CALL get_val(this%grid%proj, proj_type=ptype, &
575  longitude_south_pole=lonsp, latitude_south_pole=latsp)
576  SELECT CASE(ptype)
577  CASE('rotated_ll','stretched_rotated_ll') ! use origin of rotated system
578  IF (latsp < 0.0d0) THEN
579  lon = lonsp
580  IF (PRESENT(lonref)) THEN
581  CALL long_reset_to_cart_closest(lon, lonref)
582  CALL set_val(this%grid%proj, longitude_south_pole=lonref)
583 ! now reset rotated coordinates around zero
584  IF (c_e(this%grid%grid%xmin) .AND. c_e(this%grid%grid%xmax)) THEN
585  lonrot = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
586  ENDIF
587  londelta = lonrot
588  CALL long_reset_to_cart_closest(londelta, 0.0d0)
589  londelta = londelta - lonrot
590  this%grid%grid%xmin = this%grid%grid%xmin + londelta
591  this%grid%grid%xmax = this%grid%grid%xmax + londelta
592  ENDIF
593  ELSE ! this part to be checked
594  lon = modulo(lonsp + 180.0d0, 360.0d0)
595 ! IF (PRESENT(lonref)) THEN
596 ! CALL long_reset_to_cart_closest(lov, lonref)
597 ! CALL set_val(this%grid%proj, longitude_south_pole=lonref)
598 ! ENDIF
599  ENDIF
600  CASE default ! use real grid limits
601  IF (c_e(this%grid%grid%xmin) .AND. c_e(this%grid%grid%xmax)) THEN
602  lon = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
603  ENDIF
604  IF (PRESENT(lonref)) THEN
605  londelta = lon
606  CALL long_reset_to_cart_closest(londelta, lonref)
607  londelta = londelta - lon
608  this%grid%grid%xmin = this%grid%grid%xmin + londelta
609  this%grid%grid%xmax = this%grid%grid%xmax + londelta
610  ENDIF
611  END SELECT
612 ENDIF
613 
614 END SUBROUTINE griddim_pistola_central_lon
615 
616 
620 SUBROUTINE griddim_gen_coord(this, x, y)
621 TYPE(griddim_def),INTENT(in) :: this
622 DOUBLE PRECISION,INTENT(out) :: x(:,:)
623 DOUBLE PRECISION,INTENT(out) :: y(:,:)
624 
625 
626 CALL grid_rect_coordinates(this%grid%grid, x, y)
627 
628 END SUBROUTINE griddim_gen_coord
629 
630 
632 SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
633 TYPE(griddim_def), INTENT(in) :: this
634 INTEGER,INTENT(in) :: nx
635 INTEGER,INTENT(in) :: ny
636 DOUBLE PRECISION,INTENT(out) :: dx
637 DOUBLE PRECISION,INTENT(out) :: dy
638 
639 CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
640 
641 END SUBROUTINE griddim_steps
642 
643 
645 SUBROUTINE griddim_setsteps(this)
646 TYPE(griddim_def), INTENT(inout) :: this
647 
648 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
649 
650 END SUBROUTINE griddim_setsteps
651 
652 
653 ! TODO
654 ! bisogna sviluppare gli altri operatori
655 ELEMENTAL FUNCTION grid_eq(this, that) RESULT(res)
656 TYPE(grid_def),INTENT(IN) :: this, that
657 LOGICAL :: res
658 
659 res = this%proj == that%proj .AND. &
660  this%grid == that%grid
661 
662 END FUNCTION grid_eq
663 
664 
665 ELEMENTAL FUNCTION griddim_eq(this, that) RESULT(res)
666 TYPE(griddim_def),INTENT(IN) :: this, that
667 LOGICAL :: res
668 
669 res = this%grid == that%grid .AND. &
670  this%dim == that%dim
671 
672 END FUNCTION griddim_eq
673 
674 
675 ELEMENTAL FUNCTION grid_ne(this, that) RESULT(res)
676 TYPE(grid_def),INTENT(IN) :: this, that
677 LOGICAL :: res
678 
679 res = .NOT.(this == that)
680 
681 END FUNCTION grid_ne
682 
683 
684 ELEMENTAL FUNCTION griddim_ne(this, that) RESULT(res)
685 TYPE(griddim_def),INTENT(IN) :: this, that
686 LOGICAL :: res
687 
688 res = .NOT.(this == that)
689 
690 END FUNCTION griddim_ne
691 
692 
698 SUBROUTINE griddim_import_grid_id(this, ingrid_id)
699 #ifdef HAVE_LIBGDAL
700 USE gdal
701 #endif
702 TYPE(griddim_def),INTENT(inout) :: this
703 TYPE(grid_id),INTENT(in) :: ingrid_id
704 
705 #ifdef HAVE_LIBGRIBAPI
706 INTEGER :: gaid
707 #endif
708 #ifdef HAVE_LIBGDAL
709 TYPE(gdalrasterbandh) :: gdalid
710 #endif
711 CALL init(this)
712 
713 #ifdef HAVE_LIBGRIBAPI
714 gaid = grid_id_get_gaid(ingrid_id)
715 IF (c_e(gaid)) CALL griddim_import_gribapi(this, gaid)
716 #endif
717 #ifdef HAVE_LIBGDAL
718 gdalid = grid_id_get_gdalid(ingrid_id)
719 IF (gdalassociated(gdalid)) CALL griddim_import_gdal(this, gdalid, &
720  grid_id_get_gdal_options(ingrid_id))
721 #endif
722 
723 END SUBROUTINE griddim_import_grid_id
724 
725 
730 SUBROUTINE griddim_export_grid_id(this, outgrid_id)
731 #ifdef HAVE_LIBGDAL
732 USE gdal
733 #endif
734 TYPE(griddim_def),INTENT(in) :: this
735 TYPE(grid_id),INTENT(inout) :: outgrid_id
736 
737 #ifdef HAVE_LIBGRIBAPI
738 INTEGER :: gaid
739 #endif
740 #ifdef HAVE_LIBGDAL
741 TYPE(gdalrasterbandh) :: gdalid
742 #endif
743 
744 #ifdef HAVE_LIBGRIBAPI
745 gaid = grid_id_get_gaid(outgrid_id)
746 IF (c_e(gaid)) CALL griddim_export_gribapi(this, gaid)
747 #endif
748 #ifdef HAVE_LIBGDAL
749 gdalid = grid_id_get_gdalid(outgrid_id)
750 !IF (gdalassociated(gdalid)
751 ! export for gdal not implemented, log?
752 #endif
753 
754 END SUBROUTINE griddim_export_grid_id
755 
756 
757 #ifdef HAVE_LIBGRIBAPI
758 ! grib_api driver
759 SUBROUTINE griddim_import_gribapi(this, gaid)
760 USE grib_api
761 TYPE(griddim_def),INTENT(inout) :: this ! griddim object
762 INTEGER, INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
763 
764 DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, x1, y1, orient
765 INTEGER :: EditionNumber, iScansNegatively, jScansPositively, zone, datum, &
766  reflon, ierr
767 
768 ! Generic keys
769 CALL grib_get(gaid, 'typeOfGrid', this%grid%proj%proj_type)
770 #ifdef DEBUG
771 call l4f_category_log(this%category,l4f_debug, &
772  "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type))
773 #endif
774 CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
775 
776 IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
777  CALL grib_get(gaid,'numberOfDataPoints', this%dim%nx)
778  this%dim%ny = 1
779  this%grid%grid%component_flag = 0
780  CALL griddim_import_ellipsoid(this, gaid)
781  RETURN
782 ENDIF
783 
784 ! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
785 CALL grib_get(gaid, 'Ni', this%dim%nx)
786 CALL grib_get(gaid, 'Nj', this%dim%ny)
787 CALL griddim_import_ellipsoid(this, gaid) ! placed here, not valid for utm datum /= 1
788 
789 CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
790 CALL grib_get(gaid,'jScansPositively',jscanspositively)
791 CALL grib_get(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
792 
793 ! Keys for rotated grids (checked through missing values)
794 CALL grib_get_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
795  this%grid%proj%rotated%longitude_south_pole)
796 CALL grib_get_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
797  this%grid%proj%rotated%latitude_south_pole)
798 CALL grib_get_dmiss(gaid,'angleOfRotationInDegrees', &
799  this%grid%proj%rotated%angle_rotation)
800 
801 ! Keys for stretched grids (checked through missing values)
802 ! units must be verified, still experimental in grib_api
803 ! # TODO: Is it a float? Is it signed?
804 IF (editionnumber == 1) THEN
805  CALL grib_get_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
806  this%grid%proj%stretched%longitude_stretch_pole)
807  CALL grib_get_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
808  this%grid%proj%stretched%latitude_stretch_pole)
809  CALL grib_get_dmiss(gaid,'stretchingFactor', &
810  this%grid%proj%stretched%stretch_factor)
811 ELSE IF (editionnumber == 2) THEN
812  CALL grib_get_dmiss(gaid,'longitudeOfThePoleOfStretching', &
813  this%grid%proj%stretched%longitude_stretch_pole)
814  CALL grib_get_dmiss(gaid,'latitudeOfThePoleOfStretching', &
815  this%grid%proj%stretched%latitude_stretch_pole)
816  CALL grib_get_dmiss(gaid,'stretchingFactor', &
817  this%grid%proj%stretched%stretch_factor)
818  IF (c_e(this%grid%proj%stretched%stretch_factor)) &
819  this%grid%proj%stretched%stretch_factor = this%grid%proj%stretched%stretch_factor*1.0d-6
820 ENDIF
821 
822 ! Projection-dependent keys
823 SELECT CASE (this%grid%proj%proj_type)
824 
825 ! Keys for spherical coordinate systems
826 CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
827 
828  CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
829  CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
830  CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
831  CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
832 
833 ! longitudes are sometimes wrongly coded even in grib2 and even by the
834 ! Metoffice!
835 ! longitudeOfFirstGridPointInDegrees = 354.911;
836 ! longitudeOfLastGridPointInDegrees = 363.311;
837  CALL long_reset_0_360(lofirst)
838  CALL long_reset_0_360(lolast)
839 
840  IF (iscansnegatively == 0) THEN
841  this%grid%grid%xmin = lofirst
842  this%grid%grid%xmax = lolast
843  ELSE
844  this%grid%grid%xmax = lofirst
845  this%grid%grid%xmin = lolast
846  ENDIF
847  IF (jscanspositively == 0) THEN
848  this%grid%grid%ymax = lafirst
849  this%grid%grid%ymin = lalast
850  ELSE
851  this%grid%grid%ymin = lafirst
852  this%grid%grid%ymax = lalast
853  ENDIF
854 
855 ! reset longitudes in order to have a Cartesian plane
856  IF (this%grid%grid%xmax < this%grid%grid%xmin) &
857  this%grid%grid%xmin = this%grid%grid%xmin - 360.d0
858 
859 ! compute dx and dy (should we get them from grib?)
860  CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
861 
862 ! Keys for polar projections
863 CASE ('polar_stereographic', 'lambert', 'albers')
864 
865  CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
866  CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
867 ! latin1/latin2 may be missing (e.g. stereographic)
868  CALL grib_get_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
869  CALL grib_get_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
870 #ifdef DEBUG
871  CALL l4f_category_log(this%category,l4f_debug, &
872  "griddim_import_gribapi, latin1/2 "// &
873  trim(to_char(this%grid%proj%polar%latin1))//" "// &
874  trim(to_char(this%grid%proj%polar%latin2)))
875 #endif
876 ! projection center flag, aka hemisphere
877  CALL grib_get(gaid,'projectionCenterFlag',&
878  this%grid%proj%polar%projection_center_flag, ierr)
879  IF (ierr /= grib_success) THEN ! try center/centre
880  CALL grib_get(gaid,'projectionCentreFlag',&
881  this%grid%proj%polar%projection_center_flag)
882  ENDIF
883 
884  IF (iand(this%grid%proj%polar%projection_center_flag,64) == 1) THEN
885  CALL l4f_category_log(this%category,l4f_error, &
886  "griddim_import_gribapi, bi-polar projections not supported")
887  CALL raise_error()
888  ENDIF
889 ! line of view, aka central meridian
890  CALL grib_get(gaid,'LoVInDegrees',this%grid%proj%lov)
891 #ifdef DEBUG
892  CALL l4f_category_log(this%category,l4f_debug, &
893  "griddim_import_gribapi, central meridian "//trim(to_char(this%grid%proj%lov)))
894 #endif
895 
896 ! latitude at which dx and dy are valid
897  IF (editionnumber == 1) THEN
898 ! ECMWF (gribex/grib_api) says: Grid lengths are in metres, at the
899 ! 60-degree parallel nearest to the pole on the projection plane.
900 ! IF (IAND(this%projection_center_flag, 128) == 0) THEN
901 ! this%grid%proj%polar%lad = 60.D0
902 ! ELSE
903 ! this%grid%proj%polar%lad = -60.D0
904 ! ENDIF
905 ! WMO says: Grid lengths are in units of metres, at the secant cone
906 ! intersection parallel nearest to the pole on the projection plane.
907  this%grid%proj%polar%lad = this%grid%proj%polar%latin1
908  ELSE IF (editionnumber == 2) THEN
909  CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
910  ENDIF
911 #ifdef DEBUG
912  CALL l4f_category_log(this%category,l4f_debug, &
913  "griddim_import_gribapi, lad "//trim(to_char(this%grid%proj%polar%lad)))
914 #endif
915 
916 ! compute projected extremes from lon and lat of first point
917  CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
918  CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
919  CALL long_reset_0_360(lofirst)
920  CALL long_reset_to_cart_closest(this%grid%proj%lov, lofirst)
921 #ifdef DEBUG
922  CALL l4f_category_log(this%category,l4f_debug, &
923  "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
924  CALL l4f_category_log(this%category,l4f_debug, &
925  "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
926 #endif
927 
928  CALL proj(this, lofirst, lafirst, x1, y1)
929  IF (iscansnegatively == 0) THEN
930  this%grid%grid%xmin = x1
931  this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
932  ELSE
933  this%grid%grid%xmax = x1
934  this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
935  ENDIF
936  IF (jscanspositively == 0) THEN
937  this%grid%grid%ymax = y1
938  this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
939  ELSE
940  this%grid%grid%ymin = y1
941  this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
942  ENDIF
943 ! keep these values for personal pleasure
944  this%grid%proj%polar%lon1 = lofirst
945  this%grid%proj%polar%lat1 = lafirst
946 
947 ! Keys for equatorial projections
948 CASE ('mercator')
949 
950  CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
951  CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
952  CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
953  this%grid%proj%lov = 0.0d0 ! not in grib
954 
955  CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
956  CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
957 
958  CALL proj(this, lofirst, lafirst, x1, y1)
959  IF (iscansnegatively == 0) THEN
960  this%grid%grid%xmin = x1
961  this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
962  ELSE
963  this%grid%grid%xmax = x1
964  this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
965  ENDIF
966  IF (jscanspositively == 0) THEN
967  this%grid%grid%ymax = y1
968  this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
969  ELSE
970  this%grid%grid%ymin = y1
971  this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
972  ENDIF
973 
974  IF (editionnumber == 2) THEN
975  CALL grib_get(gaid,'orientationOfTheGridInDegrees',orient)
976  IF (orient /= 0.0d0) THEN
977  CALL l4f_category_log(this%category,l4f_error, &
978  "griddim_import_gribapi, Mercator grid orientation != 0 not supported")
979  CALL raise_error()
980  ENDIF
981  ENDIF
982 
983 #ifdef DEBUG
984  CALL unproj(this, x1, y1, lofirst, lafirst)
985  CALL l4f_category_log(this%category,l4f_debug, &
986  "griddim_import_gribapi, unprojected first point "//t2c(lofirst)//" "// &
987  t2c(lafirst))
988 
989  CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
990  CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
991  CALL proj(this, lolast, lalast, x1, y1)
992  CALL l4f_category_log(this%category,l4f_debug, &
993  "griddim_import_gribapi, extremes from grib "//t2c(x1)//" "//t2c(y1))
994  CALL l4f_category_log(this%category,l4f_debug, &
995  "griddim_import_gribapi, extremes from proj "//t2c(this%grid%grid%xmin)// &
996  " "//t2c(this%grid%grid%ymin)//" "//t2c(this%grid%grid%xmax)//" "// &
997  t2c(this%grid%grid%ymax))
998 #endif
999 
1000 CASE ('UTM')
1001 
1002  CALL grib_get(gaid,'zone',zone)
1003 
1004  CALL grib_get(gaid,'datum',datum)
1005  IF (datum == 0) THEN
1006  CALL grib_get(gaid,'referenceLongitude',reflon)
1007  CALL grib_get(gaid,'falseEasting',this%grid%proj%xoff)
1008  CALL grib_get(gaid,'falseNorthing',this%grid%proj%yoff)
1009  CALL set_val(this%grid%proj, zone=zone, lov=reflon/1.0d6)
1010  ELSE
1011  CALL l4f_category_log(this%category,l4f_error,'only datum 0 supported')
1012  CALL raise_fatal_error()
1013  ENDIF
1014 
1015  CALL grib_get(gaid,'eastingOfFirstGridPoint',lofirst)
1016  CALL grib_get(gaid,'eastingOfLastGridPoint',lolast)
1017  CALL grib_get(gaid,'northingOfFirstGridPoint',lafirst)
1018  CALL grib_get(gaid,'northingOfLastGridPoint',lalast)
1019 
1020  IF (iscansnegatively == 0) THEN
1021  this%grid%grid%xmin = lofirst
1022  this%grid%grid%xmax = lolast
1023  ELSE
1024  this%grid%grid%xmax = lofirst
1025  this%grid%grid%xmin = lolast
1026  ENDIF
1027  IF (jscanspositively == 0) THEN
1028  this%grid%grid%ymax = lafirst
1029  this%grid%grid%ymin = lalast
1030  ELSE
1031  this%grid%grid%ymin = lafirst
1032  this%grid%grid%ymax = lalast
1033  ENDIF
1034 
1035 ! compute dx and dy (should we get them from grib?)
1036  CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
1037 
1038 CASE default
1039  CALL l4f_category_log(this%category,l4f_error, &
1040  "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
1041  CALL raise_error()
1042 
1043 END SELECT
1044 
1045 CONTAINS
1046 ! utilities routines for grib_api, is there a better place?
1047 SUBROUTINE grib_get_dmiss(gaid, key, value)
1048 INTEGER,INTENT(in) :: gaid
1049 CHARACTER(len=*),INTENT(in) :: key
1050 DOUBLE PRECISION,INTENT(out) :: value
1051 
1052 INTEGER :: ierr
1053 
1054 CALL grib_get(gaid, key, value, ierr)
1055 IF (ierr /= grib_success) value = dmiss
1056 
1057 END SUBROUTINE grib_get_dmiss
1058 
1059 SUBROUTINE grib_get_imiss(gaid, key, value)
1060 INTEGER,INTENT(in) :: gaid
1061 CHARACTER(len=*),INTENT(in) :: key
1062 INTEGER,INTENT(out) :: value
1063 
1064 INTEGER :: ierr
1065 
1066 CALL grib_get(gaid, key, value, ierr)
1067 IF (ierr /= grib_success) value = imiss
1068 
1069 END SUBROUTINE grib_get_imiss
1070 
1071 
1072 SUBROUTINE griddim_import_ellipsoid(this, gaid)
1073 TYPE(griddim_def),INTENT(inout) :: this
1074 INTEGER,INTENT(in) :: gaid
1075 
1076 INTEGER :: shapeofearth, iv, is
1077 DOUBLE PRECISION :: r1, r2
1078 
1079 IF (editionnumber == 2) THEN
1080  CALL grib_get(gaid, 'shapeOfTheEarth', shapeofearth)
1081  SELECT CASE(shapeofearth)
1082  CASE(0) ! spherical
1083  CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
1084  CASE(1) ! spherical generic
1085  CALL grib_get(gaid, 'scaleFactorOfRadiusOfSphericalEarth', is)
1086  CALL grib_get(gaid, 'scaledValueOfRadiusOfSphericalEarth', iv)
1087  r1 = dble(iv) / 10**is
1088  CALL set_val(this, ellips_smaj_axis=r1, ellips_flatt=0.0d0)
1089  CASE(2) ! iau65
1090  CALL set_val(this, ellips_smaj_axis=6378160.0d0, ellips_flatt=1.0d0/297.0d0)
1091  CASE(3,7) ! ellipsoidal generic
1092  CALL grib_get(gaid, 'scaleFactorOfEarthMajorAxis', is)
1093  CALL grib_get(gaid, 'scaledValueOfEarthMajorAxis', iv)
1094  r1 = dble(iv) / 10**is
1095  CALL grib_get(gaid, 'scaleFactorOfEarthMinorAxis', is)
1096  CALL grib_get(gaid, 'scaledValueOfEarthMinorAxis', iv)
1097  r2 = dble(iv) / 10**is
1098  IF (shapeofearth == 3) THEN ! km->m
1099  r1 = r1*1000.0d0
1100  r2 = r2*1000.0d0
1101  ENDIF
1102  IF (abs(r1) < 1.0d-6) THEN ! suspicious data read from grib
1103  CALL l4f_category_log(this%category,l4f_warn,'zero Earth major axis '// &
1104  'read from grib, going on with spherical Earth but the results may be wrong')
1105  CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
1106  ELSE
1107  CALL set_val(this, ellips_smaj_axis=r1, ellips_flatt=(r1-r2)/r1)
1108  ENDIF
1109  CASE(4) ! iag-grs80
1110  CALL set_val(this, ellips_type=ellips_grs80)
1111  CASE(5) ! wgs84
1112  CALL set_val(this, ellips_type=ellips_wgs84)
1113  CASE(6) ! spherical
1114  CALL set_val(this, ellips_smaj_axis=6371229.0d0, ellips_flatt=0.0d0)
1115 ! CASE(7) ! google earth-like?
1116  CASE default
1117  CALL l4f_category_log(this%category,l4f_error,'shapeOfTheEarth '// &
1118  t2c(shapeofearth)//' not supported in grib2')
1119  CALL raise_fatal_error()
1120 
1121  END SELECT
1122 
1123 ELSE
1124 
1125  CALL grib_get(gaid, 'earthIsOblate', shapeofearth)
1126  IF (shapeofearth == 0) THEN ! spherical
1127  CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
1128  ELSE ! iau65
1129  CALL set_val(this, ellips_smaj_axis=6378160.0d0, ellips_flatt=1.0d0/297.0d0)
1130  ENDIF
1131 
1132 ENDIF
1133 
1134 END SUBROUTINE griddim_import_ellipsoid
1135 
1136 
1137 END SUBROUTINE griddim_import_gribapi
1138 
1139 
1140 ! grib_api driver
1141 SUBROUTINE griddim_export_gribapi(this, gaid)
1142 USE grib_api
1143 TYPE(griddim_def),INTENT(in) :: this ! griddim object
1144 INTEGER, INTENT(inout) :: gaid ! grib_api id of the grib loaded in memory to export
1145 
1146 INTEGER :: EditionNumber, iScansNegatively, jScansPositively, nv, pvl, zone, ierr
1147 DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, reflon
1148 DOUBLE PRECISION :: sdx, sdy, ratio, tol, ex, ey
1149 
1150 
1151 ! Generic keys
1152 CALL grib_get(gaid,'GRIBEditionNumber', editionnumber)
1153 ! the following required since eccodes
1154 IF (editionnumber == 2) CALL grib_set(gaid,'shapeOfTheEarth', 6)
1155 CALL grib_set(gaid,'typeOfGrid', this%grid%proj%proj_type)
1156 #ifdef DEBUG
1157 CALL l4f_category_log(this%category,l4f_debug, &
1158  "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
1159 #endif
1160 
1161 IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
1162  CALL grib_set(gaid,'numberOfDataPoints', this%dim%nx)
1163 ! reenable when possible
1164 ! CALL griddim_export_ellipsoid(this, gaid)
1165  RETURN
1166 ENDIF
1167 
1168 
1169 ! Edition dependent setup
1170 IF (editionnumber == 1) THEN
1171  ratio = 1.d3
1172 ELSE IF (editionnumber == 2) THEN
1173  ratio = 1.d6
1174 ELSE
1175  ratio = 0.0d0 ! signal error?!
1176 ENDIF
1177 
1178 ! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
1179 CALL griddim_export_ellipsoid(this, gaid)
1180 CALL grib_set(gaid, 'Ni', this%dim%nx)
1181 CALL grib_set(gaid, 'Nj', this%dim%ny)
1182 CALL grib_set(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
1183 
1184 CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
1185 CALL grib_get(gaid,'jScansPositively',jscanspositively)
1186 
1187 ! Keys for rotated grids (checked through missing values and/or error code)
1188 !SELECT CASE (this%grid%proj%proj_type)
1189 !CASE ('rotated_ll', 'stretched_rotated_ll', 'polar_stereographic', 'lambert', 'albers')
1190 CALL grib_set_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
1191  this%grid%proj%rotated%longitude_south_pole, 0.0d0)
1192 CALL grib_set_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
1193  this%grid%proj%rotated%latitude_south_pole, -90.0d0)
1194 IF (editionnumber == 1) THEN
1195  CALL grib_set_dmiss(gaid,'angleOfRotationInDegrees', &
1196  this%grid%proj%rotated%angle_rotation, 0.0d0)
1197 ELSE IF (editionnumber == 2)THEN
1198  CALL grib_set_dmiss(gaid,'angleOfRotationOfProjectionInDegrees', &
1199  this%grid%proj%rotated%angle_rotation, 0.0d0)
1200 ENDIF
1201 
1202 ! Keys for stretched grids (checked through missing values and/or error code)
1203 ! units must be verified, still experimental in grib_api
1204 ! # TODO: Is it a float? Is it signed?
1205 IF (editionnumber == 1) THEN
1206  CALL grib_set_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
1207  this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
1208  CALL grib_set_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
1209  this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
1210  CALL grib_set_dmiss(gaid,'stretchingFactor', &
1211  this%grid%proj%stretched%stretch_factor, 1.0d0)
1212 ELSE IF (editionnumber == 2) THEN
1213  CALL grib_set_dmiss(gaid,'longitudeOfThePoleOfStretching', &
1214  this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
1215  CALL grib_set_dmiss(gaid,'latitudeOfThePoleOfStretching', &
1216  this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
1217  CALL grib_set_dmiss(gaid,'stretchingFactor', &
1218  this%grid%proj%stretched%stretch_factor, 1.0d6, 1.0d6)
1219 ENDIF
1220 
1221 ! Projection-dependent keys
1222 SELECT CASE (this%grid%proj%proj_type)
1223 
1224 ! Keys for sphaerical coordinate systems
1225 CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
1226 
1227  IF (iscansnegatively == 0) THEN
1228  lofirst = this%grid%grid%xmin
1229  lolast = this%grid%grid%xmax
1230  ELSE
1231  lofirst = this%grid%grid%xmax
1232  lolast = this%grid%grid%xmin
1233  ENDIF
1234  IF (jscanspositively == 0) THEN
1235  lafirst = this%grid%grid%ymax
1236  lalast = this%grid%grid%ymin
1237  ELSE
1238  lafirst = this%grid%grid%ymin
1239  lalast = this%grid%grid%ymax
1240  ENDIF
1241 
1242 ! reset lon in standard grib 2 definition [0,360]
1243  IF (editionnumber == 1) THEN
1244  CALL long_reset_m180_360(lofirst)
1245  CALL long_reset_m180_360(lolast)
1246  ELSE IF (editionnumber == 2) THEN
1247  CALL long_reset_0_360(lofirst)
1248  CALL long_reset_0_360(lolast)
1249  ENDIF
1250 
1251  CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
1252  CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
1253  CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
1254  CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
1255 
1256 ! test relative coordinate truncation error with respect to tol
1257 ! tol should be tuned
1258  sdx = this%grid%grid%dx*ratio
1259  sdy = this%grid%grid%dy*ratio
1260 ! error is computed relatively to the whole grid extension
1261  ex = abs(((sdx - nint(sdx))/ratio*this%dim%nx)/ &
1262  (this%grid%grid%xmax-this%grid%grid%xmin))
1263  ey = abs(((sdy - nint(sdy))/ratio*this%dim%ny)/ &
1264  (this%grid%grid%ymax-this%grid%grid%ymin))
1265  tol = 1.0d-3
1266  IF (ex > tol .OR. ey > tol) THEN
1267 #ifdef DEBUG
1268  CALL l4f_category_log(this%category,l4f_debug, &
1269  "griddim_export_gribapi, error on x "//&
1270  trim(to_char(ex))//"/"//t2c(tol))
1271  CALL l4f_category_log(this%category,l4f_debug, &
1272  "griddim_export_gribapi, error on y "//&
1273  trim(to_char(ey))//"/"//t2c(tol))
1274 #endif
1275 ! previous frmula relative to a single grid step
1276 ! tol = 1.0d0/ratio
1277 ! IF (ABS(NINT(sdx)/sdx - 1.0d0) > tol .OR. ABS(NINT(sdy)/sdy - 1.0d0) > tol) THEN
1278 !#ifdef DEBUG
1279 ! CALL l4f_category_log(this%category,L4F_DEBUG, &
1280 ! "griddim_export_gribapi, dlon relative error: "//&
1281 ! TRIM(to_char(ABS(NINT(sdx)/sdx - 1.0d0)))//">"//TRIM(to_char(tol)))
1282 ! CALL l4f_category_log(this%category,L4F_DEBUG, &
1283 ! "griddim_export_gribapi, dlat relative error: "//&
1284 ! TRIM(to_char(ABS(NINT(sdy)/sdy - 1.0d0)))//">"//TRIM(to_char(tol)))
1285 !#endif
1286  CALL l4f_category_log(this%category,l4f_info, &
1287  "griddim_export_gribapi, increments not given: inaccurate!")
1288  CALL grib_set_missing(gaid,'Di')
1289  CALL grib_set_missing(gaid,'Dj')
1290  CALL grib_set(gaid,'ijDirectionIncrementGiven',0)
1291  ELSE
1292 #ifdef DEBUG
1293  CALL l4f_category_log(this%category,l4f_debug, &
1294  "griddim_export_gribapi, setting increments: "// &
1295  trim(to_char(this%grid%grid%dx))//' '//trim(to_char(this%grid%grid%dy)))
1296 #endif
1297  CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
1298  CALL grib_set(gaid,'iDirectionIncrement',nint(this%grid%grid%dx*ratio))
1299  CALL grib_set(gaid,'jDirectionIncrement',nint(this%grid%grid%dy*ratio))
1300 ! this does not work in grib_set
1301 ! CALL grib_set(gaid,'iDirectionIncrementInDegrees',this%grid%grid%dx)
1302 ! CALL grib_set(gaid,'jDirectionIncrementInDegrees',this%grid%grid%dy)
1303  ENDIF
1304 
1305 ! Keys for polar projections
1306 CASE ('polar_stereographic', 'lambert', 'albers')
1307 ! increments are required
1308  CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
1309  CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
1310  CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
1311 ! latin1/latin2 may be missing (e.g. stereographic)
1312  CALL grib_set_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
1313  CALL grib_set_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
1314 ! projection center flag, aka hemisphere
1315  CALL grib_set(gaid,'projectionCenterFlag',&
1316  this%grid%proj%polar%projection_center_flag, ierr)
1317  IF (ierr /= grib_success) THEN ! try center/centre
1318  CALL grib_set(gaid,'projectionCentreFlag',&
1319  this%grid%proj%polar%projection_center_flag)
1320  ENDIF
1321 
1322 
1323 ! line of view, aka central meridian
1324  CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
1325 ! latitude at which dx and dy are valid
1326  IF (editionnumber == 2) THEN
1327  CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
1328  ENDIF
1329 
1330 ! compute lon and lat of first point from projected extremes
1331  CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
1332  lofirst, lafirst)
1333 ! reset lon in standard grib 2 definition [0,360]
1334  IF (editionnumber == 1) THEN
1335  CALL long_reset_m180_360(lofirst)
1336  ELSE IF (editionnumber == 2) THEN
1337  CALL long_reset_0_360(lofirst)
1338  ENDIF
1339  CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
1340  CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
1341 
1342 ! Keys for equatorial projections
1343 CASE ('mercator')
1344 
1345 ! increments are required
1346  CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
1347  CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
1348  CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
1349 
1350 ! line of view, aka central meridian, not in grib
1351 ! CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
1352 ! latitude at which dx and dy are valid
1353  CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
1354 
1355 ! compute lon and lat of first and last points from projected extremes
1356  CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
1357  lofirst, lafirst)
1358  CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
1359  CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
1360  CALL get_unproj_last(this, iscansnegatively, jscanspositively, &
1361  lolast, lalast)
1362  CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
1363  CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
1364 
1365  IF (editionnumber == 2) THEN
1366  CALL grib_set(gaid,'orientationOfTheGridInDegrees',0.)
1367  ENDIF
1368 
1369 CASE ('UTM')
1370 
1371  CALL grib_set(gaid,'datum',0)
1372  CALL get_val(this, zone=zone, lov=reflon)
1373  CALL grib_set(gaid,'referenceLongitude',nint(reflon/1.0d6))
1374  CALL grib_set(gaid,'falseEasting',this%grid%proj%xoff)
1375  CALL grib_set(gaid,'falseNorthing',this%grid%proj%yoff)
1376 
1377  CALL grib_set(gaid,'iDirectionIncrement',this%grid%grid%dx)
1378  CALL grib_set(gaid,'jDirectionIncrement',this%grid%grid%dy)
1379 
1380 !res/scann ??
1381 
1382  CALL grib_set(gaid,'zone',zone)
1383 
1384  IF (iscansnegatively == 0) THEN
1385  lofirst = this%grid%grid%xmin
1386  lolast = this%grid%grid%xmax
1387  ELSE
1388  lofirst = this%grid%grid%xmax
1389  lolast = this%grid%grid%xmin
1390  ENDIF
1391  IF (jscanspositively == 0) THEN
1392  lafirst = this%grid%grid%ymax
1393  lalast = this%grid%grid%ymin
1394  ELSE
1395  lafirst = this%grid%grid%ymin
1396  lalast = this%grid%grid%ymax
1397  ENDIF
1398 
1399  CALL grib_set(gaid,'eastingOfFirstGridPoint',lofirst)
1400  CALL grib_set(gaid,'eastingOfLastGridPoint',lolast)
1401  CALL grib_set(gaid,'northingOfFirstGridPoint',lafirst)
1402  CALL grib_set(gaid,'northingOfLastGridPoint',lalast)
1403 
1404 CASE default
1405  CALL l4f_category_log(this%category,l4f_error, &
1406  "griddim_export_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
1407  CALL raise_error()
1408 
1409 END SELECT
1410 
1411 ! hack for position of vertical coordinate parameters
1412 ! buggy in grib_api
1413 IF (editionnumber == 1) THEN
1414 ! CALL grib_get(gaid,"PVPresent",pvp) ! alias, probably useless
1415  CALL grib_get(gaid,"NV",nv)
1416 #ifdef DEBUG
1417  CALL l4f_category_log(this%category,l4f_debug,"griddim_export_gribapi, coding "// &
1418  trim(to_char(nv))//" vertical coordinate parameters")
1419 #endif
1420 
1421  IF (nv == 0) THEN
1422  pvl = 255
1423  ELSE
1424  SELECT CASE (this%grid%proj%proj_type)
1425  CASE ('regular_ll') ! check whether "29-32 Set to zero (reserved)" are required
1426  pvl = 33
1427  CASE ('polar_stereographic')
1428  pvl = 33
1429  CASE ('rotated_ll', 'stretched_ll', 'lambert', 'albers')
1430  pvl = 43
1431  CASE ('stretched_rotated_ll')
1432  pvl = 43
1433  CASE DEFAULT
1434  pvl = 43 !?
1435  END SELECT
1436  ENDIF
1437 
1438  CALL grib_set(gaid,"pvlLocation",pvl)
1439 #ifdef DEBUG
1440  CALL l4f_category_log(this%category,l4f_debug,"griddim_export_gribapi, coding "// &
1441  trim(to_char(pvl))//" as vertical coordinate parameter location")
1442 #endif
1443 
1444 ENDIF
1445 
1446 
1447 CONTAINS
1448 ! utilities routines for grib_api, is there a better place?
1449 SUBROUTINE grib_set_dmiss(gaid, key, val, default, factor)
1450 INTEGER,INTENT(in) :: gaid
1451 CHARACTER(len=*),INTENT(in) :: key
1452 DOUBLE PRECISION,INTENT(in) :: val
1453 DOUBLE PRECISION,INTENT(in),OPTIONAL :: default
1454 DOUBLE PRECISION,INTENT(in),OPTIONAL :: factor
1455 
1456 INTEGER :: ierr
1457 
1458 IF (c_e(val)) THEN
1459  IF (PRESENT(factor)) THEN
1460  CALL grib_set(gaid, key, val*factor, ierr)
1461  ELSE
1462  CALL grib_set(gaid, key, val, ierr)
1463  ENDIF
1464 ELSE IF (PRESENT(default)) THEN
1465  CALL grib_set(gaid, key, default, ierr)
1466 ENDIF
1467 
1468 END SUBROUTINE grib_set_dmiss
1469 
1470 SUBROUTINE grib_set_imiss(gaid, key, value, default)
1471 INTEGER,INTENT(in) :: gaid
1472 CHARACTER(len=*),INTENT(in) :: key
1473 INTEGER,INTENT(in) :: value
1474 INTEGER,INTENT(in),OPTIONAL :: default
1475 
1476 INTEGER :: ierr
1477 
1478 IF (c_e(value)) THEN
1479  CALL grib_set(gaid, key, value, ierr)
1480 ELSE IF (PRESENT(default)) THEN
1481  CALL grib_set(gaid, key, default, ierr)
1482 ENDIF
1483 
1484 END SUBROUTINE grib_set_imiss
1485 
1486 SUBROUTINE griddim_export_ellipsoid(this, gaid)
1487 TYPE(griddim_def),INTENT(in) :: this
1488 INTEGER,INTENT(in) :: gaid
1489 
1490 INTEGER :: ellips_type, ierr
1491 DOUBLE PRECISION :: r1, r2, f
1492 
1493 CALL get_val(this, ellips_smaj_axis=r1, ellips_flatt=f, ellips_type=ellips_type)
1494 
1495 IF (editionnumber == 2) THEN
1496 
1497 ! why does it produce a message even with ierr?
1498  CALL grib_set_missing(gaid, 'scaleFactorOfRadiusOfSphericalEarth', ierr)
1499  CALL grib_set_missing(gaid, 'scaledValueOfRadiusOfSphericalEarth', ierr)
1500  CALL grib_set_missing(gaid, 'scaleFactorOfEarthMajorAxis', ierr)
1501  CALL grib_set_missing(gaid, 'scaledValueOfEarthMajorAxis', ierr)
1502  CALL grib_set_missing(gaid, 'scaleFactorOfEarthMinorAxis', ierr)
1503  CALL grib_set_missing(gaid, 'scaledValueOfEarthMinorAxis', ierr)
1504 
1505  SELECT CASE(ellips_type)
1506  CASE(ellips_grs80) ! iag-grs80
1507  CALL grib_set(gaid, 'shapeOfTheEarth', 4)
1508  CASE(ellips_wgs84) ! wgs84
1509  CALL grib_set(gaid, 'shapeOfTheEarth', 5)
1510  CASE default
1511  IF (f == 0.0d0) THEN ! spherical Earth
1512  IF (r1 == 6367470.0d0) THEN ! spherical
1513  CALL grib_set(gaid, 'shapeOfTheEarth', 0)
1514  ELSE IF (r1 == 6371229.0d0) THEN ! spherical
1515  CALL grib_set(gaid, 'shapeOfTheEarth', 6)
1516  ELSE ! spherical generic
1517  CALL grib_set(gaid, 'shapeOfTheEarth', 1)
1518  CALL grib_set(gaid, 'scaleFactorOfRadiusOfSphericalEarth', 2)
1519  CALL grib_set(gaid, 'scaledValueOfRadiusOfSphericalEarth', int(r1*100.0d0))
1520  ENDIF
1521  ELSE ! ellipsoidal
1522  IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
1523  CALL grib_set(gaid, 'shapeOfTheEarth', 2)
1524  ELSE ! ellipsoidal generic
1525  CALL grib_set(gaid, 'shapeOfTheEarth', 3)
1526  r2 = r1*(1.0d0 - f)
1527  CALL grib_set(gaid, 'scaleFactorOfEarthMajorAxis', 5)
1528  CALL grib_set(gaid, 'scaledValueOfEarthMajorAxis', &
1529  int(r1*100.0d0))
1530  CALL grib_set(gaid, 'scaleFactorOfEarthMinorAxis', 5)
1531  CALL grib_set(gaid, 'scaledValueOfEarthMinorAxis', &
1532  int(r2*100.0d0))
1533  ENDIF
1534  ENDIF
1535  END SELECT
1536 
1537 ELSE
1538 
1539  IF (r1 == 6367470.0d0 .AND. f == 0.0d0) THEN ! spherical
1540  CALL grib_set(gaid, 'earthIsOblate', 0)
1541  ELSE IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
1542  CALL grib_set(gaid, 'earthIsOblate', 1)
1543  ELSE IF (f == 0.0d0) THEN ! generic spherical
1544  CALL grib_set(gaid, 'earthIsOblate', 0)
1545  CALL l4f_category_log(this%category,l4f_warn,'desired spherical Earth radius &
1546  &not supported in grib 1, coding with default radius of 6367470 m')
1547  ELSE ! generic ellipsoidal
1548  CALL grib_set(gaid, 'earthIsOblate', 1)
1549  CALL l4f_category_log(this%category,l4f_warn,'desired Earth ellipsoid &
1550  &not supported in grib 1, coding with default iau65 ellipsoid')
1551  ENDIF
1552 
1553 ENDIF
1554 
1555 END SUBROUTINE griddim_export_ellipsoid
1556 
1557 SUBROUTINE get_unproj_first(this, iScansNegatively, jScansPositively, &
1558  loFirst, laFirst)
1559 TYPE(griddim_def),INTENT(in) :: this ! griddim object
1560 INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
1561 DOUBLE PRECISION,INTENT(out) :: loFirst, laFirst
1562 
1563 ! compute lon and lat of first point from projected extremes
1564 IF (iscansnegatively == 0) THEN
1565  IF (jscanspositively == 0) THEN
1566  CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymax, lofirst, lafirst)
1567  ELSE
1568  CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymin, lofirst, lafirst)
1569  ENDIF
1570 ELSE
1571  IF (jscanspositively == 0) THEN
1572  CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymax, lofirst, lafirst)
1573  ELSE
1574  CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymin, lofirst, lafirst)
1575  ENDIF
1576 ENDIF
1577 ! use the values kept for personal pleasure ?
1578 ! loFirst = this%grid%proj%polar%lon1
1579 ! laFirst = this%grid%proj%polar%lat1
1580 END SUBROUTINE get_unproj_first
1581 
1582 SUBROUTINE get_unproj_last(this, iScansNegatively, jScansPositively, &
1583  loLast, laLast)
1584 TYPE(griddim_def),INTENT(in) :: this ! griddim object
1585 INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
1586 DOUBLE PRECISION,INTENT(out) :: loLast, laLast
1587 
1588 ! compute lon and lat of last point from projected extremes
1589 IF (iscansnegatively == 0) THEN
1590  IF (jscanspositively == 0) THEN
1591  CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymin, lolast, lalast)
1592  ELSE
1593  CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymax, lolast, lalast)
1594  ENDIF
1595 ELSE
1596  IF (jscanspositively == 0) THEN
1597  CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymin, lolast, lalast)
1598  ELSE
1599  CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymax, lolast, lalast)
1600  ENDIF
1601 ENDIF
1602 ! use the values kept for personal pleasure ?
1603 ! loLast = this%grid%proj%polar%lon?
1604 ! laLast = this%grid%proj%polar%lat?
1605 END SUBROUTINE get_unproj_last
1606 
1607 END SUBROUTINE griddim_export_gribapi
1608 #endif
1609 
1610 
1611 #ifdef HAVE_LIBGDAL
1612 ! gdal driver
1613 SUBROUTINE griddim_import_gdal(this, gdalid, gdal_options)
1614 USE gdal
1615 TYPE(griddim_def),INTENT(inout) :: this ! griddim object
1616 TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal rasterband pointer
1617 TYPE(gdal_file_id_options),INTENT(in) :: gdal_options
1618 
1619 TYPE(gdaldataseth) :: hds
1620 REAL(kind=c_double) :: geotrans(6) !, invgeotrans(6), x1, y1, x2, y2
1621 INTEGER(kind=c_int) :: offsetx, offsety
1622 INTEGER :: ier
1623 
1624 hds = gdalgetbanddataset(gdalid) ! go back to dataset
1625 ier = gdalgetgeotransform(hds, geotrans)
1626 
1627 IF (ier /= 0) THEN
1628  CALL l4f_category_log(this%category, l4f_error, &
1629  'griddim_import_gdal: error in accessing gdal dataset')
1630  CALL raise_error()
1631  RETURN
1632 ENDIF
1633 IF (geotrans(3) /= 0.0_c_double .OR. geotrans(5) /= 0.0_c_double) THEN ! dataset has not a diagonal transformation
1634  CALL l4f_category_log(this%category, l4f_error, &
1635  'griddim_import_gdal: dataset has a non-diagonal transformation matrix, unsupported')
1636  CALL raise_error()
1637  RETURN
1638 ENDIF
1639 
1640 CALL gdaldatasetbbsize_f(hds, gdal_options%xmin, gdal_options%ymin, &
1641  gdal_options%xmax, gdal_options%ymax, &
1642  this%dim%nx, this%dim%ny, offsetx, offsety, &
1643  this%grid%grid%xmin, this%grid%grid%ymin, this%grid%grid%xmax, this%grid%grid%ymax)
1644 
1645 IF (this%dim%nx == 0 .OR. this%dim%ny == 0) THEN
1646  CALL l4f_category_log(this%category, l4f_warn, &
1647  'griddim_import_gdal: requested bounding box '//t2c(gdal_options%xmin)//','// &
1648  t2c(gdal_options%ymin)//','//t2c(gdal_options%xmax)//','//&
1649  t2c(gdal_options%ymax))
1650  CALL l4f_category_log(this%category, l4f_warn, &
1651  'determines an empty dataset '//t2c(this%dim%nx)//'x'//t2c(this%dim%ny))
1652 ENDIF
1653 
1654 ! get grid corners
1655 !CALL gdalapplygeotransform(geotrans, 0.5_c_double, 0.5_c_double, x1, y1)
1656 !CALL gdalapplygeotransform(geotrans, gdalgetrasterbandxsize(gdalid)-0.5_c_double, &
1657 ! gdalgetrasterbandysize(gdalid)-0.5_c_double, x2, y2)
1658 
1659 !IF (geotrans(3) == 0.0_c_double .AND. geotrans(5) == 0.0_c_double) THEN ! transformation is diagonal, no transposing
1660 ! this%dim%nx = gdalgetrasterbandxsize(gdalid)
1661 ! this%dim%ny = gdalgetrasterbandysize(gdalid)
1662 ! this%grid%grid%xmin = MIN(x1, x2)
1663 ! this%grid%grid%xmax = MAX(x1, x2)
1664 ! this%grid%grid%ymin = MIN(y1, y2)
1665 ! this%grid%grid%ymax = MAX(y1, y2)
1666 !ELSE IF (geotrans(2) == 0.0_c_double .AND. geotrans(6) == 0.0_c_double) THEN ! transformation is anti-diagonal, transposing will have to be done
1667 !
1668 ! this%dim%nx = gdalgetrasterbandysize(gdalid)
1669 ! this%dim%ny = gdalgetrasterbandxsize(gdalid)
1670 ! this%grid%grid%xmin = MIN(y1, y2)
1671 ! this%grid%grid%xmax = MAX(y1, y2)
1672 ! this%grid%grid%ymin = MIN(x1, x2)
1673 ! this%grid%grid%ymax = MAX(x1, x2)
1674 !
1675 !ELSE ! transformation is a rotation, not supported
1676 !ENDIF
1677 
1678 this%grid%proj%proj_type = 'regular_ll' ! forced, only one supported (check?)
1679 
1680 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
1681 this%grid%grid%component_flag = 0
1682 
1683 END SUBROUTINE griddim_import_gdal
1684 #endif
1685 
1686 
1688 SUBROUTINE griddim_display(this)
1689 TYPE(griddim_def),INTENT(in) :: this
1690 
1691 print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
1692 
1693 CALL display(this%grid%proj)
1694 CALL display(this%grid%grid)
1695 CALL display(this%dim)
1696 
1697 print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
1698 
1699 END SUBROUTINE griddim_display
1700 
1701 
1702 #define VOL7D_POLY_TYPE TYPE(grid_def)
1703 #define VOL7D_POLY_TYPES _grid
1704 #include "array_utilities_inc.F90"
1705 #undef VOL7D_POLY_TYPE
1706 #undef VOL7D_POLY_TYPES
1707 
1708 #define VOL7D_POLY_TYPE TYPE(griddim_def)
1709 #define VOL7D_POLY_TYPES _griddim
1710 #include "array_utilities_inc.F90"
1711 #undef VOL7D_POLY_TYPE
1712 #undef VOL7D_POLY_TYPES
1713 
1714 
1726 SUBROUTINE griddim_wind_unrot(this, rot_mat)
1728 TYPE(griddim_def), INTENT(in) :: this
1729 DOUBLE PRECISION, POINTER :: rot_mat(:,:,:)
1730 
1731 DOUBLE PRECISION :: dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j
1732 DOUBLE PRECISION :: lat_factor
1733 INTEGER :: i, j, ip1, im1, jp1, jm1;
1734 
1735 IF (this%dim%nx <= 0 .OR. this%dim%ny <= 0 .OR. &
1736  .NOT. ASSOCIATED(this%dim%lon) .OR. .NOT. ASSOCIATED(this%dim%lat)) THEN
1737  CALL l4f_category_log(this%category, l4f_error, 'rotate_uv must be called after correct init of griddim object')
1738  NULLIFY(rot_mat)
1739  RETURN
1740 ENDIF
1741 
1742 ALLOCATE(rot_mat(this%dim%nx, this%dim%ny, 4))
1743 
1744 DO j = 1, this%dim%ny
1745  jp1 = min(j+1, this%dim%ny)
1746  jm1 = max(j-1, 1)
1747  DO i = 1, this%dim%nx
1748  ip1 = min(i+1, this%dim%nx)
1749  im1 = max(i-1, 1)
1750 
1751  dlat_i = this%dim%lat(ip1,j) - this%dim%lat(im1,j)
1752  dlat_j = this%dim%lat(i,jp1) - this%dim%lat(i,jm1)
1753 
1754  dlon_i = this%dim%lon(ip1,j) - this%dim%lon(im1,j)
1755 ! if ( dlon_i > pi ) dlon_i -= 2*pi;
1756 ! if ( dlon_i < (-pi) ) dlon_i += 2*pi;
1757  dlon_j = this%dim%lon(i,jp1) - this%dim%lon(i,jm1)
1758 ! if ( dlon_j > pi ) dlon_j -= 2*pi;
1759 ! if ( dlon_j < (-pi) ) dlon_j += 2*pi;
1760 
1761 ! check whether this is really necessary !!!!
1762  lat_factor = cos(degrad*this%dim%lat(i,j))
1763  dlon_i = dlon_i * lat_factor
1764  dlon_j = dlon_j * lat_factor
1765 
1766  dist_i = sqrt(dlon_i*dlon_i + dlat_i*dlat_i)
1767  dist_j = sqrt(dlon_j*dlon_j + dlat_j*dlat_j)
1768 
1769  IF (dist_i > 0.d0) THEN
1770  rot_mat(i,j,1) = dlon_i/dist_i
1771  rot_mat(i,j,2) = dlat_i/dist_i
1772  ELSE
1773  rot_mat(i,j,1) = 0.0d0
1774  rot_mat(i,j,2) = 0.0d0
1775  ENDIF
1776  IF (dist_j > 0.d0) THEN
1777  rot_mat(i,j,3) = dlon_j/dist_j
1778  rot_mat(i,j,4) = dlat_j/dist_j
1779  ELSE
1780  rot_mat(i,j,3) = 0.0d0
1781  rot_mat(i,j,4) = 0.0d0
1782  ENDIF
1783 
1784  ENDDO
1785 ENDDO
1786 
1787 END SUBROUTINE griddim_wind_unrot
1788 
1789 
1790 ! compute zoom indices from geographical zoom coordinates
1791 SUBROUTINE griddim_zoom_coord(this, ilon, ilat, flon, flat, ix, iy, fx, fy)
1792 TYPE(griddim_def),INTENT(in) :: this
1793 DOUBLE PRECISION,INTENT(in) :: ilon, ilat, flon, flat
1794 INTEGER,INTENT(out) :: ix, iy, fx, fy
1795 
1796 DOUBLE PRECISION :: ix1, iy1, fx1, fy1
1797 
1798 ! compute projected coordinates of vertices of desired lonlat rectangle
1799 CALL proj(this, ilon, ilat, ix1, iy1)
1800 CALL proj(this, flon, flat, fx1, fy1)
1801 
1802 CALL griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
1803 
1804 END SUBROUTINE griddim_zoom_coord
1805 
1806 
1807 ! compute zoom indices from projected zoom coordinates
1808 SUBROUTINE griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
1809 TYPE(griddim_def),INTENT(in) :: this
1810 DOUBLE PRECISION,INTENT(in) :: ix1, iy1, fx1, fy1
1811 INTEGER,INTENT(out) :: ix, iy, fx, fy
1812 
1813 INTEGER :: lix, liy, lfx, lfy
1814 
1815 ! compute projected indices
1816 lix = nint((ix1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
1817 liy = nint((iy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
1818 lfx = nint((fx1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
1819 lfy = nint((fy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
1820 ! swap projected indices, in case grid is "strongly rotated"
1821 ix = min(lix, lfx)
1822 fx = max(lix, lfx)
1823 iy = min(liy, lfy)
1824 fy = max(liy, lfy)
1825 
1826 END SUBROUTINE griddim_zoom_projcoord
1827 
1828 
1832 SUBROUTINE long_reset_0_360(lon)
1833 DOUBLE PRECISION,INTENT(inout) :: lon
1834 
1835 IF (.NOT.c_e(lon)) RETURN
1836 DO WHILE(lon < 0.0d0)
1837  lon = lon + 360.0d0
1838 END DO
1839 DO WHILE(lon >= 360.0d0)
1840  lon = lon - 360.0d0
1841 END DO
1842 
1843 END SUBROUTINE long_reset_0_360
1844 
1845 
1849 SUBROUTINE long_reset_m180_360(lon)
1850 DOUBLE PRECISION,INTENT(inout) :: lon
1851 
1852 IF (.NOT.c_e(lon)) RETURN
1853 DO WHILE(lon < -180.0d0)
1854  lon = lon + 360.0d0
1855 END DO
1856 DO WHILE(lon >= 360.0d0)
1857  lon = lon - 360.0d0
1858 END DO
1859 
1860 END SUBROUTINE long_reset_m180_360
1861 
1862 
1866 !SUBROUTINE long_reset_m90_270(lon)
1867 !DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
1868 !
1869 !IF (.NOT.c_e(lon)) RETURN
1870 !DO WHILE(lon < -90.0D0)
1871 ! lon = lon + 360.0D0
1872 !END DO
1873 !DO WHILE(lon >= 270.0D0)
1874 ! lon = lon - 360.0D0
1875 !END DO
1876 !
1877 !END SUBROUTINE long_reset_m90_270
1878 
1879 
1883 SUBROUTINE long_reset_m180_180(lon)
1884 DOUBLE PRECISION,INTENT(inout) :: lon
1885 
1886 IF (.NOT.c_e(lon)) RETURN
1887 DO WHILE(lon < -180.0d0)
1888  lon = lon + 360.0d0
1889 END DO
1890 DO WHILE(lon >= 180.0d0)
1891  lon = lon - 360.0d0
1892 END DO
1893 
1894 END SUBROUTINE long_reset_m180_180
1895 
1896 
1897 SUBROUTINE long_reset_to_cart_closest(lon, lonref)
1898 DOUBLE PRECISION,INTENT(inout) :: lon
1899 DOUBLE PRECISION,INTENT(in) :: lonref
1900 
1901 IF (.NOT.c_e(lon) .OR. .NOT.c_e(lonref)) RETURN
1902 IF (abs(lon-lonref) < 180.0d0) RETURN ! nothing to do
1903 lon = lon - nint((lon-lonref)/360.0d0)*360.0d0 ! se non e` vera e` ben trovata
1904 
1905 END SUBROUTINE long_reset_to_cart_closest
1906 
1907 
1908 END MODULE grid_class
1909 
Method for testing the existence of the object.
Copy an object, creating a fully new instance.
Definition: grid_class.F90:308
Destructors of the corresponding objects.
Definition: grid_class.F90:303
Print a brief description on stdout.
Definition: grid_class.F90:355
Export griddim object to grid_id.
Definition: grid_class.F90:350
Method for returning the contents of the object.
Definition: grid_class.F90:325
Import griddim object from grid_id.
Definition: grid_class.F90:345
Constructors of the corresponding objects.
Definition: grid_class.F90:298
Compute forward coordinate transformation from geographical system to projected system.
Definition: grid_class.F90:314
Read the object from a formatted or unformatted file.
Definition: grid_class.F90:340
Method for setting the contents of the object.
Definition: grid_class.F90:330
Compute backward coordinate transformation from projected system to geographical system.
Definition: grid_class.F90:320
Write the object on a formatted or unformatted file.
Definition: grid_class.F90:335
Index method.
Emit log message for a category with specific priority.
Costanti fisiche (DOUBLEPRECISION).
Definition: phys_const.f90:72
Gestione degli errori.
Module for describing geographically referenced regular grids.
Definition: grid_class.F90:243
Module for defining the extension and coordinates of a rectangular georeferenced grid.
This module defines an abstract interface to different drivers for access to files containing gridded...
classe per la gestione del logging
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
This object, mainly for internal use, describes a grid on a geographical projection,...
Definition: grid_class.F90:263
This object completely describes a grid on a geographic projection.
Definition: grid_class.F90:276
Derived type describing the extension of a grid and the geographical coordinates of each point.
Derived type containing driver-specific options for gdal.
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...

Generated with Doxygen.