libsim Versione 7.2.6
grid_class.F90
1! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2! authors:
3! Davide Cesari <dcesari@arpa.emr.it>
4! Paolo Patruno <ppatruno@arpa.emr.it>
5
6! This program is free software; you can redistribute it and/or
7! modify it under the terms of the GNU General Public License as
8! published by the Free Software Foundation; either version 2 of
9! the License, or (at your option) any later version.
10
11! This program is distributed in the hope that it will be useful,
12! but WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14! GNU General Public License for more details.
15
16! You should have received a copy of the GNU General Public License
17! along with this program. If not, see <http://www.gnu.org/licenses/>.
18#include "config.h"
19!> Module for describing geographically referenced regular grids.
20!! This module defines classes and methods describing rectangular
21!! georeferenced grids in different geographical projections. The grid
22!! and projection definition can be specified explicitely by the
23!! caller or they can be entirely imported from a grib file (through
24!! grib_api) or from a file format supported by gdal (through gdal
25!! fortran interface).
26!!
27!! The projection is internally stored following the WMO grib
28!! conventions (gridType in grib_api). The projections currently
29!! supported or for which support is planned are:
30!! - For grib edition 1 and 2:
31!! - regular_ll (works)
32!! - mercator (to be done)
33!! - lambert (works, to be completed)
34!! - polar_stereographic (to be tested)
35!! - albers (to be done)
36!! - rotated_ll (works)
37!! - stretched_ll (to be completed and tested)
38!! - stretched_rotated_ll (to be completed and tested)
39!! - For grib edition 2 only:
40!! - UTM (ARPA-SIM extension)
41!! - equatorial_azimuthal_equidistant (to be done)
42!! - For gdal-supported formats:
43!! - regular_ll
44!!
45!!
46!! See the example program \include example_vg6d_1.f90
47!!
48!!\ingroup volgrid6d
49MODULE grid_class
50use geo_proj_class
52use grid_rect_class
58implicit none
59
60
61character (len=255),parameter:: subcategory="grid_class"
62
63
64!> This object, mainly for internal use, describes a grid on
65!! a geographical projection, except the grid dimensions.
66!! The object is opaque, thus all its members have to be set and
67!! accessed through the constructor and the ::get_val and ::set_val
68!! methods.
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
77!> This object completely describes a grid on a geographic projection.
78!! It is the main public object of this module. The grid definition \a
79!! grid is separated from the definition of the grid dimensions \a dim
80!! in order to make members of \a grid \a PRIVATE while maintaining
81!! free access to the members of \a dim.
82type griddim_def
83 type(grid_def) :: grid !< grid and projection definition
84 type(grid_dim) :: dim !< grid dimensions definition
85 integer :: category = 0 !< category for log4fortran
86end type griddim_def
87
88
89!> Logical equality operators for objects of the classes \a grid_def,
90!! and \a griddim_def. They are all defined as \c
91!! ELEMENTAL thus work also on arrays of any shape.
92INTERFACE OPERATOR (==)
93 MODULE PROCEDURE grid_eq, griddim_eq
94END INTERFACE
95
96!> Logical inequality operators for objects of the classes \a grid_def,
97!! and \a griddim_def. They are all defined as \c
98!! ELEMENTAL thus work also on arrays of any shape.
99INTERFACE OPERATOR (/=)
100 MODULE PROCEDURE grid_ne, griddim_ne
101END INTERFACE
102
103!> Constructors of the corresponding objects.
104INTERFACE init
105 MODULE PROCEDURE griddim_init
106END INTERFACE
107
108!> Destructors of the corresponding objects.
109INTERFACE delete
110 MODULE PROCEDURE griddim_delete
111END INTERFACE
112
113!> Copy an object, creating a fully new instance.
114INTERFACE copy
115 MODULE PROCEDURE griddim_copy
116END INTERFACE
117
118!> Compute forward coordinate transformation from geographical system to
119!! projected system.
120INTERFACE proj
121 MODULE PROCEDURE griddim_coord_proj!, griddim_proj
122END INTERFACE
123
124!> Compute backward coordinate transformation from projected system
125!! to geographical system.
126INTERFACE unproj
127 MODULE PROCEDURE griddim_coord_unproj, griddim_unproj
128END INTERFACE
129
130!> Method for returning the contents of the object.
131INTERFACE get_val
132 MODULE PROCEDURE griddim_get_val
133END INTERFACE
134
135!> Method for setting the contents of the object.
136INTERFACE set_val
137 MODULE PROCEDURE griddim_set_val
138END INTERFACE
139
140!> Write the object on a formatted or unformatted file.
141INTERFACE write_unit
142 MODULE PROCEDURE griddim_write_unit
143END INTERFACE
144
145!> Read the object from a formatted or unformatted file.
146INTERFACE read_unit
147 MODULE PROCEDURE griddim_read_unit
148END INTERFACE
149
150!> Import griddim object from grid_id.
151INTERFACE import
152 MODULE PROCEDURE griddim_import_grid_id
153END INTERFACE
154
155!> Export griddim object to grid_id.
156INTERFACE export
157 MODULE PROCEDURE griddim_export_grid_id
158END INTERFACE
159
160!> Print a brief description on stdout.
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
196!> Constructor for a \a griddim_def object.
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 !< object to be created
206INTEGER,INTENT(in),OPTIONAL :: nx !< number of points along the x axis
207INTEGER,INTENT(in),OPTIONAL :: ny !< number of points along the y axis
208DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin !< lower bound for x coordinate on grid in projection units (degrees or meters depending on the projection type)
209DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmax !< upper bound for x coordinate
210DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymin !< lower bound for y coordinate
211DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymax !< upper bound for y coordinate
212DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx !< grid step in x direction
213DOUBLE PRECISION,INTENT(in),OPTIONAL :: dy !< grid step in y direction
214!> Resolved u- and v- components of vector quantities relative to 0=the easterly and northerly directions
215!! 1=the defined grid in the direction of increasing x and y (or i and j) coordinates respectively (0=north, 128=south)
216INTEGER,INTENT(in),OPTIONAL :: component_flag
217CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type !< type of projection
218DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov !< line of view, also known as reference longitude or orientation of the grid (polar projections)
219INTEGER,INTENT(in),OPTIONAL :: zone !< Earth zone (mainly for UTM), sets lov to the correct zone central meridian
220DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff !< offset on x axis (false easting)
221DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff !< offset on y axis (false northing)
222DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole !< longitude of the southern pole of projection
223DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole !< latitude of the southern pole of projection
224DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation !< angle of rotation of projection
225DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole !< longitude of the pole of stretching
226DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole !< latitude of the pole of stretching
227DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor !< stretching factor
228DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1 !< first standard latitude from main pole (Lambert)
229DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2 !< second standard latitude from main pole (Lambert)
230DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad !< latitude at which dx and dy (in m) are specified (Lambert, grib2 only)
231INTEGER,INTENT(in),OPTIONAL :: projection_center_flag !< flag indicating which pole is represented
232DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis !< Earth semi-major axis
233DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt !< Earth flattening
234INTEGER,INTENT(in),OPTIONAL :: ellips_type !< number in the interval [1,nellips] indicating a predefined ellipsoid, alternative to the previous arguments
235CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
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
269!> Destroy a \a griddim_def object.
270SUBROUTINE griddim_delete(this)
271TYPE(griddim_def),INTENT(inout) :: this !< object to be destroyed
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
282!> Create an independent copy of a \a griddim_def object.
283SUBROUTINE griddim_copy(this, that, categoryappend)
284TYPE(griddim_def),INTENT(in) :: this !< object to be copied
285TYPE(griddim_def),INTENT(out) :: that !< copied object
286CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
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
308!> Computes and returns coordinates in the projected system given the
309!! geographical coordinates.
310ELEMENTAL SUBROUTINE griddim_coord_proj(this, lon, lat, x, y)
311TYPE(griddim_def),INTENT(in) :: this !< definition of projection
312!> geographical coordinates
313DOUBLE PRECISION,INTENT(in) :: lon, lat
314!> projected coordinates
315DOUBLE PRECISION,INTENT(out) :: x, y
316
317CALL proj(this%grid%proj, lon, lat, x, y)
318
319END SUBROUTINE griddim_coord_proj
320
321
322!> Computes and returns geographical coordinates given the coordinates
323!! in the projected system.
324ELEMENTAL SUBROUTINE griddim_coord_unproj(this, x, y, lon, lat)
325TYPE(griddim_def),INTENT(in) :: this !< definition of projection
326!> projected coordinates
327DOUBLE PRECISION,INTENT(in) :: x, y
328!> geographical coordinates
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
351!> Computes the geographical coordinates of all the grid points in the
352!! \a griddim_def object and stores them in the object itself. The \a
353!! this::dim::lon and \a this::dim:lat arrays are allocated, and upon
354!! return they will contain the coordinates' values. These arrays can
355!! be directly accessed by the user, and their association status
356!! should be checked with the \a ASSOCIATED() intrinsic before using
357!! them.
358SUBROUTINE griddim_unproj(this)
359TYPE(griddim_def),INTENT(inout) :: this !< definition of projection and grid
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
379!> Query the object content.
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 !< object to be queried
388INTEGER,INTENT(out),OPTIONAL :: nx !< number of points along the x axis
389INTEGER,INTENT(out),OPTIONAL :: ny !< number of points along the y axis
390!> longitudini e latitudini minime e massime
391DOUBLE PRECISION,INTENT(out),OPTIONAL :: xmin, xmax, ymin, ymax !< grid extremes in projection units (degrees or meters depending on the projection type)
392!> grid steps in x and y directions
393DOUBLE PRECISION,INTENT(out),OPTIONAL :: dx, dy
394!> Resolved u- and v- components of vector quantities relative to 0=the easterly and northerly directions
395!! 1=the defined grid in the direction of increasing x and y (or i and j) coordinates respectively (0=north, 128=south)
396INTEGER,INTENT(out),OPTIONAL :: component_flag
397TYPE(geo_proj),INTENT(out),OPTIONAL :: proj !< the complete projection object associated
398CHARACTER(len=*),INTENT(out),OPTIONAL :: proj_type !< type of projection
399DOUBLE PRECISION,INTENT(out),OPTIONAL :: lov !< line of view, also known as reference longitude or orientation of the grid (polar projections)
400INTEGER,INTENT(out),OPTIONAL :: zone !< Earth zone (mainly for UTM), sets lov to the correct zone central meridian
401DOUBLE PRECISION,INTENT(out),OPTIONAL :: xoff !< offset on x axis (false easting)
402DOUBLE PRECISION,INTENT(out),OPTIONAL :: yoff !< offset on y axis (false northing)
403DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_south_pole !< longitude of the southern pole of projection
404DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_south_pole !< latitude of the southern pole of projection
405DOUBLE PRECISION,INTENT(out),OPTIONAL :: angle_rotation !< angle of rotation of projection
406DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_stretch_pole !< longitude of the pole of stretching
407DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_stretch_pole !< latitude of the pole of stretching
408DOUBLE PRECISION,INTENT(out),OPTIONAL :: stretch_factor !< stretching factor
409DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin1 !< first standard latitude from main pole (Lambert)
410DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin2 !< second standard latitude from main pole (Lambert)
411DOUBLE PRECISION,INTENT(out),OPTIONAL :: lad !< latitude at which dx and dy (in m) are specified (Lambert, grib2 only)
412INTEGER,INTENT(out),OPTIONAL :: projection_center_flag !< flag indicating which pole is represented
413DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_smaj_axis !< Earth semi-major axis
414DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_flatt !< Earth flattening
415INTEGER,INTENT(out),OPTIONAL :: ellips_type !< number in the interval [1,nellips] indicating a predefined ellipsoid, alternative to the previous arguments
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
439!> Set the object content.
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 !< object to be queried
448INTEGER,INTENT(in),OPTIONAL :: nx !< number of points along the x axis
449INTEGER,INTENT(in),OPTIONAL :: ny !< number of points along the y axis
450!> longitudini e latitudini minime e massime
451DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin, xmax, ymin, ymax !< grid extremes in projection units (degrees or meters depending on the projection type)
452!> grid steps in x and y directions
453DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx, dy
454!> Resolved u- and v- components of vector quantities relative to 0=the easterly and northerly directions
455!! 1=the defined grid in the direction of increasing x and y (or i and j) coordinates respectively (0=north, 128=south)
456INTEGER,INTENT(in),OPTIONAL :: component_flag
457CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type !< type of projection
458DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov !< line of view, also known as reference longitude or orientation of the grid (polar projections)
459INTEGER,INTENT(in),OPTIONAL :: zone !< Earth zone (mainly for UTM), sets lov to the correct zone central meridian
460DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff !< offset on x axis (false easting)
461DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff !< offset on y axis (false northing)
462DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole !< longitude of the southern pole of projection
463DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole !< latitude of the southern pole of projection
464DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation !< angle of rotation of projection
465DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole !< longitude of the pole of stretching
466DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole !< latitude of the pole of stretching
467DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor !< stretching factor
468DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1 !< first standard latitude from main pole (Lambert)
469DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2 !< second standard latitude from main pole (Lambert)
470DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad !< latitude at which dx and dy (in m) are specified (Lambert, grib2 only)
471INTEGER,INTENT(in),OPTIONAL :: projection_center_flag !< flag indicating which pole is represented
472DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis !< Earth semi-major axis
473DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt !< Earth flattening
474INTEGER,INTENT(in),OPTIONAL :: ellips_type !< number in the interval [1,nellips] indicating a predefined ellipsoid, alternative to the previous arguments
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
495!> This method reads from a Fortran file unit the contents of the
496!! object \a this. The record to be read must have been written with
497!! the ::write_unit method. The method works both on formatted and
498!! unformatted files.
499SUBROUTINE griddim_read_unit(this, unit)
500TYPE(griddim_def),INTENT(out) :: this !< object to be read
501INTEGER, INTENT(in) :: unit !< unit from which to read, it must be an opened Fortran file 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
511!> This method writes on a Fortran file unit the contents of the
512!! object \a this. The record can successively be read by the
513!! ::read_unit method. The method works both on formatted and
514!! unformatted files.
515SUBROUTINE griddim_write_unit(this, unit)
516TYPE(griddim_def),INTENT(in) :: this !< object to be written
517INTEGER, INTENT(in) :: unit !< unit where to write, it must be an opened Fortran file 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
527!> Euristically determine the approximate central longitude of the
528!! grid in degrees.
529!! The method depends on the projection used.
530FUNCTION griddim_central_lon(this) RESULT(lon)
531TYPE(griddim_def),INTENT(inout) :: this !< grid descriptor
532
533DOUBLE PRECISION :: lon
534
535CALL griddim_pistola_central_lon(this, lon)
536
537END FUNCTION griddim_central_lon
538
539
540!> Euristically reset the approximate central longitude of the
541!! grid to a value compatible to the provided longitude \a lonref.
542!! The method depends on the projection used.
543SUBROUTINE griddim_set_central_lon(this, lonref)
544TYPE(griddim_def),INTENT(inout) :: this !< grid descriptor
545DOUBLE PRECISION,INTENT(in) :: lonref !< reference longitude
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
617!> Generates coordinates of every point of a generic grid from the
618!! grid description. The number of grid points along both direction is
619!! guessed from the shape of x and y arrays, which must be conformal.
620SUBROUTINE griddim_gen_coord(this, x, y)
621TYPE(griddim_def),INTENT(in) :: this !< generic grid descriptor
622DOUBLE PRECISION,INTENT(out) :: x(:,:) !< x coordinate of every point, linearly computed between grid extremes
623DOUBLE PRECISION,INTENT(out) :: y(:,:) !< y coordinate of every point, linearly computed between grid extremes, it should have the same shape as x(:,:)
624
625
626CALL grid_rect_coordinates(this%grid%grid, x, y)
627
628END SUBROUTINE griddim_gen_coord
629
631!> Compute and return grid steps.
632SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
633TYPE(griddim_def), INTENT(in) :: this !< generic grid descriptor
634INTEGER,INTENT(in) :: nx !< number of points along x direction
635INTEGER,INTENT(in) :: ny !< number of points along y direction
636DOUBLE PRECISION,INTENT(out) :: dx !< grid step along x direction
637DOUBLE PRECISION,INTENT(out) :: dy !< grid step along y direction
638
639CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
640
641END SUBROUTINE griddim_steps
642
643
644!> Compute and set grid steps.
645SUBROUTINE griddim_setsteps(this)
646TYPE(griddim_def), INTENT(inout) :: this !< generic grid descriptor
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
693!> Import a griddim object from a grid_id object associated to a
694!! supported gridded dataset driver (typically a grib message from
695!! grib_api or a raster band from gdal). The griddim object is
696!! populated with all the grid information (size, projection, etc.)
697!! carried by the grid_id object provided.
698SUBROUTINE griddim_import_grid_id(this, ingrid_id)
699#ifdef HAVE_LIBGDAL
700USE gdal
701#endif
702TYPE(griddim_def),INTENT(inout) :: this !< griddim object
703TYPE(grid_id),INTENT(in) :: ingrid_id !< grid_id object with information about the grid
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
726!> Export a griddim object to a grid_id object associated to a
727!! supported gridded dataset driver (typically a grib message from
728!! grib_api). All the grid information (size, projection, etc.)
729!! contained in the griddim object is exported to the grid_id object.
730SUBROUTINE griddim_export_grid_id(this, outgrid_id)
731#ifdef HAVE_LIBGDAL
732USE gdal
733#endif
734TYPE(griddim_def),INTENT(in) :: this !< griddim object
735TYPE(grid_id),INTENT(inout) :: outgrid_id !< grid_id object which will contain information about the grid
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
1687!> Display on the screen a brief content of the object.
1688SUBROUTINE griddim_display(this)
1689TYPE(griddim_def),INTENT(in) :: this !< griddim object to display
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
1715!> Compute rotation matrix for wind unrotation. It allocates and
1716!! computes a matrix suitable for recomputing wind components in the
1717!! geographical system from the components in the projected
1718!! system. The rotation matrix \a rot_mat has to be passed as a
1719!! pointer and successively deallocated by the caller; it is a
1720!! 3-dimensional array where the first two dimensions are lon and lat
1721!! and the third, with extension 4, contains the packed rotation
1722!! matrix for the given grid point. It should work for every
1723!! projection. In order for the method to work, the \a griddim_unproj
1724!! method must have already been called for the \a griddim_def object.
1725!! \todo Check the algorithm and add some orthogonality tests.
1726SUBROUTINE griddim_wind_unrot(this, rot_mat)
1728TYPE(griddim_def), INTENT(in) :: this !< object describing the grid
1729DOUBLE PRECISION, POINTER :: rot_mat(:,:,:) !< rotation matrix for every grid point, to be deallocated by the caller, if \c .NOT. \c ASSOCIATED() an error occurred
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
1829!> Reset a longitude value in the interval [0-360[.
1830!! The value is reset in place. This is usually useful in connection
1831!! with grib2 coding/decoding.
1832SUBROUTINE long_reset_0_360(lon)
1833DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
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
1846!> Reset a longitude value in the interval [-180-360[.
1847!! The value is reset in place. This is usually useful in connection
1848!! with grib1 coding/decoding.
1849SUBROUTINE long_reset_m180_360(lon)
1850DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
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
1863!> Reset a longitude value in the interval [-90-270[.
1864!! The value is reset in place. This is usually useful in connection
1865!! with grib2 coding/decoding.
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
1880!> Reset a longitude value in the interval [-180-180[.
1881!! The value is reset in place. This is usually useful in connection
1882!! with grib2 coding/decoding.
1883SUBROUTINE long_reset_m180_180(lon)
1884DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
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 !< the longitude to reset
1899DOUBLE PRECISION,INTENT(in) :: lonref !< the longitude to compare
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.