libsim  Versione6.3.0
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)
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)
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
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
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(lov, lonref)
582  CALL set_val(this%grid%proj, longitude_south_pole=lonref)
583  ENDIF
584  ELSE
585  lon = modulo(lonsp + 180.0d0, 360.0d0)
586 ! IF (PRESENT(lonref)) THEN
587 ! CALL long_reset_to_cart_closest(lov, lonref)
588 ! CALL set_val(this%grid%proj, longitude_south_pole=lonref)
589 ! ENDIF
590  ENDIF
591  CASE default ! use real grid limits
592  IF (c_e(this%grid%grid%xmin) .AND. c_e(this%grid%grid%xmin)) THEN
593  lon = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
594  ENDIF
595  IF (PRESENT(lonref)) THEN
596  londelta = lon
597  CALL long_reset_to_cart_closest(londelta, lonref)
598  londelta = londelta - lon
599  this%grid%grid%xmin = this%grid%grid%xmin + londelta
600  this%grid%grid%xmax = this%grid%grid%xmax + londelta
601  ENDIF
602  END SELECT
603 ENDIF
604 
605 END SUBROUTINE griddim_pistola_central_lon
606 
607 
611 SUBROUTINE griddim_gen_coord(this, x, y)
612 TYPE(griddim_def),INTENT(in) :: this
613 DOUBLE PRECISION,INTENT(out) :: x(:,:)
614 DOUBLE PRECISION,INTENT(out) :: y(:,:)
615 
616 
617 CALL grid_rect_coordinates(this%grid%grid, x, y)
618 
619 END SUBROUTINE griddim_gen_coord
620 
621 
623 SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
624 TYPE(griddim_def), INTENT(in) :: this
625 INTEGER,INTENT(in) :: nx
626 INTEGER,INTENT(in) :: ny
627 DOUBLE PRECISION,INTENT(out) :: dx
628 DOUBLE PRECISION,INTENT(out) :: dy
629 
630 CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
631 
632 END SUBROUTINE griddim_steps
633 
634 
636 SUBROUTINE griddim_setsteps(this)
637 TYPE(griddim_def), INTENT(inout) :: this
638 
639 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
640 
641 END SUBROUTINE griddim_setsteps
642 
643 
644 ! TODO
645 ! bisogna sviluppare gli altri operatori
646 ELEMENTAL FUNCTION grid_eq(this, that) RESULT(res)
647 TYPE(grid_def),INTENT(IN) :: this, that
648 LOGICAL :: res
649 
650 res = this%proj == that%proj .AND. &
651  this%grid == that%grid
652 
653 END FUNCTION grid_eq
654 
655 
656 ELEMENTAL FUNCTION griddim_eq(this, that) RESULT(res)
657 TYPE(griddim_def),INTENT(IN) :: this, that
658 LOGICAL :: res
659 
660 res = this%grid == that%grid .AND. &
661  this%dim == that%dim
662 
663 END FUNCTION griddim_eq
664 
665 
666 ELEMENTAL FUNCTION grid_ne(this, that) RESULT(res)
667 TYPE(grid_def),INTENT(IN) :: this, that
668 LOGICAL :: res
669 
670 res = .NOT.(this == that)
671 
672 END FUNCTION grid_ne
673 
674 
675 ELEMENTAL FUNCTION griddim_ne(this, that) RESULT(res)
676 TYPE(griddim_def),INTENT(IN) :: this, that
677 LOGICAL :: res
678 
679 res = .NOT.(this == that)
680 
681 END FUNCTION griddim_ne
682 
683 
689 SUBROUTINE griddim_import_grid_id(this, ingrid_id)
690 #ifdef HAVE_LIBGDAL
691 USE gdal
692 #endif
693 TYPE(griddim_def),INTENT(inout) :: this
694 TYPE(grid_id),INTENT(in) :: ingrid_id
695 
696 #ifdef HAVE_LIBGRIBAPI
697 INTEGER :: gaid
698 #endif
699 #ifdef HAVE_LIBGDAL
700 TYPE(gdalrasterbandh) :: gdalid
701 #endif
702 CALL init(this)
703 
704 #ifdef HAVE_LIBGRIBAPI
705 gaid = grid_id_get_gaid(ingrid_id)
706 IF (c_e(gaid)) CALL griddim_import_gribapi(this, gaid)
707 #endif
708 #ifdef HAVE_LIBGDAL
709 gdalid = grid_id_get_gdalid(ingrid_id)
710 IF (gdalassociated(gdalid)) CALL griddim_import_gdal(this, gdalid, &
711  grid_id_get_gdal_options(ingrid_id))
712 #endif
713 
714 END SUBROUTINE griddim_import_grid_id
715 
716 
721 SUBROUTINE griddim_export_grid_id(this, outgrid_id)
722 #ifdef HAVE_LIBGDAL
723 USE gdal
724 #endif
725 TYPE(griddim_def),INTENT(in) :: this
726 TYPE(grid_id),INTENT(inout) :: outgrid_id
727 
728 #ifdef HAVE_LIBGRIBAPI
729 INTEGER :: gaid
730 #endif
731 #ifdef HAVE_LIBGDAL
732 TYPE(gdalrasterbandh) :: gdalid
733 #endif
734 
735 #ifdef HAVE_LIBGRIBAPI
736 gaid = grid_id_get_gaid(outgrid_id)
737 IF (c_e(gaid)) CALL griddim_export_gribapi(this, gaid)
738 #endif
739 #ifdef HAVE_LIBGDAL
740 gdalid = grid_id_get_gdalid(outgrid_id)
741 !IF (gdalassociated(gdalid)
742 ! export for gdal not implemented, log?
743 #endif
745 END SUBROUTINE griddim_export_grid_id
746 
747 
748 #ifdef HAVE_LIBGRIBAPI
749 ! grib_api driver
750 SUBROUTINE griddim_import_gribapi(this, gaid)
751 USE grib_api
752 TYPE(griddim_def),INTENT(inout) :: this ! griddim object
753 INTEGER, INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
754 
755 DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, x1, y1
756 INTEGER :: EditionNumber, iScansNegatively, jScansPositively, zone, datum, &
757  reflon, ierr
758 
759 ! Generic keys
760 CALL grib_get(gaid, 'typeOfGrid', this%grid%proj%proj_type)
761 #ifdef DEBUG
762 call l4f_category_log(this%category,l4f_debug, &
763  "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type))
764 #endif
765 CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
766 
767 ! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
768 CALL grib_get(gaid, 'Ni', this%dim%nx)
769 CALL grib_get(gaid, 'Nj', this%dim%ny)
770 CALL griddim_import_ellipsoid(this, gaid) ! placed here, not valid for utm datum /= 1
771 
772 CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
773 CALL grib_get(gaid,'jScansPositively',jscanspositively)
774 CALL grib_get(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
775 
776 ! Keys for rotated grids (checked through missing values)
777 CALL grib_get_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
778  this%grid%proj%rotated%longitude_south_pole)
779 CALL grib_get_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
780  this%grid%proj%rotated%latitude_south_pole)
781 CALL grib_get_dmiss(gaid,'angleOfRotationInDegrees', &
782  this%grid%proj%rotated%angle_rotation)
783 
784 ! Keys for stretched grids (checked through missing values)
785 ! units must be verified, still experimental in grib_api
786 ! # TODO: Is it a float? Is it signed?
787 IF (editionnumber == 1) THEN
788  CALL grib_get_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
789  this%grid%proj%stretched%longitude_stretch_pole)
790  CALL grib_get_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
791  this%grid%proj%stretched%latitude_stretch_pole)
792  CALL grib_get_dmiss(gaid,'stretchingFactor', &
793  this%grid%proj%stretched%stretch_factor)
794 ELSE IF (editionnumber == 2) THEN
795  CALL grib_get_dmiss(gaid,'longitudeOfThePoleOfStretching', &
796  this%grid%proj%stretched%longitude_stretch_pole)
797  CALL grib_get_dmiss(gaid,'latitudeOfThePoleOfStretching', &
798  this%grid%proj%stretched%latitude_stretch_pole)
799  CALL grib_get_dmiss(gaid,'stretchingFactor', &
800  this%grid%proj%stretched%stretch_factor)
801  IF (c_e(this%grid%proj%stretched%stretch_factor)) &
802  this%grid%proj%stretched%stretch_factor = this%grid%proj%stretched%stretch_factor*1.0d-6
803 ENDIF
804 
805 ! Projection-dependent keys
806 SELECT CASE (this%grid%proj%proj_type)
807 
808 ! Keys for spherical coordinate systems
809 CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
810 
811  CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
812  CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
813  CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
814  CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
815 
816 ! longitudes are sometimes wrongly coded even in grib2 and even by the
817 ! Metoffice!
818 ! longitudeOfFirstGridPointInDegrees = 354.911;
819 ! longitudeOfLastGridPointInDegrees = 363.311;
820  CALL long_reset_0_360(lofirst)
821  CALL long_reset_0_360(lolast)
822 
823  IF (iscansnegatively == 0) THEN
824  this%grid%grid%xmin = lofirst
825  this%grid%grid%xmax = lolast
826  ELSE
827  this%grid%grid%xmax = lofirst
828  this%grid%grid%xmin = lolast
829  ENDIF
830  IF (jscanspositively == 0) THEN
831  this%grid%grid%ymax = lafirst
832  this%grid%grid%ymin = lalast
833  ELSE
834  this%grid%grid%ymin = lafirst
835  this%grid%grid%ymax = lalast
836  ENDIF
837 
838 ! reset longitudes in order to have a Cartesian plane
839  IF (this%grid%grid%xmax < this%grid%grid%xmin) &
840  this%grid%grid%xmin = this%grid%grid%xmin - 360.d0
842 ! compute dx and dy (should we get them from grib?)
843  CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
844 
845 ! Keys for polar projections
846 CASE ('polar_stereographic', 'lambert', 'albers')
847 
848  CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
849  CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
850 ! latin1/latin2 may be missing (e.g. stereographic)
851  CALL grib_get_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
852  CALL grib_get_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
853 #ifdef DEBUG
854  CALL l4f_category_log(this%category,l4f_debug, &
855  "griddim_import_gribapi, latin1/2 "// &
856  trim(to_char(this%grid%proj%polar%latin1))//" "// &
857  trim(to_char(this%grid%proj%polar%latin2)))
858 #endif
859 ! projection center flag, aka hemisphere
860  CALL grib_get(gaid,'projectionCenterFlag',&
861  this%grid%proj%polar%projection_center_flag, ierr)
862  IF (ierr /= grib_success) THEN ! try center/centre
863  CALL grib_get(gaid,'projectionCentreFlag',&
864  this%grid%proj%polar%projection_center_flag)
865  ENDIF
866 
867  IF (iand(this%grid%proj%polar%projection_center_flag,64) == 1) THEN
868  CALL l4f_category_log(this%category,l4f_error, &
869  "griddim_import_gribapi, bi-polar projections not supported")
870  CALL raise_error()
871  ENDIF
872 ! line of view, aka central meridian
873  CALL grib_get(gaid,'LoVInDegrees',this%grid%proj%lov)
874 #ifdef DEBUG
875  CALL l4f_category_log(this%category,l4f_debug, &
876  "griddim_import_gribapi, central meridian "//trim(to_char(this%grid%proj%lov)))
877 #endif
878 
879 ! latitude at which dx and dy are valid
880  IF (editionnumber == 1) THEN
881 ! ECMWF (gribex/grib_api) says: Grid lengths are in metres, at the
882 ! 60-degree parallel nearest to the pole on the projection plane.
883 ! IF (IAND(this%projection_center_flag, 128) == 0) THEN
884 ! this%grid%proj%polar%lad = 60.D0
885 ! ELSE
886 ! this%grid%proj%polar%lad = -60.D0
887 ! ENDIF
888 ! WMO says: Grid lengths are in units of metres, at the secant cone
889 ! intersection parallel nearest to the pole on the projection plane.
890  this%grid%proj%polar%lad = this%grid%proj%polar%latin1
891  ELSE IF (editionnumber == 2) THEN
892  CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
893  ENDIF
894 #ifdef DEBUG
895  CALL l4f_category_log(this%category,l4f_debug, &
896  "griddim_import_gribapi, lad "//trim(to_char(this%grid%proj%polar%lad)))
897 #endif
898 
899 ! compute projected extremes from lon and lat of first point
900  CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
901  CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
902  CALL long_reset_0_360(lofirst)
903  CALL long_reset_to_cart_closest(this%grid%proj%lov, lofirst)
904 #ifdef DEBUG
905  CALL l4f_category_log(this%category,l4f_debug, &
906  "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
907  CALL l4f_category_log(this%category,l4f_debug, &
908  "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
909 #endif
910 
911 
912  CALL proj(this, lofirst, lafirst, x1, y1)
913  IF (iscansnegatively == 0) THEN
914  this%grid%grid%xmin = x1
915  this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
916  ELSE
917  this%grid%grid%xmax = x1
918  this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
919  ENDIF
920  IF (jscanspositively == 0) THEN
921  this%grid%grid%ymax = y1
922  this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
923  ELSE
924  this%grid%grid%ymin = y1
925  this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
926  ENDIF
927 ! keep these values for personal pleasure
928  this%grid%proj%polar%lon1 = lofirst
929  this%grid%proj%polar%lat1 = lafirst
930 
931 CASE ('UTM')
932 
933  CALL grib_get(gaid,'zone',zone)
935  CALL grib_get(gaid,'datum',datum)
936  IF (datum == 0) THEN
937  CALL grib_get(gaid,'referenceLongitude',reflon)
938  CALL grib_get(gaid,'falseEasting',this%grid%proj%xoff)
939  CALL grib_get(gaid,'falseNorthing',this%grid%proj%yoff)
940  CALL set_val(this%grid%proj, zone=zone, lov=reflon/1.0d6)
941  ELSE
942  CALL l4f_category_log(this%category,l4f_error,'only datum 0 supported')
943  CALL raise_fatal_error()
944  ENDIF
945 
946  CALL grib_get(gaid,'eastingOfFirstGridPoint',lofirst)
947  CALL grib_get(gaid,'eastingOfLastGridPoint',lolast)
948  CALL grib_get(gaid,'northingOfFirstGridPoint',lafirst)
949  CALL grib_get(gaid,'northingOfLastGridPoint',lalast)
950 
951  IF (iscansnegatively == 0) THEN
952  this%grid%grid%xmin = lofirst
953  this%grid%grid%xmax = lolast
954  ELSE
955  this%grid%grid%xmax = lofirst
956  this%grid%grid%xmin = lolast
957  ENDIF
958  IF (jscanspositively == 0) THEN
959  this%grid%grid%ymax = lafirst
960  this%grid%grid%ymin = lalast
961  ELSE
962  this%grid%grid%ymin = lafirst
963  this%grid%grid%ymax = lalast
964  ENDIF
965 
966 ! compute dx and dy (should we get them from grib?)
967  CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
968 
969 CASE default
970  CALL l4f_category_log(this%category,l4f_error, &
971  "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
972  CALL raise_error()
973 
974 END SELECT
975 
976 CONTAINS
977 ! utilities routines for grib_api, is there a better place?
978 SUBROUTINE grib_get_dmiss(gaid, key, value)
979 INTEGER,INTENT(in) :: gaid
980 CHARACTER(len=*),INTENT(in) :: key
981 DOUBLE PRECISION,INTENT(out) :: value
982 
983 INTEGER :: ierr
984 
985 CALL grib_get(gaid, key, value, ierr)
986 IF (ierr /= grib_success) value = dmiss
988 END SUBROUTINE grib_get_dmiss
989 
990 SUBROUTINE grib_get_imiss(gaid, key, value)
991 INTEGER,INTENT(in) :: gaid
992 CHARACTER(len=*),INTENT(in) :: key
993 INTEGER,INTENT(out) :: value
994 
995 INTEGER :: ierr
996 
997 CALL grib_get(gaid, key, value, ierr)
998 IF (ierr /= grib_success) value = imiss
999 
1000 END SUBROUTINE grib_get_imiss
1001 
1002 
1003 SUBROUTINE griddim_import_ellipsoid(this, gaid)
1004 TYPE(griddim_def),INTENT(inout) :: this
1005 INTEGER,INTENT(in) :: gaid
1006 
1007 INTEGER :: shapeofearth, iv, is
1008 DOUBLE PRECISION :: r1, r2
1009 
1010 IF (editionnumber == 2) THEN
1011  CALL grib_get(gaid, 'shapeOfTheEarth', shapeofearth)
1012  SELECT CASE(shapeofearth)
1013  CASE(0) ! spherical
1014  CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
1015  CASE(1) ! spherical generic
1016  CALL grib_get(gaid, 'scaleFactorOfRadiusOfSphericalEarth', is)
1017  CALL grib_get(gaid, 'scaledValueOfRadiusOfSphericalEarth', iv)
1018  r1 = dble(iv) / 10**is
1019  CALL set_val(this, ellips_smaj_axis=r1, ellips_flatt=0.0d0)
1020  CASE(2) ! iau65
1021  CALL set_val(this, ellips_smaj_axis=6378160.0d0, ellips_flatt=1.0d0/297.0d0)
1022  CASE(3,7) ! ellipsoidal generic
1023  CALL grib_get(gaid, 'scaleFactorOfEarthMajorAxis', is)
1024  CALL grib_get(gaid, 'scaledValueOfEarthMajorAxis', iv)
1025  r1 = dble(iv) / 10**is
1026  CALL grib_get(gaid, 'scaleFactorOfEarthMinorAxis', is)
1027  CALL grib_get(gaid, 'scaledValueOfEarthMinorAxis', iv)
1028  r2 = dble(iv) / 10**is
1029  IF (shapeofearth == 3) THEN ! km->m
1030  r1 = r1*1000.0d0
1031  r2 = r2*1000.0d0
1032  ENDIF
1033  IF (abs(r1) < 1.0d-6) THEN ! suspicious data read from grib
1034  CALL l4f_category_log(this%category,l4f_warn,'zero Earth major axis '// &
1035  'read from grib, going on with spherical Earth but the results may be wrong')
1036  CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
1037  ELSE
1038  CALL set_val(this, ellips_smaj_axis=r1, ellips_flatt=(r1-r2)/r1)
1039  ENDIF
1040  CASE(4) ! iag-grs80
1041  CALL set_val(this, ellips_type=ellips_grs80)
1042  CASE(5) ! wgs84
1043  CALL set_val(this, ellips_type=ellips_wgs84)
1044  CASE(6) ! spherical
1045  CALL set_val(this, ellips_smaj_axis=6371229.0d0, ellips_flatt=0.0d0)
1046 ! CASE(7) ! google earth-like?
1047  CASE default
1048  CALL l4f_category_log(this%category,l4f_error,'shapeOfTheEarth '// &
1049  t2c(shapeofearth)//' not supported in grib2')
1050  CALL raise_fatal_error()
1051 
1052  END SELECT
1053 
1054 ELSE
1055 
1056  CALL grib_get(gaid, 'earthIsOblate', shapeofearth)
1057  IF (shapeofearth == 0) THEN ! spherical
1058  CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
1059  ELSE ! iau65
1060  CALL set_val(this, ellips_smaj_axis=6378160.0d0, ellips_flatt=1.0d0/297.0d0)
1061  ENDIF
1062 
1063 ENDIF
1064 
1065 END SUBROUTINE griddim_import_ellipsoid
1066 
1067 
1068 END SUBROUTINE griddim_import_gribapi
1069 
1070 
1071 ! grib_api driver
1072 SUBROUTINE griddim_export_gribapi(this, gaid)
1073 USE grib_api
1074 TYPE(griddim_def),INTENT(in) :: this ! griddim object
1075 INTEGER, INTENT(inout) :: gaid ! grib_api id of the grib loaded in memory to export
1076 
1077 INTEGER :: EditionNumber, iScansNegatively, jScansPositively, nv, pvl, zone, ierr
1078 DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, reflon
1079 DOUBLE PRECISION :: sdx, sdy, ratio, tol
1080 
1081 
1082 ! Generic keys
1083 CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
1084 CALL grib_set(gaid,'typeOfGrid' ,this%grid%proj%proj_type)
1085 #ifdef DEBUG
1086 CALL l4f_category_log(this%category,l4f_debug, &
1087  "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
1088 #endif
1089 
1090 ! Edition dependent setup
1091 IF (editionnumber == 1) THEN
1092  ratio = 1.d3
1093 ELSE IF (editionnumber == 2) THEN
1094  ratio = 1.d6
1095 ELSE
1096  ratio = 0.0d0 ! signal error?!
1097 ENDIF
1098 
1099 ! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
1100 CALL griddim_export_ellipsoid(this, gaid)
1101 CALL grib_set(gaid, 'Ni', this%dim%nx)
1102 CALL grib_set(gaid, 'Nj', this%dim%ny)
1103 CALL grib_set(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
1104 
1105 CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
1106 CALL grib_get(gaid,'jScansPositively',jscanspositively)
1107 
1108 ! Keys for rotated grids (checked through missing values and/or error code)
1109 !SELECT CASE (this%grid%proj%proj_type)
1110 !CASE ('rotated_ll', 'stretched_rotated_ll', 'polar_stereographic', 'lambert', 'albers')
1111 CALL grib_set_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
1112  this%grid%proj%rotated%longitude_south_pole, 0.0d0)
1113 CALL grib_set_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
1114  this%grid%proj%rotated%latitude_south_pole, -90.0d0)
1115 IF (editionnumber == 1) THEN
1116  CALL grib_set_dmiss(gaid,'angleOfRotationInDegrees', &
1117  this%grid%proj%rotated%angle_rotation, 0.0d0)
1118 ELSE IF (editionnumber == 2)THEN
1119  CALL grib_set_dmiss(gaid,'angleOfRotationOfProjectionInDegrees', &
1120  this%grid%proj%rotated%angle_rotation, 0.0d0)
1121 ENDIF
1122 
1123 ! Keys for stretched grids (checked through missing values and/or error code)
1124 ! units must be verified, still experimental in grib_api
1125 ! # TODO: Is it a float? Is it signed?
1126 IF (editionnumber == 1) THEN
1127  CALL grib_set_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
1128  this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
1129  CALL grib_set_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
1130  this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
1131  CALL grib_set_dmiss(gaid,'stretchingFactor', &
1132  this%grid%proj%stretched%stretch_factor, 1.0d0)
1133 ELSE IF (editionnumber == 2) THEN
1134  CALL grib_set_dmiss(gaid,'longitudeOfThePoleOfStretching', &
1135  this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
1136  CALL grib_set_dmiss(gaid,'latitudeOfThePoleOfStretching', &
1137  this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
1138  CALL grib_set_dmiss(gaid,'stretchingFactor', &
1139  this%grid%proj%stretched%stretch_factor, 1.0d6, 1.0d6)
1140 ENDIF
1141 
1142 ! Projection-dependent keys
1143 SELECT CASE (this%grid%proj%proj_type)
1144 
1145 ! Keys for sphaerical coordinate systems
1146 CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
1147 
1148  IF (iscansnegatively == 0) THEN
1149  lofirst = this%grid%grid%xmin
1150  lolast = this%grid%grid%xmax
1151  ELSE
1152  lofirst = this%grid%grid%xmax
1153  lolast = this%grid%grid%xmin
1154  ENDIF
1155  IF (jscanspositively == 0) THEN
1156  lafirst = this%grid%grid%ymax
1157  lalast = this%grid%grid%ymin
1158  ELSE
1159  lafirst = this%grid%grid%ymin
1160  lalast = this%grid%grid%ymax
1161  ENDIF
1162 
1163 ! reset lon in standard grib 2 definition [0,360]
1164  IF (editionnumber == 1) THEN
1165  CALL long_reset_m180_360(lofirst)
1166  CALL long_reset_m180_360(lolast)
1167  ELSE IF (editionnumber == 2) THEN
1168  CALL long_reset_0_360(lofirst)
1169  CALL long_reset_0_360(lolast)
1170  ENDIF
1171 
1172  CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
1173  CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
1174  CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
1175  CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
1176 
1177 ! test relative coordinate truncation error with respect to tol
1178 ! tol should be tuned
1179  sdx = this%grid%grid%dx*ratio
1180  sdy = this%grid%grid%dy*ratio
1181  tol = 1.0d0/ratio
1182 
1183  IF (abs(nint(sdx)/sdx - 1.0d0) > tol .OR. abs(nint(sdy)/sdy - 1.0d0) > tol) THEN
1184  CALL l4f_category_log(this%category,l4f_info, &
1185  "griddim_export_gribapi, increments not given: inaccurate!")
1186 #ifdef DEBUG
1187  CALL l4f_category_log(this%category,l4f_debug, &
1188  "griddim_export_gribapi, dlon relative error: "//&
1189  trim(to_char(abs(nint(sdx)/sdx - 1.0d0)))//">"//trim(to_char(tol)))
1190  CALL l4f_category_log(this%category,l4f_debug, &
1191  "griddim_export_gribapi, dlat relative error: "//&
1192  trim(to_char(abs(nint(sdy)/sdy - 1.0d0)))//">"//trim(to_char(tol)))
1193 #endif
1194  CALL grib_set_missing(gaid,'Di')
1195  CALL grib_set_missing(gaid,'Dj')
1196  CALL grib_set(gaid,'ijDirectionIncrementGiven',0)
1197  ELSE
1198 #ifdef DEBUG
1199  CALL l4f_category_log(this%category,l4f_debug, &
1200  "griddim_export_gribapi, setting increments: "// &
1201  trim(to_char(this%grid%grid%dx))//' '//trim(to_char(this%grid%grid%dy)))
1202 #endif
1203  CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
1204  CALL grib_set(gaid,'iDirectionIncrement',nint(this%grid%grid%dx*ratio))
1205  CALL grib_set(gaid,'jDirectionIncrement',nint(this%grid%grid%dy*ratio))
1206 ! this does not work in grib_set
1207 ! CALL grib_set(gaid,'iDirectionIncrementInDegrees',this%grid%grid%dx)
1208 ! CALL grib_set(gaid,'jDirectionIncrementInDegrees',this%grid%grid%dy)
1209  ENDIF
1210 
1211 ! Keys for polar projections
1212 CASE ('polar_stereographic', 'lambert', 'albers')
1213 ! increments are required
1214  CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
1215  CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
1216  CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
1217 ! latin1/latin2 may be missing (e.g. stereographic)
1218  CALL grib_set_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
1219  CALL grib_set_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
1220 ! projection center flag, aka hemisphere
1221  CALL grib_set(gaid,'projectionCenterFlag',&
1222  this%grid%proj%polar%projection_center_flag, ierr)
1223  IF (ierr /= grib_success) THEN ! try center/centre
1224  CALL grib_set(gaid,'projectionCentreFlag',&
1225  this%grid%proj%polar%projection_center_flag)
1226  ENDIF
1227 
1228 
1229 ! line of view, aka central meridian
1230  CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
1231 ! latitude at which dx and dy are valid
1232  IF (editionnumber == 2) THEN
1233  CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
1234  ENDIF
1235 
1236 ! compute lon and lat of first point from projected extremes
1237  IF (iscansnegatively == 0) THEN
1238  IF (jscanspositively == 0) THEN
1239  CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymax, lofirst, lafirst)
1240  ELSE
1241  CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymin, lofirst, lafirst)
1242  ENDIF
1243  ELSE
1244  IF (jscanspositively == 0) THEN
1245  CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymax, lofirst, lafirst)
1246  ELSE
1247  CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymin, lofirst, lafirst)
1248  ENDIF
1249  ENDIF
1250 ! use the values kept for personal pleasure ?
1251 ! loFirst = this%grid%proj%polar%lon1
1252 ! laFirst = this%grid%proj%polar%lat1
1253 ! reset lon in standard grib 2 definition [0,360]
1254  IF (editionnumber == 1) THEN
1255  CALL long_reset_m180_360(lofirst)
1256  ELSE IF (editionnumber == 2) THEN
1257  CALL long_reset_0_360(lofirst)
1258  ENDIF
1259  CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
1260  CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
1261 
1262 CASE ('UTM')
1263 
1264  CALL grib_set(gaid,'datum',0)
1265  CALL get_val(this, zone=zone, lov=reflon)
1266  CALL grib_set(gaid,'referenceLongitude',nint(reflon/1.0d6))
1267  CALL grib_set(gaid,'falseEasting',this%grid%proj%xoff)
1268  CALL grib_set(gaid,'falseNorthing',this%grid%proj%yoff)
1269 
1270  CALL grib_set(gaid,'iDirectionIncrement',this%grid%grid%dx)
1271  CALL grib_set(gaid,'jDirectionIncrement',this%grid%grid%dy)
1272 
1273 !res/scann ??
1274 
1275  CALL grib_set(gaid,'zone',zone)
1276 
1277  IF (iscansnegatively == 0) THEN
1278  lofirst = this%grid%grid%xmin
1279  lolast = this%grid%grid%xmax
1280  ELSE
1281  lofirst = this%grid%grid%xmax
1282  lolast = this%grid%grid%xmin
1283  ENDIF
1284  IF (jscanspositively == 0) THEN
1285  lafirst = this%grid%grid%ymax
1286  lalast = this%grid%grid%ymin
1287  ELSE
1288  lafirst = this%grid%grid%ymin
1289  lalast = this%grid%grid%ymax
1290  ENDIF
1291 
1292  CALL grib_set(gaid,'eastingOfFirstGridPoint',lofirst)
1293  CALL grib_set(gaid,'eastingOfLastGridPoint',lolast)
1294  CALL grib_set(gaid,'northingOfFirstGridPoint',lafirst)
1295  CALL grib_set(gaid,'northingOfLastGridPoint',lalast)
1296 
1297 CASE default
1298  CALL l4f_category_log(this%category,l4f_error, &
1299  "griddim_export_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
1300  CALL raise_error()
1301 
1302 END SELECT
1303 
1304 ! hack for position of vertical coordinate parameters
1305 ! buggy in grib_api
1306 IF (editionnumber == 1) THEN
1307 ! CALL grib_get(gaid,"PVPresent",pvp) ! alias, probably useless
1308  CALL grib_get(gaid,"NV",nv)
1309 #ifdef DEBUG
1310  CALL l4f_category_log(this%category,l4f_debug,"griddim_export_gribapi, coding "// &
1311  trim(to_char(nv))//" vertical coordinate parameters")
1312 #endif
1313 
1314  IF (nv == 0) THEN
1315  pvl = 255
1316  ELSE
1317  SELECT CASE (this%grid%proj%proj_type)
1318  CASE ('regular_ll') ! check whether "29-32 Set to zero (reserved)" are required
1319  pvl = 33
1320  CASE ('polar_stereographic')
1321  pvl = 33
1322  CASE ('rotated_ll', 'stretched_ll', 'lambert', 'albers')
1323  pvl = 43
1324  CASE ('stretched_rotated_ll')
1325  pvl = 43
1326  CASE DEFAULT
1327  pvl = 43 !?
1328  END SELECT
1329  ENDIF
1330 
1331  CALL grib_set(gaid,"pvlLocation",pvl)
1332 #ifdef DEBUG
1333  CALL l4f_category_log(this%category,l4f_debug,"griddim_export_gribapi, coding "// &
1334  trim(to_char(pvl))//" as vertical coordinate parameter location")
1335 #endif
1336 
1337 ENDIF
1338 
1339 
1340 CONTAINS
1341 ! utilities routines for grib_api, is there a better place?
1342 SUBROUTINE grib_set_dmiss(gaid, key, val, default, factor)
1343 INTEGER,INTENT(in) :: gaid
1344 CHARACTER(len=*),INTENT(in) :: key
1345 DOUBLE PRECISION,INTENT(in) :: val
1346 DOUBLE PRECISION,INTENT(in),OPTIONAL :: default
1347 DOUBLE PRECISION,INTENT(in),OPTIONAL :: factor
1348 
1349 INTEGER :: ierr
1350 
1351 IF (c_e(val)) THEN
1352  IF (PRESENT(factor)) THEN
1353  CALL grib_set(gaid, key, val*factor, ierr)
1354  ELSE
1355  CALL grib_set(gaid, key, val, ierr)
1356  ENDIF
1357 ELSE IF (PRESENT(default)) THEN
1358  CALL grib_set(gaid, key, default, ierr)
1359 ENDIF
1360 
1361 END SUBROUTINE grib_set_dmiss
1362 
1363 SUBROUTINE grib_set_imiss(gaid, key, value, default)
1364 INTEGER,INTENT(in) :: gaid
1365 CHARACTER(len=*),INTENT(in) :: key
1366 INTEGER,INTENT(in) :: value
1367 INTEGER,INTENT(in),OPTIONAL :: default
1368 
1369 INTEGER :: ierr
1370 
1371 IF (c_e(value)) THEN
1372  CALL grib_set(gaid, key, value, ierr)
1373 ELSE IF (PRESENT(default)) THEN
1374  CALL grib_set(gaid, key, default, ierr)
1375 ENDIF
1376 
1377 END SUBROUTINE grib_set_imiss
1378 
1379 SUBROUTINE griddim_export_ellipsoid(this, gaid)
1380 TYPE(griddim_def),INTENT(in) :: this
1381 INTEGER,INTENT(in) :: gaid
1382 
1383 INTEGER :: ellips_type
1384 DOUBLE PRECISION :: r1, r2, f
1385 
1386 CALL get_val(this, ellips_smaj_axis=r1, ellips_flatt=f, ellips_type=ellips_type)
1387 
1388 IF (editionnumber == 2) THEN
1389 
1390  CALL grib_set_missing(gaid, 'scaleFactorOfRadiusOfSphericalEarth')
1391  CALL grib_set_missing(gaid, 'scaledValueOfRadiusOfSphericalEarth')
1392  CALL grib_set_missing(gaid, 'scaleFactorOfEarthMajorAxis')
1393  CALL grib_set_missing(gaid, 'scaledValueOfEarthMajorAxis')
1394  CALL grib_set_missing(gaid, 'scaleFactorOfEarthMinorAxis')
1395  CALL grib_set_missing(gaid, 'scaledValueOfEarthMinorAxis')
1396 
1397  SELECT CASE(ellips_type)
1398  CASE(ellips_grs80) ! iag-grs80
1399  CALL grib_set(gaid, 'shapeOfTheEarth', 4)
1400  CASE(ellips_wgs84) ! wgs84
1401  CALL grib_set(gaid, 'shapeOfTheEarth', 5)
1402  CASE default
1403  IF (f == 0.0d0) THEN ! spherical Earth
1404  IF (r1 == 6367470.0d0) THEN ! spherical
1405  CALL grib_set(gaid, 'shapeOfTheEarth', 0)
1406  ELSE IF (r1 == 6371229.0d0) THEN ! spherical
1407  CALL grib_set(gaid, 'shapeOfTheEarth', 6)
1408  ELSE ! spherical generic
1409  CALL grib_set(gaid, 'shapeOfTheEarth', 1)
1410  CALL grib_set(gaid, 'scaleFactorOfRadiusOfSphericalEarth', 2)
1411  CALL grib_set(gaid, 'scaledValueOfRadiusOfSphericalEarth', int(r1*100.0d0))
1412  ENDIF
1413  ELSE ! ellipsoidal
1414  IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
1415  CALL grib_set(gaid, 'shapeOfTheEarth', 2)
1416  ELSE ! ellipsoidal generic
1417  CALL grib_set(gaid, 'shapeOfTheEarth', 3)
1418  r2 = r1*(1.0d0 - f)
1419  CALL grib_set(gaid, 'scaleFactorOfEarthMajorAxis', 5)
1420  CALL grib_set(gaid, 'scaledValueOfEarthMajorAxis', &
1421  int(r1*100.0d0))
1422  CALL grib_set(gaid, 'scaleFactorOfEarthMinorAxis', 5)
1423  CALL grib_set(gaid, 'scaledValueOfEarthMinorAxis', &
1424  int(r2*100.0d0))
1425  ENDIF
1426  ENDIF
1427  END SELECT
1428 
1429 ELSE
1430 
1431  IF (r1 == 6367470.0d0 .AND. f == 0.0d0) THEN ! spherical
1432  CALL grib_set(gaid, 'earthIsOblate', 0)
1433  ELSE IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
1434  CALL grib_set(gaid, 'earthIsOblate', 1)
1435  ELSE IF (f == 0.0d0) THEN ! generic spherical
1436  CALL grib_set(gaid, 'earthIsOblate', 0)
1437  CALL l4f_category_log(this%category,l4f_warn,'desired spherical Earth radius &
1438  &not supported in grib 1, coding with default radius of 6367470 m')
1439  ELSE ! generic ellipsoidal
1440  CALL grib_set(gaid, 'earthIsOblate', 1)
1441  CALL l4f_category_log(this%category,l4f_warn,'desired Earth ellipsoid &
1442  &not supported in grib 1, coding with default iau65 ellipsoid')
1443  ENDIF
1444 
1445 ENDIF
1446 
1447 END SUBROUTINE griddim_export_ellipsoid
1448 
1449 END SUBROUTINE griddim_export_gribapi
1450 #endif
1451 
1452 
1453 #ifdef HAVE_LIBGDAL
1454 ! gdal driver
1455 SUBROUTINE griddim_import_gdal(this, gdalid, gdal_options)
1456 USE gdal
1457 TYPE(griddim_def),INTENT(inout) :: this ! griddim object
1458 TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal rasterband pointer
1459 TYPE(gdal_file_id_options),INTENT(in) :: gdal_options
1460 
1461 TYPE(gdaldataseth) :: hds
1462 REAL(kind=c_double) :: geotrans(6) !, invgeotrans(6), x1, y1, x2, y2
1463 INTEGER(kind=c_int) :: offsetx, offsety
1464 INTEGER :: ier
1465 
1466 hds = gdalgetbanddataset(gdalid) ! go back to dataset
1467 ier = gdalgetgeotransform(hds, geotrans)
1468 
1469 IF (ier /= 0) THEN
1470  CALL l4f_category_log(this%category, l4f_error, &
1471  'griddim_import_gdal: error in accessing gdal dataset')
1472  CALL raise_error()
1473  RETURN
1474 ENDIF
1475 IF (geotrans(3) /= 0.0_c_double .OR. geotrans(5) /= 0.0_c_double) THEN ! dataset has not a diagonal transformation
1476  CALL l4f_category_log(this%category, l4f_error, &
1477  'griddim_import_gdal: dataset has a non-diagonal transformation matrix, unsupported')
1478  CALL raise_error()
1479  RETURN
1480 ENDIF
1481 
1482 CALL gdaldatasetbbsize_f(hds, gdal_options%xmin, gdal_options%ymin, &
1483  gdal_options%xmax, gdal_options%ymax, &
1484  this%dim%nx, this%dim%ny, offsetx, offsety, &
1485  this%grid%grid%xmin, this%grid%grid%ymin, this%grid%grid%xmax, this%grid%grid%ymax)
1486 
1487 IF (this%dim%nx == 0 .OR. this%dim%ny == 0) THEN
1488  CALL l4f_category_log(this%category, l4f_warn, &
1489  'griddim_import_gdal: requested bounding box '//t2c(gdal_options%xmin)//','// &
1490  t2c(gdal_options%ymin)//','//t2c(gdal_options%xmax)//','//&
1491  t2c(gdal_options%ymax))
1492  CALL l4f_category_log(this%category, l4f_warn, &
1493  'determines an empty dataset '//t2c(this%dim%nx)//'x'//t2c(this%dim%ny))
1494 ENDIF
1495 
1496 ! get grid corners
1497 !CALL gdalapplygeotransform(geotrans, 0.5_c_double, 0.5_c_double, x1, y1)
1498 !CALL gdalapplygeotransform(geotrans, gdalgetrasterbandxsize(gdalid)-0.5_c_double, &
1499 ! gdalgetrasterbandysize(gdalid)-0.5_c_double, x2, y2)
1500 
1501 !IF (geotrans(3) == 0.0_c_double .AND. geotrans(5) == 0.0_c_double) THEN ! transformation is diagonal, no transposing
1502 ! this%dim%nx = gdalgetrasterbandxsize(gdalid)
1503 ! this%dim%ny = gdalgetrasterbandysize(gdalid)
1504 ! this%grid%grid%xmin = MIN(x1, x2)
1505 ! this%grid%grid%xmax = MAX(x1, x2)
1506 ! this%grid%grid%ymin = MIN(y1, y2)
1507 ! this%grid%grid%ymax = MAX(y1, y2)
1508 !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
1509 !
1510 ! this%dim%nx = gdalgetrasterbandysize(gdalid)
1511 ! this%dim%ny = gdalgetrasterbandxsize(gdalid)
1512 ! this%grid%grid%xmin = MIN(y1, y2)
1513 ! this%grid%grid%xmax = MAX(y1, y2)
1514 ! this%grid%grid%ymin = MIN(x1, x2)
1515 ! this%grid%grid%ymax = MAX(x1, x2)
1516 !
1517 !ELSE ! transformation is a rotation, not supported
1518 !ENDIF
1519 
1520 this%grid%proj%proj_type = 'regular_ll' ! forced, only one supported (check?)
1521 
1522 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
1523 this%grid%grid%component_flag = 0
1524 
1525 END SUBROUTINE griddim_import_gdal
1526 #endif
1527 
1528 
1530 SUBROUTINE griddim_display(this)
1531 TYPE(griddim_def),INTENT(in) :: this
1532 
1533 print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
1534 
1535 CALL display(this%grid%proj)
1536 CALL display(this%grid%grid)
1537 CALL display(this%dim)
1538 
1539 print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
1540 
1541 END SUBROUTINE griddim_display
1542 
1543 
1544 #define VOL7D_POLY_TYPE TYPE(grid_def)
1545 #define VOL7D_POLY_TYPES _grid
1546 #include "array_utilities_inc.F90"
1547 #undef VOL7D_POLY_TYPE
1548 #undef VOL7D_POLY_TYPES
1549 
1550 #define VOL7D_POLY_TYPE TYPE(griddim_def)
1551 #define VOL7D_POLY_TYPES _griddim
1552 #include "array_utilities_inc.F90"
1553 #undef VOL7D_POLY_TYPE
1554 #undef VOL7D_POLY_TYPES
1555 
1556 
1568 SUBROUTINE griddim_wind_unrot(this, rot_mat)
1570 TYPE(griddim_def), INTENT(in) :: this
1571 DOUBLE PRECISION, POINTER :: rot_mat(:,:,:)
1572 
1573 DOUBLE PRECISION :: dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j
1574 DOUBLE PRECISION :: lat_factor
1575 INTEGER :: i, j, ip1, im1, jp1, jm1;
1576 
1577 IF (this%dim%nx <= 0 .OR. this%dim%ny <= 0 .OR. &
1578  .NOT. ASSOCIATED(this%dim%lon) .OR. .NOT. ASSOCIATED(this%dim%lat)) THEN
1579  CALL l4f_category_log(this%category, l4f_error, 'rotate_uv must be called after correct init of griddim object')
1580  NULLIFY(rot_mat)
1581  RETURN
1582 ENDIF
1583 
1584 ALLOCATE(rot_mat(this%dim%nx, this%dim%ny, 4))
1585 
1586 DO j = 1, this%dim%ny
1587  jp1 = min(j+1, this%dim%ny)
1588  jm1 = max(j-1, 1)
1589  DO i = 1, this%dim%nx
1590  ip1 = min(i+1, this%dim%nx)
1591  im1 = max(i-1, 1)
1592 
1593  dlat_i = this%dim%lat(ip1,j) - this%dim%lat(im1,j)
1594  dlat_j = this%dim%lat(i,jp1) - this%dim%lat(i,jm1)
1595 
1596  dlon_i = this%dim%lon(ip1,j) - this%dim%lon(im1,j)
1597 ! if ( dlon_i > pi ) dlon_i -= 2*pi;
1598 ! if ( dlon_i < (-pi) ) dlon_i += 2*pi;
1599  dlon_j = this%dim%lon(i,jp1) - this%dim%lon(i,jm1)
1600 ! if ( dlon_j > pi ) dlon_j -= 2*pi;
1601 ! if ( dlon_j < (-pi) ) dlon_j += 2*pi;
1602 
1603 ! check whether this is really necessary !!!!
1604  lat_factor = cos(degrad*this%dim%lat(i,j))
1605  dlon_i = dlon_i * lat_factor
1606  dlon_j = dlon_j * lat_factor
1607 
1608  dist_i = sqrt(dlon_i*dlon_i + dlat_i*dlat_i)
1609  dist_j = sqrt(dlon_j*dlon_j + dlat_j*dlat_j)
1610 
1611  IF (dist_i > 0.d0) THEN
1612  rot_mat(i,j,1) = dlon_i/dist_i
1613  rot_mat(i,j,2) = dlat_i/dist_i
1614  ELSE
1615  rot_mat(i,j,1) = 0.0d0
1616  rot_mat(i,j,2) = 0.0d0
1617  ENDIF
1618  IF (dist_j > 0.d0) THEN
1619  rot_mat(i,j,3) = dlon_j/dist_j
1620  rot_mat(i,j,4) = dlat_j/dist_j
1621  ELSE
1622  rot_mat(i,j,3) = 0.0d0
1623  rot_mat(i,j,4) = 0.0d0
1624  ENDIF
1625 
1626  ENDDO
1627 ENDDO
1628 
1629 END SUBROUTINE griddim_wind_unrot
1630 
1631 
1632 ! compute zoom indices from geographical zoom coordinates
1633 SUBROUTINE griddim_zoom_coord(this, ilon, ilat, flon, flat, ix, iy, fx, fy)
1634 TYPE(griddim_def),INTENT(in) :: this
1635 DOUBLE PRECISION,INTENT(in) :: ilon, ilat, flon, flat
1636 INTEGER,INTENT(out) :: ix, iy, fx, fy
1637 
1638 DOUBLE PRECISION :: ix1, iy1, fx1, fy1
1639 
1640 ! compute projected coordinates of vertices of desired lonlat rectangle
1641 CALL proj(this, ilon, ilat, ix1, iy1)
1642 CALL proj(this, flon, flat, fx1, fy1)
1643 
1644 CALL griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
1645 
1646 END SUBROUTINE griddim_zoom_coord
1647 
1648 
1649 ! compute zoom indices from projected zoom coordinates
1650 SUBROUTINE griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
1651 TYPE(griddim_def),INTENT(in) :: this
1652 DOUBLE PRECISION,INTENT(in) :: ix1, iy1, fx1, fy1
1653 INTEGER,INTENT(out) :: ix, iy, fx, fy
1654 
1655 INTEGER :: lix, liy, lfx, lfy
1656 
1657 ! compute projected indices
1658 lix = nint((ix1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
1659 liy = nint((iy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
1660 lfx = nint((fx1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
1661 lfy = nint((fy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
1662 ! swap projected indices, in case grid is "strongly rotated"
1663 ix = min(lix, lfx)
1664 fx = max(lix, lfx)
1665 iy = min(liy, lfy)
1666 fy = max(liy, lfy)
1667 
1668 END SUBROUTINE griddim_zoom_projcoord
1669 
1670 
1674 SUBROUTINE long_reset_0_360(lon)
1675 DOUBLE PRECISION,INTENT(inout) :: lon
1676 
1677 IF (.NOT.c_e(lon)) RETURN
1678 DO WHILE(lon < 0.0d0)
1679  lon = lon + 360.0d0
1680 END DO
1681 DO WHILE(lon >= 360.0d0)
1682  lon = lon - 360.0d0
1683 END DO
1684 
1685 END SUBROUTINE long_reset_0_360
1686 
1687 
1691 SUBROUTINE long_reset_m180_360(lon)
1692 DOUBLE PRECISION,INTENT(inout) :: lon
1693 
1694 IF (.NOT.c_e(lon)) RETURN
1695 DO WHILE(lon < -180.0d0)
1696  lon = lon + 360.0d0
1697 END DO
1698 DO WHILE(lon >= 360.0d0)
1699  lon = lon - 360.0d0
1700 END DO
1701 
1702 END SUBROUTINE long_reset_m180_360
1703 
1704 
1708 !SUBROUTINE long_reset_m90_270(lon)
1709 !DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
1710 !
1711 !IF (.NOT.c_e(lon)) RETURN
1712 !DO WHILE(lon < -90.0D0)
1713 ! lon = lon + 360.0D0
1714 !END DO
1715 !DO WHILE(lon >= 270.0D0)
1716 ! lon = lon - 360.0D0
1717 !END DO
1718 !
1719 !END SUBROUTINE long_reset_m90_270
1720 
1721 
1725 SUBROUTINE long_reset_m180_180(lon)
1726 DOUBLE PRECISION,INTENT(inout) :: lon
1727 
1728 IF (.NOT.c_e(lon)) RETURN
1729 DO WHILE(lon < -180.0d0)
1730  lon = lon + 360.0d0
1731 END DO
1732 DO WHILE(lon >= 180.0d0)
1733  lon = lon - 360.0d0
1734 END DO
1735 
1736 END SUBROUTINE long_reset_m180_180
1737 
1738 
1739 SUBROUTINE long_reset_to_cart_closest(lon, lonref)
1740 DOUBLE PRECISION,INTENT(inout) :: lon
1741 DOUBLE PRECISION,INTENT(in) :: lonref
1742 
1743 IF (.NOT.c_e(lon) .OR. .NOT.c_e(lonref)) RETURN
1744 IF (abs(lon-lonref) < 180.0d0) RETURN ! nothing to do
1745 lon = lon - nint((lon-lonref)/360.0d0)*360.0d0 ! se non e` vera e` ben trovata
1746 
1747 END SUBROUTINE long_reset_to_cart_closest
1748 
1749 
1750 END MODULE grid_class
1751 
Export griddim object to grid_id.
Definition: grid_class.F90:369
Method for setting the contents of the object.
Definition: grid_class.F90:349
This object, mainly for internal use, describes a grid on a geographical projection, except the grid dimensions.
Definition: grid_class.F90:282
Compute forward coordinate transformation from geographical system to projected system.
Definition: grid_class.F90:333
Constructors of the corresponding objects.
Definition: grid_class.F90:317
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Print a brief description on stdout.
Definition: grid_class.F90:374
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:262
Write the object on a formatted or unformatted file.
Definition: grid_class.F90:354
Compute backward coordinate transformation from projected system to geographical system.
Definition: grid_class.F90:339
This object completely describes a grid on a geographic projection.
Definition: grid_class.F90:295
Method for returning the contents of the object.
Definition: grid_class.F90:344
Index method.
Gestione degli errori.
Definitions of constants and functions for working with missing values.
Costanti fisiche (DOUBLEPRECISION).
Definition: phys_const.f90:67
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:327
Read the object from a formatted or unformatted file.
Definition: grid_class.F90:359
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:322
Import griddim object from grid_id.
Definition: grid_class.F90:364

Generated with Doxygen.