libsim  Versione7.2.6
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
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
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 
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)
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
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 
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
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 
Export griddim object to grid_id.
Definition: grid_class.F90:354
Method for setting the contents of the object.
Definition: grid_class.F90:334
This object, mainly for internal use, describes a grid on a geographical projection, except the grid dimensions.
Definition: grid_class.F90:267
Compute forward coordinate transformation from geographical system to projected system.
Definition: grid_class.F90:318
Constructors of the corresponding objects.
Definition: grid_class.F90:302
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Print a brief description on stdout.
Definition: grid_class.F90:359
This module defines an abstract interface to different drivers for access to files containing gridded...
Module for describing geographically referenced regular grids.
Definition: grid_class.F90:247
Write the object on a formatted or unformatted file.
Definition: grid_class.F90:339
Compute backward coordinate transformation from projected system to geographical system.
Definition: grid_class.F90:324
This object completely describes a grid on a geographic projection.
Definition: grid_class.F90:280
Method for returning the contents of the object.
Definition: grid_class.F90:329
Index method.
Gestione degli errori.
Definitions of constants and functions for working with missing values.
Costanti fisiche (DOUBLEPRECISION).
Definition: phys_const.f90:72
Module for defining the extension and coordinates of a rectangular georeferenced grid.
classe per la gestione del logging
Method for testing the existence of the object.
Copy an object, creating a fully new instance.
Definition: grid_class.F90:312
Read the object from a formatted or unformatted file.
Definition: grid_class.F90:344
Emit log message for a category with specific priority.
Derived type describing the extension of a grid and the geographical coordinates of each point...
Destructors of the corresponding objects.
Definition: grid_class.F90:307
Import griddim object from grid_id.
Definition: grid_class.F90:349

Generated with Doxygen.