libsim Versione 7.2.4
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
49MODULE grid_class
50use geo_proj_class
52use grid_rect_class
58implicit none
59
60
61character (len=255),parameter:: subcategory="grid_class"
62
63
69type grid_def
70 private
71 type(geo_proj) :: proj
72 type(grid_rect) :: grid
73 integer :: category = 0
74end type grid_def
75
76
82type griddim_def
83 type(grid_def) :: grid
84 type(grid_dim) :: dim
85 integer :: category = 0
86end type griddim_def
87
88
92INTERFACE OPERATOR (==)
93 MODULE PROCEDURE grid_eq, griddim_eq
94END INTERFACE
95
99INTERFACE OPERATOR (/=)
100 MODULE PROCEDURE grid_ne, griddim_ne
101END INTERFACE
102
104INTERFACE init
105 MODULE PROCEDURE griddim_init
106END INTERFACE
107
109INTERFACE delete
110 MODULE PROCEDURE griddim_delete
111END INTERFACE
112
114INTERFACE copy
115 MODULE PROCEDURE griddim_copy
116END INTERFACE
117
120INTERFACE proj
121 MODULE PROCEDURE griddim_coord_proj!, griddim_proj
122END INTERFACE
123
126INTERFACE unproj
127 MODULE PROCEDURE griddim_coord_unproj, griddim_unproj
128END INTERFACE
129
131INTERFACE get_val
132 MODULE PROCEDURE griddim_get_val
133END INTERFACE
134
136INTERFACE set_val
137 MODULE PROCEDURE griddim_set_val
138END INTERFACE
139
141INTERFACE write_unit
142 MODULE PROCEDURE griddim_write_unit
143END INTERFACE
144
146INTERFACE read_unit
147 MODULE PROCEDURE griddim_read_unit
148END INTERFACE
149
151INTERFACE import
152 MODULE PROCEDURE griddim_import_grid_id
153END INTERFACE
154
156INTERFACE export
157 MODULE PROCEDURE griddim_export_grid_id
158END INTERFACE
159
161INTERFACE display
162 MODULE PROCEDURE griddim_display
163END 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
177INTERFACE wind_unrot
178 MODULE PROCEDURE griddim_wind_unrot
179END INTERFACE
180
181
182PRIVATE
183
184PUBLIC proj, unproj, griddim_unproj, griddim_gen_coord, &
185 griddim_zoom_coord, griddim_zoom_projcoord, &
186 griddim_setsteps, griddim_def, grid_def, grid_dim
187PUBLIC init, delete, copy
189PUBLIC OPERATOR(==),OPERATOR(/=)
190PUBLIC count_distinct, pack_distinct, count_and_pack_distinct, &
191 map_distinct, map_inv_distinct,index
192PUBLIC wind_unrot, import, export
193PUBLIC griddim_central_lon, griddim_set_central_lon
194CONTAINS
195
197SUBROUTINE 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)
205TYPE(griddim_def),INTENT(inout) :: this
206INTEGER,INTENT(in),OPTIONAL :: nx
207INTEGER,INTENT(in),OPTIONAL :: ny
208DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin
209DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmax
210DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymin
211DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymax
212DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx
213DOUBLE PRECISION,INTENT(in),OPTIONAL :: dy
216INTEGER,INTENT(in),OPTIONAL :: component_flag
217CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
218DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
219INTEGER,INTENT(in),OPTIONAL :: zone
220DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
221DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
222DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
223DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
224DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
225DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
226DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
227DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
228DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
229DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
230DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
231INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
232DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
233DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
234INTEGER,INTENT(in),OPTIONAL :: ellips_type
235CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
236
237CHARACTER(len=512) :: a_name
238
239IF (PRESENT(categoryappend)) THEN
240 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
241 trim(categoryappend))
242ELSE
243 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
244ENDIF
245this%category=l4f_category_get(a_name)
246
247! geographical projection
248this%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
257this%grid%grid = grid_rect_new( &
258 xmin, xmax, ymin, ymax, dx, dy, component_flag)
259! grid size
260this%dim = grid_dim_new(nx, ny)
261
262#ifdef DEBUG
263call l4f_category_log(this%category,l4f_debug,"init gtype: "//this%grid%proj%proj_type )
264#endif
265
266END SUBROUTINE griddim_init
267
268
270SUBROUTINE griddim_delete(this)
271TYPE(griddim_def),INTENT(inout) :: this
273CALL delete(this%dim)
274CALL delete(this%grid%proj)
275CALL delete(this%grid%grid)
276
277call l4f_category_delete(this%category)
278
279END SUBROUTINE griddim_delete
281
283SUBROUTINE griddim_copy(this, that, categoryappend)
284TYPE(griddim_def),INTENT(in) :: this
285TYPE(griddim_def),INTENT(out) :: that
286CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
288CHARACTER(len=512) :: a_name
289
290CALL init(that)
291
292CALL copy(this%grid%proj, that%grid%proj)
293CALL copy(this%grid%grid, that%grid%grid)
294CALL copy(this%dim, that%dim)
295
296! new category
297IF (PRESENT(categoryappend)) THEN
298 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
299 trim(categoryappend))
300ELSE
301 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
302ENDIF
303that%category=l4f_category_get(a_name)
304
305END SUBROUTINE griddim_copy
306
307
310ELEMENTAL SUBROUTINE griddim_coord_proj(this, lon, lat, x, y)
311TYPE(griddim_def),INTENT(in) :: this
313DOUBLE PRECISION,INTENT(in) :: lon, lat
315DOUBLE PRECISION,INTENT(out) :: x, y
316
317CALL proj(this%grid%proj, lon, lat, x, y)
318
319END SUBROUTINE griddim_coord_proj
320
321
324ELEMENTAL SUBROUTINE griddim_coord_unproj(this, x, y, lon, lat)
325TYPE(griddim_def),INTENT(in) :: this
327DOUBLE PRECISION,INTENT(in) :: x, y
329DOUBLE PRECISION,INTENT(out) :: lon, lat
330
331CALL unproj(this%grid%proj, x, y, lon, lat)
332
333END 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)
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
358SUBROUTINE griddim_unproj(this)
359TYPE(griddim_def),INTENT(inout) :: this
360
361IF (.NOT.c_e(this%dim%nx) .OR. .NOT.c_e(this%dim%ny)) RETURN
362CALL alloc(this%dim)
363CALL griddim_unproj_internal(this)
364
365END SUBROUTINE griddim_unproj
366
367! internal subroutine needed for allocating automatic arrays
368SUBROUTINE griddim_unproj_internal(this)
369TYPE(griddim_def),INTENT(inout) ::this ! definition of projection and grid
370
371DOUBLE PRECISION :: x(this%dim%nx,this%dim%ny), y(this%dim%nx,this%dim%ny)
372
373CALL grid_rect_coordinates(this%grid%grid, x, y)
374CALL unproj(this, x, y, this%dim%lon, this%dim%lat)
375
376END SUBROUTINE griddim_unproj_internal
377
378
380SUBROUTINE 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)
387TYPE(griddim_def),INTENT(in) :: this
388INTEGER,INTENT(out),OPTIONAL :: nx
389INTEGER,INTENT(out),OPTIONAL :: ny
391DOUBLE PRECISION,INTENT(out),OPTIONAL :: xmin, xmax, ymin, ymax
393DOUBLE PRECISION,INTENT(out),OPTIONAL :: dx, dy
396INTEGER,INTENT(out),OPTIONAL :: component_flag
397TYPE(geo_proj),INTENT(out),OPTIONAL :: proj
398CHARACTER(len=*),INTENT(out),OPTIONAL :: proj_type
399DOUBLE PRECISION,INTENT(out),OPTIONAL :: lov
400INTEGER,INTENT(out),OPTIONAL :: zone
401DOUBLE PRECISION,INTENT(out),OPTIONAL :: xoff
402DOUBLE PRECISION,INTENT(out),OPTIONAL :: yoff
403DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_south_pole
404DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_south_pole
405DOUBLE PRECISION,INTENT(out),OPTIONAL :: angle_rotation
406DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_stretch_pole
407DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_stretch_pole
408DOUBLE PRECISION,INTENT(out),OPTIONAL :: stretch_factor
409DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin1
410DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin2
411DOUBLE PRECISION,INTENT(out),OPTIONAL :: lad
412INTEGER,INTENT(out),OPTIONAL :: projection_center_flag
413DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_smaj_axis
414DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_flatt
415INTEGER,INTENT(out),OPTIONAL :: ellips_type
416
417IF (PRESENT(nx)) nx = this%dim%nx
418IF (PRESENT(ny)) ny = this%dim%ny
419
420IF (PRESENT(proj)) proj = this%grid%proj
421
422CALL 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
433CALL get_val(this%grid%grid, &
434 xmin, xmax, ymin, ymax, dx, dy, component_flag)
435
436END SUBROUTINE griddim_get_val
437
438
440SUBROUTINE 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)
447TYPE(griddim_def),INTENT(inout) :: this
448INTEGER,INTENT(in),OPTIONAL :: nx
449INTEGER,INTENT(in),OPTIONAL :: ny
451DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin, xmax, ymin, ymax
453DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx, dy
456INTEGER,INTENT(in),OPTIONAL :: component_flag
457CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
458DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
459INTEGER,INTENT(in),OPTIONAL :: zone
460DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
461DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
462DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
463DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
464DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
465DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
466DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
467DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
468DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
469DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
470DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
471INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
472DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
473DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
474INTEGER,INTENT(in),OPTIONAL :: ellips_type
475
476IF (PRESENT(nx)) this%dim%nx = nx
477IF (PRESENT(ny)) this%dim%ny = ny
478
479CALL 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
489CALL set_val(this%grid%grid, &
490 xmin, xmax, ymin, ymax, dx, dy, component_flag)
491
492END SUBROUTINE griddim_set_val
493
494
499SUBROUTINE griddim_read_unit(this, unit)
500TYPE(griddim_def),INTENT(out) :: this
501INTEGER, INTENT(in) :: unit
502
503
504CALL read_unit(this%dim, unit)
505CALL read_unit(this%grid%proj, unit)
506CALL read_unit(this%grid%grid, unit)
507
508END SUBROUTINE griddim_read_unit
509
510
515SUBROUTINE griddim_write_unit(this, unit)
516TYPE(griddim_def),INTENT(in) :: this
517INTEGER, INTENT(in) :: unit
518
519
520CALL write_unit(this%dim, unit)
521CALL write_unit(this%grid%proj, unit)
522CALL write_unit(this%grid%grid, unit)
523
524END SUBROUTINE griddim_write_unit
525
526
530FUNCTION griddim_central_lon(this) RESULT(lon)
531TYPE(griddim_def),INTENT(inout) :: this
532
533DOUBLE PRECISION :: lon
534
535CALL griddim_pistola_central_lon(this, lon)
536
537END FUNCTION griddim_central_lon
538
539
543SUBROUTINE griddim_set_central_lon(this, lonref)
544TYPE(griddim_def),INTENT(inout) :: this
545DOUBLE PRECISION,INTENT(in) :: lonref
546
547DOUBLE PRECISION :: lon
548
549CALL griddim_pistola_central_lon(this, lon, lonref)
550
551END SUBROUTINE griddim_set_central_lon
552
553
554! internal subroutine for performing tasks common to the prevous two
555SUBROUTINE griddim_pistola_central_lon(this, lon, lonref)
556TYPE(griddim_def),INTENT(inout) :: this ! grid descriptor
557DOUBLE PRECISION,INTENT(inout) :: lon ! central longitude
558DOUBLE PRECISION,INTENT(in),OPTIONAL :: lonref ! reference longitude
559
560INTEGER :: unit
561DOUBLE PRECISION :: lonsp, latsp, londelta, lov, lonrot
562CHARACTER(len=80) :: ptype
563
564lon = dmiss
565CALL get_val(this%grid%proj, unit=unit)
566IF (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
573ELSE 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
612ENDIF
613
614END SUBROUTINE griddim_pistola_central_lon
615
616
620SUBROUTINE griddim_gen_coord(this, x, y)
621TYPE(griddim_def),INTENT(in) :: this
622DOUBLE PRECISION,INTENT(out) :: x(:,:)
623DOUBLE PRECISION,INTENT(out) :: y(:,:)
624
625
626CALL grid_rect_coordinates(this%grid%grid, x, y)
627
628END SUBROUTINE griddim_gen_coord
629
632SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
633TYPE(griddim_def), INTENT(in) :: this
634INTEGER,INTENT(in) :: nx
635INTEGER,INTENT(in) :: ny
636DOUBLE PRECISION,INTENT(out) :: dx
637DOUBLE PRECISION,INTENT(out) :: dy
638
639CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
640
641END SUBROUTINE griddim_steps
642
643
645SUBROUTINE griddim_setsteps(this)
646TYPE(griddim_def), INTENT(inout) :: this
647
648CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
649
650END SUBROUTINE griddim_setsteps
651
653! TODO
654! bisogna sviluppare gli altri operatori
655ELEMENTAL FUNCTION grid_eq(this, that) RESULT(res)
656TYPE(grid_def),INTENT(IN) :: this, that
657LOGICAL :: res
658
659res = this%proj == that%proj .AND. &
660 this%grid == that%grid
661
662END FUNCTION grid_eq
663
664
665ELEMENTAL FUNCTION griddim_eq(this, that) RESULT(res)
666TYPE(griddim_def),INTENT(IN) :: this, that
667LOGICAL :: res
668
669res = this%grid == that%grid .AND. &
670 this%dim == that%dim
671
672END FUNCTION griddim_eq
673
674
675ELEMENTAL FUNCTION grid_ne(this, that) RESULT(res)
676TYPE(grid_def),INTENT(IN) :: this, that
677LOGICAL :: res
678
679res = .NOT.(this == that)
680
681END FUNCTION grid_ne
682
683
684ELEMENTAL FUNCTION griddim_ne(this, that) RESULT(res)
685TYPE(griddim_def),INTENT(IN) :: this, that
686LOGICAL :: res
687
688res = .NOT.(this == that)
689
690END FUNCTION griddim_ne
691
692
698SUBROUTINE griddim_import_grid_id(this, ingrid_id)
699#ifdef HAVE_LIBGDAL
700USE gdal
701#endif
702TYPE(griddim_def),INTENT(inout) :: this
703TYPE(grid_id),INTENT(in) :: ingrid_id
704
705#ifdef HAVE_LIBGRIBAPI
706INTEGER :: gaid
707#endif
708#ifdef HAVE_LIBGDAL
709TYPE(gdalrasterbandh) :: gdalid
710#endif
711CALL init(this)
713#ifdef HAVE_LIBGRIBAPI
714gaid = grid_id_get_gaid(ingrid_id)
715IF (c_e(gaid)) CALL griddim_import_gribapi(this, gaid)
716#endif
717#ifdef HAVE_LIBGDAL
718gdalid = grid_id_get_gdalid(ingrid_id)
719IF (gdalassociated(gdalid)) CALL griddim_import_gdal(this, gdalid, &
720 grid_id_get_gdal_options(ingrid_id))
721#endif
722
723END SUBROUTINE griddim_import_grid_id
724
725
730SUBROUTINE griddim_export_grid_id(this, outgrid_id)
731#ifdef HAVE_LIBGDAL
732USE gdal
733#endif
734TYPE(griddim_def),INTENT(in) :: this
735TYPE(grid_id),INTENT(inout) :: outgrid_id
736
737#ifdef HAVE_LIBGRIBAPI
738INTEGER :: gaid
739#endif
740#ifdef HAVE_LIBGDAL
741TYPE(gdalrasterbandh) :: gdalid
742#endif
743
744#ifdef HAVE_LIBGRIBAPI
745gaid = grid_id_get_gaid(outgrid_id)
746IF (c_e(gaid)) CALL griddim_export_gribapi(this, gaid)
747#endif
748#ifdef HAVE_LIBGDAL
749gdalid = grid_id_get_gdalid(outgrid_id)
750!IF (gdalassociated(gdalid)
751! export for gdal not implemented, log?
752#endif
753
754END SUBROUTINE griddim_export_grid_id
755
756
757#ifdef HAVE_LIBGRIBAPI
758! grib_api driver
759SUBROUTINE griddim_import_gribapi(this, gaid)
760USE grib_api
761TYPE(griddim_def),INTENT(inout) :: this ! griddim object
762INTEGER, INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
763
764DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, x1, y1, orient
765INTEGER :: EditionNumber, iScansNegatively, jScansPositively, zone, datum, &
766 reflon, ierr
767
768! Generic keys
769CALL grib_get(gaid, 'typeOfGrid', this%grid%proj%proj_type)
770#ifdef DEBUG
771call l4f_category_log(this%category,l4f_debug, &
772 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type))
773#endif
774CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
775
776IF (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
782ENDIF
783
784! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
785CALL grib_get(gaid, 'Ni', this%dim%nx)
786CALL grib_get(gaid, 'Nj', this%dim%ny)
787CALL griddim_import_ellipsoid(this, gaid) ! placed here, not valid for utm datum /= 1
788
789CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
790CALL grib_get(gaid,'jScansPositively',jscanspositively)
791CALL grib_get(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
792
793! Keys for rotated grids (checked through missing values)
794CALL grib_get_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
795 this%grid%proj%rotated%longitude_south_pole)
796CALL grib_get_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
797 this%grid%proj%rotated%latitude_south_pole)
798CALL 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?
804IF (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)
811ELSE 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
820ENDIF
821
822! Projection-dependent keys
823SELECT CASE (this%grid%proj%proj_type)
824
825! Keys for spherical coordinate systems
826CASE ('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
863CASE ('polar_stereographic', 'lambert', 'albers')
864
865 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
866 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
867! latin1/latin2 may be missing (e.g. stereographic)
868 CALL grib_get_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
869 CALL grib_get_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
870#ifdef DEBUG
871 CALL l4f_category_log(this%category,l4f_debug, &
872 "griddim_import_gribapi, latin1/2 "// &
873 trim(to_char(this%grid%proj%polar%latin1))//" "// &
874 trim(to_char(this%grid%proj%polar%latin2)))
875#endif
876! projection center flag, aka hemisphere
877 CALL grib_get(gaid,'projectionCenterFlag',&
878 this%grid%proj%polar%projection_center_flag, ierr)
879 IF (ierr /= grib_success) THEN ! try center/centre
880 CALL grib_get(gaid,'projectionCentreFlag',&
881 this%grid%proj%polar%projection_center_flag)
882 ENDIF
883
884 IF (iand(this%grid%proj%polar%projection_center_flag,64) == 1) THEN
885 CALL l4f_category_log(this%category,l4f_error, &
886 "griddim_import_gribapi, bi-polar projections not supported")
887 CALL raise_error()
888 ENDIF
889! line of view, aka central meridian
890 CALL grib_get(gaid,'LoVInDegrees',this%grid%proj%lov)
891#ifdef DEBUG
892 CALL l4f_category_log(this%category,l4f_debug, &
893 "griddim_import_gribapi, central meridian "//trim(to_char(this%grid%proj%lov)))
894#endif
895
896! latitude at which dx and dy are valid
897 IF (editionnumber == 1) THEN
898! ECMWF (gribex/grib_api) says: Grid lengths are in metres, at the
899! 60-degree parallel nearest to the pole on the projection plane.
900! IF (IAND(this%projection_center_flag, 128) == 0) THEN
901! this%grid%proj%polar%lad = 60.D0
902! ELSE
903! this%grid%proj%polar%lad = -60.D0
904! ENDIF
905! WMO says: Grid lengths are in units of metres, at the secant cone
906! intersection parallel nearest to the pole on the projection plane.
907 this%grid%proj%polar%lad = this%grid%proj%polar%latin1
908 ELSE IF (editionnumber == 2) THEN
909 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
910 ENDIF
911#ifdef DEBUG
912 CALL l4f_category_log(this%category,l4f_debug, &
913 "griddim_import_gribapi, lad "//trim(to_char(this%grid%proj%polar%lad)))
914#endif
915
916! compute projected extremes from lon and lat of first point
917 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
918 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
919 CALL long_reset_0_360(lofirst)
920 CALL long_reset_to_cart_closest(this%grid%proj%lov, lofirst)
921#ifdef DEBUG
922 CALL l4f_category_log(this%category,l4f_debug, &
923 "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
924 CALL l4f_category_log(this%category,l4f_debug, &
925 "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
926#endif
927
928 CALL proj(this, lofirst, lafirst, x1, y1)
929 IF (iscansnegatively == 0) THEN
930 this%grid%grid%xmin = x1
931 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
932 ELSE
933 this%grid%grid%xmax = x1
934 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
935 ENDIF
936 IF (jscanspositively == 0) THEN
937 this%grid%grid%ymax = y1
938 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
939 ELSE
940 this%grid%grid%ymin = y1
941 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
942 ENDIF
943! keep these values for personal pleasure
944 this%grid%proj%polar%lon1 = lofirst
945 this%grid%proj%polar%lat1 = lafirst
946
947! Keys for equatorial projections
948CASE ('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
1000CASE ('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
1038CASE 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
1043END SELECT
1044
1045CONTAINS
1046! utilities routines for grib_api, is there a better place?
1047SUBROUTINE grib_get_dmiss(gaid, key, value)
1048INTEGER,INTENT(in) :: gaid
1049CHARACTER(len=*),INTENT(in) :: key
1050DOUBLE PRECISION,INTENT(out) :: value
1051
1052INTEGER :: ierr
1053
1054CALL grib_get(gaid, key, value, ierr)
1055IF (ierr /= grib_success) value = dmiss
1056
1057END SUBROUTINE grib_get_dmiss
1058
1059SUBROUTINE grib_get_imiss(gaid, key, value)
1060INTEGER,INTENT(in) :: gaid
1061CHARACTER(len=*),INTENT(in) :: key
1062INTEGER,INTENT(out) :: value
1063
1064INTEGER :: ierr
1065
1066CALL grib_get(gaid, key, value, ierr)
1067IF (ierr /= grib_success) value = imiss
1068
1069END SUBROUTINE grib_get_imiss
1070
1071
1072SUBROUTINE griddim_import_ellipsoid(this, gaid)
1073TYPE(griddim_def),INTENT(inout) :: this
1074INTEGER,INTENT(in) :: gaid
1075
1076INTEGER :: shapeofearth, iv, is
1077DOUBLE PRECISION :: r1, r2
1078
1079IF (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
1123ELSE
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
1132ENDIF
1133
1134END SUBROUTINE griddim_import_ellipsoid
1135
1136
1137END SUBROUTINE griddim_import_gribapi
1138
1139
1140! grib_api driver
1141SUBROUTINE griddim_export_gribapi(this, gaid)
1142USE grib_api
1143TYPE(griddim_def),INTENT(in) :: this ! griddim object
1144INTEGER, INTENT(inout) :: gaid ! grib_api id of the grib loaded in memory to export
1145
1146INTEGER :: EditionNumber, iScansNegatively, jScansPositively, nv, pvl, zone, ierr
1147DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, reflon
1148DOUBLE PRECISION :: sdx, sdy, ratio, tol, ex, ey
1149
1150
1151! Generic keys
1152CALL grib_get(gaid,'GRIBEditionNumber', editionnumber)
1153! the following required since eccodes
1154IF (editionnumber == 2) CALL grib_set(gaid,'shapeOfTheEarth', 6)
1155CALL grib_set(gaid,'typeOfGrid', this%grid%proj%proj_type)
1156#ifdef DEBUG
1157CALL l4f_category_log(this%category,l4f_debug, &
1158 "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
1159#endif
1160
1161IF (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
1166ENDIF
1167
1168
1169! Edition dependent setup
1170IF (editionnumber == 1) THEN
1171 ratio = 1.d3
1172ELSE IF (editionnumber == 2) THEN
1173 ratio = 1.d6
1174ELSE
1175 ratio = 0.0d0 ! signal error?!
1176ENDIF
1177
1178! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
1179CALL griddim_export_ellipsoid(this, gaid)
1180CALL grib_set(gaid, 'Ni', this%dim%nx)
1181CALL grib_set(gaid, 'Nj', this%dim%ny)
1182CALL grib_set(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
1183
1184CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
1185CALL 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')
1190CALL grib_set_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
1191 this%grid%proj%rotated%longitude_south_pole, 0.0d0)
1192CALL grib_set_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
1193 this%grid%proj%rotated%latitude_south_pole, -90.0d0)
1194IF (editionnumber == 1) THEN
1195 CALL grib_set_dmiss(gaid,'angleOfRotationInDegrees', &
1196 this%grid%proj%rotated%angle_rotation, 0.0d0)
1197ELSE IF (editionnumber == 2)THEN
1198 CALL grib_set_dmiss(gaid,'angleOfRotationOfProjectionInDegrees', &
1199 this%grid%proj%rotated%angle_rotation, 0.0d0)
1200ENDIF
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?
1205IF (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)
1212ELSE 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)
1219ENDIF
1220
1221! Projection-dependent keys
1222SELECT CASE (this%grid%proj%proj_type)
1223
1224! Keys for sphaerical coordinate systems
1225CASE ('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
1306CASE ('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
1343CASE ('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
1369CASE ('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
1404CASE 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
1409END SELECT
1410
1411! hack for position of vertical coordinate parameters
1412! buggy in grib_api
1413IF (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
1444ENDIF
1445
1446
1447CONTAINS
1448! utilities routines for grib_api, is there a better place?
1449SUBROUTINE grib_set_dmiss(gaid, key, val, default, factor)
1450INTEGER,INTENT(in) :: gaid
1451CHARACTER(len=*),INTENT(in) :: key
1452DOUBLE PRECISION,INTENT(in) :: val
1453DOUBLE PRECISION,INTENT(in),OPTIONAL :: default
1454DOUBLE PRECISION,INTENT(in),OPTIONAL :: factor
1455
1456INTEGER :: ierr
1457
1458IF (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
1464ELSE IF (PRESENT(default)) THEN
1465 CALL grib_set(gaid, key, default, ierr)
1466ENDIF
1467
1468END SUBROUTINE grib_set_dmiss
1469
1470SUBROUTINE grib_set_imiss(gaid, key, value, default)
1471INTEGER,INTENT(in) :: gaid
1472CHARACTER(len=*),INTENT(in) :: key
1473INTEGER,INTENT(in) :: value
1474INTEGER,INTENT(in),OPTIONAL :: default
1475
1476INTEGER :: ierr
1477
1478IF (c_e(value)) THEN
1479 CALL grib_set(gaid, key, value, ierr)
1480ELSE IF (PRESENT(default)) THEN
1481 CALL grib_set(gaid, key, default, ierr)
1482ENDIF
1483
1484END SUBROUTINE grib_set_imiss
1485
1486SUBROUTINE griddim_export_ellipsoid(this, gaid)
1487TYPE(griddim_def),INTENT(in) :: this
1488INTEGER,INTENT(in) :: gaid
1489
1490INTEGER :: ellips_type, ierr
1491DOUBLE PRECISION :: r1, r2, f
1492
1493CALL get_val(this, ellips_smaj_axis=r1, ellips_flatt=f, ellips_type=ellips_type)
1494
1495IF (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
1537ELSE
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
1553ENDIF
1554
1555END SUBROUTINE griddim_export_ellipsoid
1556
1557SUBROUTINE get_unproj_first(this, iScansNegatively, jScansPositively, &
1558 loFirst, laFirst)
1559TYPE(griddim_def),INTENT(in) :: this ! griddim object
1560INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
1561DOUBLE PRECISION,INTENT(out) :: loFirst, laFirst
1562
1563! compute lon and lat of first point from projected extremes
1564IF (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
1570ELSE
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
1576ENDIF
1577! use the values kept for personal pleasure ?
1578! loFirst = this%grid%proj%polar%lon1
1579! laFirst = this%grid%proj%polar%lat1
1580END SUBROUTINE get_unproj_first
1581
1582SUBROUTINE get_unproj_last(this, iScansNegatively, jScansPositively, &
1583 loLast, laLast)
1584TYPE(griddim_def),INTENT(in) :: this ! griddim object
1585INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
1586DOUBLE PRECISION,INTENT(out) :: loLast, laLast
1587
1588! compute lon and lat of last point from projected extremes
1589IF (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
1595ELSE
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
1601ENDIF
1602! use the values kept for personal pleasure ?
1603! loLast = this%grid%proj%polar%lon?
1604! laLast = this%grid%proj%polar%lat?
1605END SUBROUTINE get_unproj_last
1606
1607END SUBROUTINE griddim_export_gribapi
1608#endif
1609
1610
1611#ifdef HAVE_LIBGDAL
1612! gdal driver
1613SUBROUTINE griddim_import_gdal(this, gdalid, gdal_options)
1614USE gdal
1615TYPE(griddim_def),INTENT(inout) :: this ! griddim object
1616TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal rasterband pointer
1617TYPE(gdal_file_id_options),INTENT(in) :: gdal_options
1618
1619TYPE(gdaldataseth) :: hds
1620REAL(kind=c_double) :: geotrans(6) !, invgeotrans(6), x1, y1, x2, y2
1621INTEGER(kind=c_int) :: offsetx, offsety
1622INTEGER :: ier
1623
1624hds = gdalgetbanddataset(gdalid) ! go back to dataset
1625ier = gdalgetgeotransform(hds, geotrans)
1626
1627IF (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
1632ENDIF
1633IF (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
1638ENDIF
1639
1640CALL 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
1645IF (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))
1652ENDIF
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
1678this%grid%proj%proj_type = 'regular_ll' ! forced, only one supported (check?)
1679
1680CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
1681this%grid%grid%component_flag = 0
1682
1683END SUBROUTINE griddim_import_gdal
1684#endif
1685
1686
1688SUBROUTINE griddim_display(this)
1689TYPE(griddim_def),INTENT(in) :: this
1690
1691print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
1692
1693CALL display(this%grid%proj)
1694CALL display(this%grid%grid)
1695CALL display(this%dim)
1696
1697print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
1698
1699END 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
1726SUBROUTINE griddim_wind_unrot(this, rot_mat)
1728TYPE(griddim_def), INTENT(in) :: this
1729DOUBLE PRECISION, POINTER :: rot_mat(:,:,:)
1730
1731DOUBLE PRECISION :: dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j
1732DOUBLE PRECISION :: lat_factor
1733INTEGER :: i, j, ip1, im1, jp1, jm1;
1734
1735IF (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
1740ENDIF
1741
1742ALLOCATE(rot_mat(this%dim%nx, this%dim%ny, 4))
1743
1744DO 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
1785ENDDO
1786
1787END SUBROUTINE griddim_wind_unrot
1788
1789
1790! compute zoom indices from geographical zoom coordinates
1791SUBROUTINE griddim_zoom_coord(this, ilon, ilat, flon, flat, ix, iy, fx, fy)
1792TYPE(griddim_def),INTENT(in) :: this
1793DOUBLE PRECISION,INTENT(in) :: ilon, ilat, flon, flat
1794INTEGER,INTENT(out) :: ix, iy, fx, fy
1795
1796DOUBLE PRECISION :: ix1, iy1, fx1, fy1
1797
1798! compute projected coordinates of vertices of desired lonlat rectangle
1799CALL proj(this, ilon, ilat, ix1, iy1)
1800CALL proj(this, flon, flat, fx1, fy1)
1801
1802CALL griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
1803
1804END SUBROUTINE griddim_zoom_coord
1805
1806
1807! compute zoom indices from projected zoom coordinates
1808SUBROUTINE griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
1809TYPE(griddim_def),INTENT(in) :: this
1810DOUBLE PRECISION,INTENT(in) :: ix1, iy1, fx1, fy1
1811INTEGER,INTENT(out) :: ix, iy, fx, fy
1812
1813INTEGER :: lix, liy, lfx, lfy
1814
1815! compute projected indices
1816lix = nint((ix1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
1817liy = nint((iy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
1818lfx = nint((fx1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
1819lfy = nint((fy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
1820! swap projected indices, in case grid is "strongly rotated"
1821ix = min(lix, lfx)
1822fx = max(lix, lfx)
1823iy = min(liy, lfy)
1824fy = max(liy, lfy)
1825
1826END SUBROUTINE griddim_zoom_projcoord
1827
1828
1832SUBROUTINE long_reset_0_360(lon)
1833DOUBLE PRECISION,INTENT(inout) :: lon
1834
1835IF (.NOT.c_e(lon)) RETURN
1836DO WHILE(lon < 0.0d0)
1837 lon = lon + 360.0d0
1838END DO
1839DO WHILE(lon >= 360.0d0)
1840 lon = lon - 360.0d0
1841END DO
1842
1843END SUBROUTINE long_reset_0_360
1844
1845
1849SUBROUTINE long_reset_m180_360(lon)
1850DOUBLE PRECISION,INTENT(inout) :: lon
1851
1852IF (.NOT.c_e(lon)) RETURN
1853DO WHILE(lon < -180.0d0)
1854 lon = lon + 360.0d0
1855END DO
1856DO WHILE(lon >= 360.0d0)
1857 lon = lon - 360.0d0
1858END DO
1859
1860END 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
1883SUBROUTINE long_reset_m180_180(lon)
1884DOUBLE PRECISION,INTENT(inout) :: lon
1885
1886IF (.NOT.c_e(lon)) RETURN
1887DO WHILE(lon < -180.0d0)
1888 lon = lon + 360.0d0
1889END DO
1890DO WHILE(lon >= 180.0d0)
1891 lon = lon - 360.0d0
1892END DO
1893
1894END SUBROUTINE long_reset_m180_180
1895
1896
1897SUBROUTINE long_reset_to_cart_closest(lon, lonref)
1898DOUBLE PRECISION,INTENT(inout) :: lon
1899DOUBLE PRECISION,INTENT(in) :: lonref
1900
1901IF (.NOT.c_e(lon) .OR. .NOT.c_e(lonref)) RETURN
1902IF (abs(lon-lonref) < 180.0d0) RETURN ! nothing to do
1903lon = lon - nint((lon-lonref)/360.0d0)*360.0d0 ! se non e` vera e` ben trovata
1904
1905END SUBROUTINE long_reset_to_cart_closest
1906
1907
1908END MODULE grid_class
1909
Method for testing the existence of the object.
Copy an object, creating a fully new instance.
Destructors of the corresponding objects.
Print a brief description on stdout.
Export griddim object to grid_id.
Method for returning the contents of the object.
Import griddim object from grid_id.
Constructors of the corresponding objects.
Compute forward coordinate transformation from geographical system to projected system.
Read the object from a formatted or unformatted file.
Method for setting the contents of the object.
Compute backward coordinate transformation from projected system to geographical system.
Write the object on a formatted or unformatted file.
Index method.
Emit log message for a category with specific priority.
Costanti fisiche (DOUBLEPRECISION).
Gestione degli errori.
Module for describing geographically referenced regular grids.
Module for defining the extension and coordinates of a rectangular georeferenced grid.
This module defines an abstract interface to different drivers for access to files containing gridded...
classe per la gestione del logging
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
This object, mainly for internal use, describes a grid on a geographical projection,...
This object completely describes a grid on a geographic projection.
Derived type describing the extension of a grid and the geographical coordinates of each point.
Derived type containing driver-specific options for gdal.
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...

Generated with Doxygen.