libsim Versione 7.2.6

◆ long_reset_m180_180()

subroutine long_reset_m180_180 ( double precision, intent(inout) lon)
private

Reset a longitude value in the interval [-90-270[.

The value is reset in place. This is usually useful in connection with grib2 coding/decoding. Reset a longitude value in the interval [-180-180[. The value is reset in place. This is usually useful in connection with grib2 coding/decoding.

Parametri
[in,out]lonthe longitude to reset

Definizione alla linea 3149 del file grid_class.F90.

3150! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
3151! authors:
3152! Davide Cesari <dcesari@arpa.emr.it>
3153! Paolo Patruno <ppatruno@arpa.emr.it>
3154
3155! This program is free software; you can redistribute it and/or
3156! modify it under the terms of the GNU General Public License as
3157! published by the Free Software Foundation; either version 2 of
3158! the License, or (at your option) any later version.
3159
3160! This program is distributed in the hope that it will be useful,
3161! but WITHOUT ANY WARRANTY; without even the implied warranty of
3162! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3163! GNU General Public License for more details.
3164
3165! You should have received a copy of the GNU General Public License
3166! along with this program. If not, see <http://www.gnu.org/licenses/>.
3167#include "config.h"
3168!> Module for describing geographically referenced regular grids.
3169!! This module defines classes and methods describing rectangular
3170!! georeferenced grids in different geographical projections. The grid
3171!! and projection definition can be specified explicitely by the
3172!! caller or they can be entirely imported from a grib file (through
3173!! grib_api) or from a file format supported by gdal (through gdal
3174!! fortran interface).
3175!!
3176!! The projection is internally stored following the WMO grib
3177!! conventions (gridType in grib_api). The projections currently
3178!! supported or for which support is planned are:
3179!! - For grib edition 1 and 2:
3180!! - regular_ll (works)
3181!! - mercator (to be done)
3182!! - lambert (works, to be completed)
3183!! - polar_stereographic (to be tested)
3184!! - albers (to be done)
3185!! - rotated_ll (works)
3186!! - stretched_ll (to be completed and tested)
3187!! - stretched_rotated_ll (to be completed and tested)
3188!! - For grib edition 2 only:
3189!! - UTM (ARPA-SIM extension)
3190!! - equatorial_azimuthal_equidistant (to be done)
3191!! - For gdal-supported formats:
3192!! - regular_ll
3193!!
3194!!
3195!! See the example program \include example_vg6d_1.f90
3196!!
3197!!\ingroup volgrid6d
3198MODULE grid_class
3199use geo_proj_class
3201use grid_rect_class
3202use grid_id_class
3203use err_handling
3206use log4fortran
3207implicit none
3208
3209
3210character (len=255),parameter:: subcategory="grid_class"
3211
3212
3213!> This object, mainly for internal use, describes a grid on
3214!! a geographical projection, except the grid dimensions.
3215!! The object is opaque, thus all its members have to be set and
3216!! accessed through the constructor and the ::get_val and ::set_val
3217!! methods.
3218type grid_def
3219 private
3220 type(geo_proj) :: proj
3221 type(grid_rect) :: grid
3222 integer :: category = 0
3223end type grid_def
3224
3225
3226!> This object completely describes a grid on a geographic projection.
3227!! It is the main public object of this module. The grid definition \a
3228!! grid is separated from the definition of the grid dimensions \a dim
3229!! in order to make members of \a grid \a PRIVATE while maintaining
3230!! free access to the members of \a dim.
3231type griddim_def
3232 type(grid_def) :: grid !< grid and projection definition
3233 type(grid_dim) :: dim !< grid dimensions definition
3234 integer :: category = 0 !< category for log4fortran
3235end type griddim_def
3236
3237
3238!> Logical equality operators for objects of the classes \a grid_def,
3239!! and \a griddim_def. They are all defined as \c
3240!! ELEMENTAL thus work also on arrays of any shape.
3241INTERFACE OPERATOR (==)
3242 MODULE PROCEDURE grid_eq, griddim_eq
3243END INTERFACE
3244
3245!> Logical inequality operators for objects of the classes \a grid_def,
3246!! and \a griddim_def. They are all defined as \c
3247!! ELEMENTAL thus work also on arrays of any shape.
3248INTERFACE OPERATOR (/=)
3249 MODULE PROCEDURE grid_ne, griddim_ne
3250END INTERFACE
3251
3252!> Constructors of the corresponding objects.
3253INTERFACE init
3254 MODULE PROCEDURE griddim_init
3255END INTERFACE
3256
3257!> Destructors of the corresponding objects.
3258INTERFACE delete
3259 MODULE PROCEDURE griddim_delete
3260END INTERFACE
3261
3262!> Copy an object, creating a fully new instance.
3263INTERFACE copy
3264 MODULE PROCEDURE griddim_copy
3265END INTERFACE
3266
3267!> Compute forward coordinate transformation from geographical system to
3268!! projected system.
3269INTERFACE proj
3270 MODULE PROCEDURE griddim_coord_proj!, griddim_proj
3271END INTERFACE
3272
3273!> Compute backward coordinate transformation from projected system
3274!! to geographical system.
3275INTERFACE unproj
3276 MODULE PROCEDURE griddim_coord_unproj, griddim_unproj
3277END INTERFACE
3278
3279!> Method for returning the contents of the object.
3280INTERFACE get_val
3281 MODULE PROCEDURE griddim_get_val
3282END INTERFACE
3283
3284!> Method for setting the contents of the object.
3285INTERFACE set_val
3286 MODULE PROCEDURE griddim_set_val
3287END INTERFACE
3288
3289!> Write the object on a formatted or unformatted file.
3290INTERFACE write_unit
3291 MODULE PROCEDURE griddim_write_unit
3292END INTERFACE
3293
3294!> Read the object from a formatted or unformatted file.
3295INTERFACE read_unit
3296 MODULE PROCEDURE griddim_read_unit
3297END INTERFACE
3298
3299!> Import griddim object from grid_id.
3300INTERFACE import
3301 MODULE PROCEDURE griddim_import_grid_id
3302END INTERFACE
3303
3304!> Export griddim object to grid_id.
3305INTERFACE export
3306 MODULE PROCEDURE griddim_export_grid_id
3307END INTERFACE
3308
3309!> Print a brief description on stdout.
3310INTERFACE display
3311 MODULE PROCEDURE griddim_display
3312END INTERFACE
3313
3314#define VOL7D_POLY_TYPE TYPE(grid_def)
3315#define VOL7D_POLY_TYPES _grid
3316#include "array_utilities_pre.F90"
3317#undef VOL7D_POLY_TYPE
3318#undef VOL7D_POLY_TYPES
3319
3320#define VOL7D_POLY_TYPE TYPE(griddim_def)
3321#define VOL7D_POLY_TYPES _griddim
3322#include "array_utilities_pre.F90"
3323#undef VOL7D_POLY_TYPE
3324#undef VOL7D_POLY_TYPES
3325
3326INTERFACE wind_unrot
3327 MODULE PROCEDURE griddim_wind_unrot
3328END INTERFACE
3329
3330
3331PRIVATE
3332
3333PUBLIC proj, unproj, griddim_unproj, griddim_gen_coord, &
3334 griddim_zoom_coord, griddim_zoom_projcoord, &
3335 griddim_setsteps, griddim_def, grid_def, grid_dim
3336PUBLIC init, delete, copy
3338PUBLIC OPERATOR(==),OPERATOR(/=)
3339PUBLIC count_distinct, pack_distinct, count_and_pack_distinct, &
3340 map_distinct, map_inv_distinct,index
3341PUBLIC wind_unrot, import, export
3342PUBLIC griddim_central_lon, griddim_set_central_lon
3343CONTAINS
3344
3345!> Constructor for a \a griddim_def object.
3346SUBROUTINE griddim_init(this, nx, ny, &
3347 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
3348 proj_type, lov, zone, xoff, yoff, &
3349 longitude_south_pole, latitude_south_pole, angle_rotation, &
3350 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
3351 latin1, latin2, lad, projection_center_flag, &
3352 ellips_smaj_axis, ellips_flatt, ellips_type, &
3353 categoryappend)
3354TYPE(griddim_def),INTENT(inout) :: this !< object to be created
3355INTEGER,INTENT(in),OPTIONAL :: nx !< number of points along the x axis
3356INTEGER,INTENT(in),OPTIONAL :: ny !< number of points along the y axis
3357DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin !< lower bound for x coordinate on grid in projection units (degrees or meters depending on the projection type)
3358DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmax !< upper bound for x coordinate
3359DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymin !< lower bound for y coordinate
3360DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymax !< upper bound for y coordinate
3361DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx !< grid step in x direction
3362DOUBLE PRECISION,INTENT(in),OPTIONAL :: dy !< grid step in y direction
3363!> Resolved u- and v- components of vector quantities relative to 0=the easterly and northerly directions
3364!! 1=the defined grid in the direction of increasing x and y (or i and j) coordinates respectively (0=north, 128=south)
3365INTEGER,INTENT(in),OPTIONAL :: component_flag
3366CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type !< type of projection
3367DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov !< line of view, also known as reference longitude or orientation of the grid (polar projections)
3368INTEGER,INTENT(in),OPTIONAL :: zone !< Earth zone (mainly for UTM), sets lov to the correct zone central meridian
3369DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff !< offset on x axis (false easting)
3370DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff !< offset on y axis (false northing)
3371DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole !< longitude of the southern pole of projection
3372DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole !< latitude of the southern pole of projection
3373DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation !< angle of rotation of projection
3374DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole !< longitude of the pole of stretching
3375DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole !< latitude of the pole of stretching
3376DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor !< stretching factor
3377DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1 !< first standard latitude from main pole (Lambert)
3378DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2 !< second standard latitude from main pole (Lambert)
3379DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad !< latitude at which dx and dy (in m) are specified (Lambert, grib2 only)
3380INTEGER,INTENT(in),OPTIONAL :: projection_center_flag !< flag indicating which pole is represented
3381DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis !< Earth semi-major axis
3382DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt !< Earth flattening
3383INTEGER,INTENT(in),OPTIONAL :: ellips_type !< number in the interval [1,nellips] indicating a predefined ellipsoid, alternative to the previous arguments
3384CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
3385
3386CHARACTER(len=512) :: a_name
3387
3388IF (PRESENT(categoryappend)) THEN
3389 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
3390 trim(categoryappend))
3391ELSE
3392 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
3393ENDIF
3394this%category=l4f_category_get(a_name)
3395
3396! geographical projection
3397this%grid%proj = geo_proj_new( &
3398 proj_type=proj_type, lov=lov, zone=zone, xoff=xoff, yoff=yoff, &
3399 longitude_south_pole=longitude_south_pole, &
3400 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
3401 longitude_stretch_pole=longitude_stretch_pole, &
3402 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
3403 latin1=latin1, latin2=latin2, lad=lad, projection_center_flag=projection_center_flag, &
3404 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, ellips_type=ellips_type)
3405! grid extension
3406this%grid%grid = grid_rect_new( &
3407 xmin, xmax, ymin, ymax, dx, dy, component_flag)
3408! grid size
3409this%dim = grid_dim_new(nx, ny)
3410
3411#ifdef DEBUG
3412call l4f_category_log(this%category,l4f_debug,"init gtype: "//this%grid%proj%proj_type )
3413#endif
3414
3415END SUBROUTINE griddim_init
3416
3417
3418!> Destroy a \a griddim_def object.
3419SUBROUTINE griddim_delete(this)
3420TYPE(griddim_def),INTENT(inout) :: this !< object to be destroyed
3421
3422CALL delete(this%dim)
3423CALL delete(this%grid%proj)
3424CALL delete(this%grid%grid)
3425
3426call l4f_category_delete(this%category)
3427
3428END SUBROUTINE griddim_delete
3429
3430
3431!> Create an independent copy of a \a griddim_def object.
3432SUBROUTINE griddim_copy(this, that, categoryappend)
3433TYPE(griddim_def),INTENT(in) :: this !< object to be copied
3434TYPE(griddim_def),INTENT(out) :: that !< copied object
3435CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
3436
3437CHARACTER(len=512) :: a_name
3438
3439CALL init(that)
3440
3441CALL copy(this%grid%proj, that%grid%proj)
3442CALL copy(this%grid%grid, that%grid%grid)
3443CALL copy(this%dim, that%dim)
3444
3445! new category
3446IF (PRESENT(categoryappend)) THEN
3447 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
3448 trim(categoryappend))
3449ELSE
3450 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
3451ENDIF
3452that%category=l4f_category_get(a_name)
3453
3454END SUBROUTINE griddim_copy
3455
3456
3457!> Computes and returns coordinates in the projected system given the
3458!! geographical coordinates.
3459ELEMENTAL SUBROUTINE griddim_coord_proj(this, lon, lat, x, y)
3460TYPE(griddim_def),INTENT(in) :: this !< definition of projection
3461!> geographical coordinates
3462DOUBLE PRECISION,INTENT(in) :: lon, lat
3463!> projected coordinates
3464DOUBLE PRECISION,INTENT(out) :: x, y
3465
3466CALL proj(this%grid%proj, lon, lat, x, y)
3467
3468END SUBROUTINE griddim_coord_proj
3469
3470
3471!> Computes and returns geographical coordinates given the coordinates
3472!! in the projected system.
3473ELEMENTAL SUBROUTINE griddim_coord_unproj(this, x, y, lon, lat)
3474TYPE(griddim_def),INTENT(in) :: this !< definition of projection
3475!> projected coordinates
3476DOUBLE PRECISION,INTENT(in) :: x, y
3477!> geographical coordinates
3478DOUBLE PRECISION,INTENT(out) :: lon, lat
3479
3480CALL unproj(this%grid%proj, x, y, lon, lat)
3481
3482END SUBROUTINE griddim_coord_unproj
3483
3484
3485! Computes and sets the grid parameters required to compute
3486! coordinates of grid points in the projected system.
3487! probably meaningless
3488!SUBROUTINE griddim_proj(this)
3489!TYPE(griddim_def),INTENT(inout) :: this !< definition of projection and grid
3490!
3491!CALL proj(this, this%dim%lon(1,1), this%dim%lat(1,1), &
3492! this%grid%grid%xmin, this%grid%grid%ymin)
3493!
3494!CALL proj(this, this%dim%lon(this%dim%nx,this%dim%ny), &
3495! this%dim%lat(this%dim%nx,this%dim%ny), &
3496! this%grid%grid%xmax, this%grid%grid%ymax)
3497!
3498!END SUBROUTINE griddim_proj
3499
3500!> Computes the geographical coordinates of all the grid points in the
3501!! \a griddim_def object and stores them in the object itself. The \a
3502!! this::dim::lon and \a this::dim:lat arrays are allocated, and upon
3503!! return they will contain the coordinates' values. These arrays can
3504!! be directly accessed by the user, and their association status
3505!! should be checked with the \a ASSOCIATED() intrinsic before using
3506!! them.
3507SUBROUTINE griddim_unproj(this)
3508TYPE(griddim_def),INTENT(inout) :: this !< definition of projection and grid
3509
3510IF (.NOT.c_e(this%dim%nx) .OR. .NOT.c_e(this%dim%ny)) RETURN
3511CALL alloc(this%dim)
3512CALL griddim_unproj_internal(this)
3513
3514END SUBROUTINE griddim_unproj
3515
3516! internal subroutine needed for allocating automatic arrays
3517SUBROUTINE griddim_unproj_internal(this)
3518TYPE(griddim_def),INTENT(inout) ::this ! definition of projection and grid
3519
3520DOUBLE PRECISION :: x(this%dim%nx,this%dim%ny), y(this%dim%nx,this%dim%ny)
3521
3522CALL grid_rect_coordinates(this%grid%grid, x, y)
3523CALL unproj(this, x, y, this%dim%lon, this%dim%lat)
3524
3525END SUBROUTINE griddim_unproj_internal
3526
3527
3528!> Query the object content.
3529SUBROUTINE griddim_get_val(this, nx, ny, &
3530 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
3531 proj, proj_type, lov, zone, xoff, yoff, &
3532 longitude_south_pole, latitude_south_pole, angle_rotation, &
3533 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
3534 latin1, latin2, lad, projection_center_flag, &
3535 ellips_smaj_axis, ellips_flatt, ellips_type)
3536TYPE(griddim_def),INTENT(in) :: this !< object to be queried
3537INTEGER,INTENT(out),OPTIONAL :: nx !< number of points along the x axis
3538INTEGER,INTENT(out),OPTIONAL :: ny !< number of points along the y axis
3539!> longitudini e latitudini minime e massime
3540DOUBLE PRECISION,INTENT(out),OPTIONAL :: xmin, xmax, ymin, ymax !< grid extremes in projection units (degrees or meters depending on the projection type)
3541!> grid steps in x and y directions
3542DOUBLE PRECISION,INTENT(out),OPTIONAL :: dx, dy
3543!> Resolved u- and v- components of vector quantities relative to 0=the easterly and northerly directions
3544!! 1=the defined grid in the direction of increasing x and y (or i and j) coordinates respectively (0=north, 128=south)
3545INTEGER,INTENT(out),OPTIONAL :: component_flag
3546TYPE(geo_proj),INTENT(out),OPTIONAL :: proj !< the complete projection object associated
3547CHARACTER(len=*),INTENT(out),OPTIONAL :: proj_type !< type of projection
3548DOUBLE PRECISION,INTENT(out),OPTIONAL :: lov !< line of view, also known as reference longitude or orientation of the grid (polar projections)
3549INTEGER,INTENT(out),OPTIONAL :: zone !< Earth zone (mainly for UTM), sets lov to the correct zone central meridian
3550DOUBLE PRECISION,INTENT(out),OPTIONAL :: xoff !< offset on x axis (false easting)
3551DOUBLE PRECISION,INTENT(out),OPTIONAL :: yoff !< offset on y axis (false northing)
3552DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_south_pole !< longitude of the southern pole of projection
3553DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_south_pole !< latitude of the southern pole of projection
3554DOUBLE PRECISION,INTENT(out),OPTIONAL :: angle_rotation !< angle of rotation of projection
3555DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_stretch_pole !< longitude of the pole of stretching
3556DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_stretch_pole !< latitude of the pole of stretching
3557DOUBLE PRECISION,INTENT(out),OPTIONAL :: stretch_factor !< stretching factor
3558DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin1 !< first standard latitude from main pole (Lambert)
3559DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin2 !< second standard latitude from main pole (Lambert)
3560DOUBLE PRECISION,INTENT(out),OPTIONAL :: lad !< latitude at which dx and dy (in m) are specified (Lambert, grib2 only)
3561INTEGER,INTENT(out),OPTIONAL :: projection_center_flag !< flag indicating which pole is represented
3562DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_smaj_axis !< Earth semi-major axis
3563DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_flatt !< Earth flattening
3564INTEGER,INTENT(out),OPTIONAL :: ellips_type !< number in the interval [1,nellips] indicating a predefined ellipsoid, alternative to the previous arguments
3565
3566IF (PRESENT(nx)) nx = this%dim%nx
3567IF (PRESENT(ny)) ny = this%dim%ny
3568
3569IF (PRESENT(proj)) proj = this%grid%proj
3570
3571CALL get_val(this%grid%proj, proj_type=proj_type, lov=lov, zone=zone, &
3572 xoff=xoff, yoff=yoff, &
3573 longitude_south_pole=longitude_south_pole, &
3574 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
3575 longitude_stretch_pole=longitude_stretch_pole, &
3576 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
3577 latin1=latin1, latin2=latin2, lad=lad, &
3578 projection_center_flag=projection_center_flag, &
3579 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
3580 ellips_type=ellips_type)
3581
3582CALL get_val(this%grid%grid, &
3583 xmin, xmax, ymin, ymax, dx, dy, component_flag)
3584
3585END SUBROUTINE griddim_get_val
3586
3587
3588!> Set the object content.
3589SUBROUTINE griddim_set_val(this, nx, ny, &
3590 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
3591 proj_type, lov, zone, xoff, yoff, &
3592 longitude_south_pole, latitude_south_pole, angle_rotation, &
3593 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
3594 latin1, latin2, lad, projection_center_flag, &
3595 ellips_smaj_axis, ellips_flatt, ellips_type)
3596TYPE(griddim_def),INTENT(inout) :: this !< object to be queried
3597INTEGER,INTENT(in),OPTIONAL :: nx !< number of points along the x axis
3598INTEGER,INTENT(in),OPTIONAL :: ny !< number of points along the y axis
3599!> longitudini e latitudini minime e massime
3600DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin, xmax, ymin, ymax !< grid extremes in projection units (degrees or meters depending on the projection type)
3601!> grid steps in x and y directions
3602DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx, dy
3603!> Resolved u- and v- components of vector quantities relative to 0=the easterly and northerly directions
3604!! 1=the defined grid in the direction of increasing x and y (or i and j) coordinates respectively (0=north, 128=south)
3605INTEGER,INTENT(in),OPTIONAL :: component_flag
3606CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type !< type of projection
3607DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov !< line of view, also known as reference longitude or orientation of the grid (polar projections)
3608INTEGER,INTENT(in),OPTIONAL :: zone !< Earth zone (mainly for UTM), sets lov to the correct zone central meridian
3609DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff !< offset on x axis (false easting)
3610DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff !< offset on y axis (false northing)
3611DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole !< longitude of the southern pole of projection
3612DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole !< latitude of the southern pole of projection
3613DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation !< angle of rotation of projection
3614DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole !< longitude of the pole of stretching
3615DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole !< latitude of the pole of stretching
3616DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor !< stretching factor
3617DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1 !< first standard latitude from main pole (Lambert)
3618DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2 !< second standard latitude from main pole (Lambert)
3619DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad !< latitude at which dx and dy (in m) are specified (Lambert, grib2 only)
3620INTEGER,INTENT(in),OPTIONAL :: projection_center_flag !< flag indicating which pole is represented
3621DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis !< Earth semi-major axis
3622DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt !< Earth flattening
3623INTEGER,INTENT(in),OPTIONAL :: ellips_type !< number in the interval [1,nellips] indicating a predefined ellipsoid, alternative to the previous arguments
3624
3625IF (PRESENT(nx)) this%dim%nx = nx
3626IF (PRESENT(ny)) this%dim%ny = ny
3627
3628CALL set_val(this%grid%proj, proj_type=proj_type, lov=lov, zone=zone, &
3629 xoff=xoff, yoff=yoff, longitude_south_pole=longitude_south_pole, &
3630 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
3631 longitude_stretch_pole=longitude_stretch_pole, &
3632 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
3633 latin1=latin1, latin2=latin2, lad=lad, &
3634 projection_center_flag=projection_center_flag, &
3635 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
3636 ellips_type=ellips_type)
3637
3638CALL set_val(this%grid%grid, &
3639 xmin, xmax, ymin, ymax, dx, dy, component_flag)
3640
3641END SUBROUTINE griddim_set_val
3642
3643
3644!> This method reads from a Fortran file unit the contents of the
3645!! object \a this. The record to be read must have been written with
3646!! the ::write_unit method. The method works both on formatted and
3647!! unformatted files.
3648SUBROUTINE griddim_read_unit(this, unit)
3649TYPE(griddim_def),INTENT(out) :: this !< object to be read
3650INTEGER, INTENT(in) :: unit !< unit from which to read, it must be an opened Fortran file unit
3651
3652
3653CALL read_unit(this%dim, unit)
3654CALL read_unit(this%grid%proj, unit)
3655CALL read_unit(this%grid%grid, unit)
3656
3657END SUBROUTINE griddim_read_unit
3658
3659
3660!> This method writes on a Fortran file unit the contents of the
3661!! object \a this. The record can successively be read by the
3662!! ::read_unit method. The method works both on formatted and
3663!! unformatted files.
3664SUBROUTINE griddim_write_unit(this, unit)
3665TYPE(griddim_def),INTENT(in) :: this !< object to be written
3666INTEGER, INTENT(in) :: unit !< unit where to write, it must be an opened Fortran file unit
3667
3668
3669CALL write_unit(this%dim, unit)
3670CALL write_unit(this%grid%proj, unit)
3671CALL write_unit(this%grid%grid, unit)
3672
3673END SUBROUTINE griddim_write_unit
3674
3675
3676!> Euristically determine the approximate central longitude of the
3677!! grid in degrees.
3678!! The method depends on the projection used.
3679FUNCTION griddim_central_lon(this) RESULT(lon)
3680TYPE(griddim_def),INTENT(inout) :: this !< grid descriptor
3681
3682DOUBLE PRECISION :: lon
3683
3684CALL griddim_pistola_central_lon(this, lon)
3685
3686END FUNCTION griddim_central_lon
3687
3688
3689!> Euristically reset the approximate central longitude of the
3690!! grid to a value compatible to the provided longitude \a lonref.
3691!! The method depends on the projection used.
3692SUBROUTINE griddim_set_central_lon(this, lonref)
3693TYPE(griddim_def),INTENT(inout) :: this !< grid descriptor
3694DOUBLE PRECISION,INTENT(in) :: lonref !< reference longitude
3695
3696DOUBLE PRECISION :: lon
3697
3698CALL griddim_pistola_central_lon(this, lon, lonref)
3699
3700END SUBROUTINE griddim_set_central_lon
3701
3702
3703! internal subroutine for performing tasks common to the prevous two
3704SUBROUTINE griddim_pistola_central_lon(this, lon, lonref)
3705TYPE(griddim_def),INTENT(inout) :: this ! grid descriptor
3706DOUBLE PRECISION,INTENT(inout) :: lon ! central longitude
3707DOUBLE PRECISION,INTENT(in),OPTIONAL :: lonref ! reference longitude
3708
3709INTEGER :: unit
3710DOUBLE PRECISION :: lonsp, latsp, londelta, lov, lonrot
3711CHARACTER(len=80) :: ptype
3712
3713lon = dmiss
3714CALL get_val(this%grid%proj, unit=unit)
3715IF (unit == geo_proj_unit_meter) THEN ! it is a plane projection
3716 CALL get_val(this%grid%proj, lov=lon)
3717 IF (PRESENT(lonref)) THEN
3718 CALL long_reset_to_cart_closest(lov, lonref)
3719 CALL set_val(this%grid%proj, lov=lon)
3720 ENDIF
3721
3722ELSE IF (unit == geo_proj_unit_degree) THEN ! it is a spheric projection
3723 CALL get_val(this%grid%proj, proj_type=ptype, &
3724 longitude_south_pole=lonsp, latitude_south_pole=latsp)
3725 SELECT CASE(ptype)
3726 CASE('rotated_ll','stretched_rotated_ll') ! use origin of rotated system
3727 IF (latsp < 0.0d0) THEN
3728 lon = lonsp
3729 IF (PRESENT(lonref)) THEN
3730 CALL long_reset_to_cart_closest(lon, lonref)
3731 CALL set_val(this%grid%proj, longitude_south_pole=lonref)
3732! now reset rotated coordinates around zero
3733 IF (c_e(this%grid%grid%xmin) .AND. c_e(this%grid%grid%xmax)) THEN
3734 lonrot = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
3735 ENDIF
3736 londelta = lonrot
3737 CALL long_reset_to_cart_closest(londelta, 0.0d0)
3738 londelta = londelta - lonrot
3739 this%grid%grid%xmin = this%grid%grid%xmin + londelta
3740 this%grid%grid%xmax = this%grid%grid%xmax + londelta
3741 ENDIF
3742 ELSE ! this part to be checked
3743 lon = modulo(lonsp + 180.0d0, 360.0d0)
3744! IF (PRESENT(lonref)) THEN
3745! CALL long_reset_to_cart_closest(lov, lonref)
3746! CALL set_val(this%grid%proj, longitude_south_pole=lonref)
3747! ENDIF
3748 ENDIF
3749 CASE default ! use real grid limits
3750 IF (c_e(this%grid%grid%xmin) .AND. c_e(this%grid%grid%xmax)) THEN
3751 lon = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
3752 ENDIF
3753 IF (PRESENT(lonref)) THEN
3754 londelta = lon
3755 CALL long_reset_to_cart_closest(londelta, lonref)
3756 londelta = londelta - lon
3757 this%grid%grid%xmin = this%grid%grid%xmin + londelta
3758 this%grid%grid%xmax = this%grid%grid%xmax + londelta
3759 ENDIF
3760 END SELECT
3761ENDIF
3762
3763END SUBROUTINE griddim_pistola_central_lon
3764
3765
3766!> Generates coordinates of every point of a generic grid from the
3767!! grid description. The number of grid points along both direction is
3768!! guessed from the shape of x and y arrays, which must be conformal.
3769SUBROUTINE griddim_gen_coord(this, x, y)
3770TYPE(griddim_def),INTENT(in) :: this !< generic grid descriptor
3771DOUBLE PRECISION,INTENT(out) :: x(:,:) !< x coordinate of every point, linearly computed between grid extremes
3772DOUBLE PRECISION,INTENT(out) :: y(:,:) !< y coordinate of every point, linearly computed between grid extremes, it should have the same shape as x(:,:)
3773
3774
3775CALL grid_rect_coordinates(this%grid%grid, x, y)
3776
3777END SUBROUTINE griddim_gen_coord
3778
3779
3780!> Compute and return grid steps.
3781SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
3782TYPE(griddim_def), INTENT(in) :: this !< generic grid descriptor
3783INTEGER,INTENT(in) :: nx !< number of points along x direction
3784INTEGER,INTENT(in) :: ny !< number of points along y direction
3785DOUBLE PRECISION,INTENT(out) :: dx !< grid step along x direction
3786DOUBLE PRECISION,INTENT(out) :: dy !< grid step along y direction
3787
3788CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
3789
3790END SUBROUTINE griddim_steps
3791
3792
3793!> Compute and set grid steps.
3794SUBROUTINE griddim_setsteps(this)
3795TYPE(griddim_def), INTENT(inout) :: this !< generic grid descriptor
3796
3797CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3798
3799END SUBROUTINE griddim_setsteps
3800
3801
3802! TODO
3803! bisogna sviluppare gli altri operatori
3804ELEMENTAL FUNCTION grid_eq(this, that) RESULT(res)
3805TYPE(grid_def),INTENT(IN) :: this, that
3806LOGICAL :: res
3807
3808res = this%proj == that%proj .AND. &
3809 this%grid == that%grid
3810
3811END FUNCTION grid_eq
3812
3813
3814ELEMENTAL FUNCTION griddim_eq(this, that) RESULT(res)
3815TYPE(griddim_def),INTENT(IN) :: this, that
3816LOGICAL :: res
3817
3818res = this%grid == that%grid .AND. &
3819 this%dim == that%dim
3820
3821END FUNCTION griddim_eq
3822
3823
3824ELEMENTAL FUNCTION grid_ne(this, that) RESULT(res)
3825TYPE(grid_def),INTENT(IN) :: this, that
3826LOGICAL :: res
3827
3828res = .NOT.(this == that)
3829
3830END FUNCTION grid_ne
3831
3832
3833ELEMENTAL FUNCTION griddim_ne(this, that) RESULT(res)
3834TYPE(griddim_def),INTENT(IN) :: this, that
3835LOGICAL :: res
3836
3837res = .NOT.(this == that)
3838
3839END FUNCTION griddim_ne
3840
3841
3842!> Import a griddim object from a grid_id object associated to a
3843!! supported gridded dataset driver (typically a grib message from
3844!! grib_api or a raster band from gdal). The griddim object is
3845!! populated with all the grid information (size, projection, etc.)
3846!! carried by the grid_id object provided.
3847SUBROUTINE griddim_import_grid_id(this, ingrid_id)
3848#ifdef HAVE_LIBGDAL
3849USE gdal
3850#endif
3851TYPE(griddim_def),INTENT(inout) :: this !< griddim object
3852TYPE(grid_id),INTENT(in) :: ingrid_id !< grid_id object with information about the grid
3853
3854#ifdef HAVE_LIBGRIBAPI
3855INTEGER :: gaid
3856#endif
3857#ifdef HAVE_LIBGDAL
3858TYPE(gdalrasterbandh) :: gdalid
3859#endif
3860CALL init(this)
3861
3862#ifdef HAVE_LIBGRIBAPI
3863gaid = grid_id_get_gaid(ingrid_id)
3864IF (c_e(gaid)) CALL griddim_import_gribapi(this, gaid)
3865#endif
3866#ifdef HAVE_LIBGDAL
3867gdalid = grid_id_get_gdalid(ingrid_id)
3868IF (gdalassociated(gdalid)) CALL griddim_import_gdal(this, gdalid, &
3869 grid_id_get_gdal_options(ingrid_id))
3870#endif
3871
3872END SUBROUTINE griddim_import_grid_id
3873
3874
3875!> Export a griddim object to a grid_id object associated to a
3876!! supported gridded dataset driver (typically a grib message from
3877!! grib_api). All the grid information (size, projection, etc.)
3878!! contained in the griddim object is exported to the grid_id object.
3879SUBROUTINE griddim_export_grid_id(this, outgrid_id)
3880#ifdef HAVE_LIBGDAL
3881USE gdal
3882#endif
3883TYPE(griddim_def),INTENT(in) :: this !< griddim object
3884TYPE(grid_id),INTENT(inout) :: outgrid_id !< grid_id object which will contain information about the grid
3885
3886#ifdef HAVE_LIBGRIBAPI
3887INTEGER :: gaid
3888#endif
3889#ifdef HAVE_LIBGDAL
3890TYPE(gdalrasterbandh) :: gdalid
3891#endif
3892
3893#ifdef HAVE_LIBGRIBAPI
3894gaid = grid_id_get_gaid(outgrid_id)
3895IF (c_e(gaid)) CALL griddim_export_gribapi(this, gaid)
3896#endif
3897#ifdef HAVE_LIBGDAL
3898gdalid = grid_id_get_gdalid(outgrid_id)
3899!IF (gdalassociated(gdalid)
3900! export for gdal not implemented, log?
3901#endif
3902
3903END SUBROUTINE griddim_export_grid_id
3904
3905
3906#ifdef HAVE_LIBGRIBAPI
3907! grib_api driver
3908SUBROUTINE griddim_import_gribapi(this, gaid)
3909USE grib_api
3910TYPE(griddim_def),INTENT(inout) :: this ! griddim object
3911INTEGER, INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
3912
3913DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, x1, y1, orient
3914INTEGER :: EditionNumber, iScansNegatively, jScansPositively, zone, datum, &
3915 reflon, ierr
3916
3917! Generic keys
3918CALL grib_get(gaid, 'typeOfGrid', this%grid%proj%proj_type)
3919#ifdef DEBUG
3920call l4f_category_log(this%category,l4f_debug, &
3921 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type))
3922#endif
3923CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
3924
3925IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
3926 CALL grib_get(gaid,'numberOfDataPoints', this%dim%nx)
3927 this%dim%ny = 1
3928 this%grid%grid%component_flag = 0
3929 CALL griddim_import_ellipsoid(this, gaid)
3930 RETURN
3931ENDIF
3932
3933! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
3934CALL grib_get(gaid, 'Ni', this%dim%nx)
3935CALL grib_get(gaid, 'Nj', this%dim%ny)
3936CALL griddim_import_ellipsoid(this, gaid) ! placed here, not valid for utm datum /= 1
3937
3938CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
3939CALL grib_get(gaid,'jScansPositively',jscanspositively)
3940CALL grib_get(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
3941
3942! Keys for rotated grids (checked through missing values)
3943CALL grib_get_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
3944 this%grid%proj%rotated%longitude_south_pole)
3945CALL grib_get_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
3946 this%grid%proj%rotated%latitude_south_pole)
3947CALL grib_get_dmiss(gaid,'angleOfRotationInDegrees', &
3948 this%grid%proj%rotated%angle_rotation)
3949
3950! Keys for stretched grids (checked through missing values)
3951! units must be verified, still experimental in grib_api
3952! # TODO: Is it a float? Is it signed?
3953IF (editionnumber == 1) THEN
3954 CALL grib_get_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
3955 this%grid%proj%stretched%longitude_stretch_pole)
3956 CALL grib_get_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
3957 this%grid%proj%stretched%latitude_stretch_pole)
3958 CALL grib_get_dmiss(gaid,'stretchingFactor', &
3959 this%grid%proj%stretched%stretch_factor)
3960ELSE IF (editionnumber == 2) THEN
3961 CALL grib_get_dmiss(gaid,'longitudeOfThePoleOfStretching', &
3962 this%grid%proj%stretched%longitude_stretch_pole)
3963 CALL grib_get_dmiss(gaid,'latitudeOfThePoleOfStretching', &
3964 this%grid%proj%stretched%latitude_stretch_pole)
3965 CALL grib_get_dmiss(gaid,'stretchingFactor', &
3966 this%grid%proj%stretched%stretch_factor)
3967 IF (c_e(this%grid%proj%stretched%stretch_factor)) &
3968 this%grid%proj%stretched%stretch_factor = this%grid%proj%stretched%stretch_factor*1.0d-6
3969ENDIF
3970
3971! Projection-dependent keys
3972SELECT CASE (this%grid%proj%proj_type)
3973
3974! Keys for spherical coordinate systems
3975CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
3976
3977 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3978 CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3979 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3980 CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3981
3982! longitudes are sometimes wrongly coded even in grib2 and even by the
3983! Metoffice!
3984! longitudeOfFirstGridPointInDegrees = 354.911;
3985! longitudeOfLastGridPointInDegrees = 363.311;
3986 CALL long_reset_0_360(lofirst)
3987 CALL long_reset_0_360(lolast)
3988
3989 IF (iscansnegatively == 0) THEN
3990 this%grid%grid%xmin = lofirst
3991 this%grid%grid%xmax = lolast
3992 ELSE
3993 this%grid%grid%xmax = lofirst
3994 this%grid%grid%xmin = lolast
3995 ENDIF
3996 IF (jscanspositively == 0) THEN
3997 this%grid%grid%ymax = lafirst
3998 this%grid%grid%ymin = lalast
3999 ELSE
4000 this%grid%grid%ymin = lafirst
4001 this%grid%grid%ymax = lalast
4002 ENDIF
4003
4004! reset longitudes in order to have a Cartesian plane
4005 IF (this%grid%grid%xmax < this%grid%grid%xmin) &
4006 this%grid%grid%xmin = this%grid%grid%xmin - 360.d0
4007
4008! compute dx and dy (should we get them from grib?)
4009 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
4010
4011! Keys for polar projections
4012CASE ('polar_stereographic', 'lambert', 'albers')
4013
4014 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
4015 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
4016! latin1/latin2 may be missing (e.g. stereographic)
4017 CALL grib_get_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
4018 CALL grib_get_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
4019#ifdef DEBUG
4020 CALL l4f_category_log(this%category,l4f_debug, &
4021 "griddim_import_gribapi, latin1/2 "// &
4022 trim(to_char(this%grid%proj%polar%latin1))//" "// &
4023 trim(to_char(this%grid%proj%polar%latin2)))
4024#endif
4025! projection center flag, aka hemisphere
4026 CALL grib_get(gaid,'projectionCenterFlag',&
4027 this%grid%proj%polar%projection_center_flag, ierr)
4028 IF (ierr /= grib_success) THEN ! try center/centre
4029 CALL grib_get(gaid,'projectionCentreFlag',&
4030 this%grid%proj%polar%projection_center_flag)
4031 ENDIF
4032
4033 IF (iand(this%grid%proj%polar%projection_center_flag,64) == 1) THEN
4034 CALL l4f_category_log(this%category,l4f_error, &
4035 "griddim_import_gribapi, bi-polar projections not supported")
4036 CALL raise_error()
4037 ENDIF
4038! line of view, aka central meridian
4039 CALL grib_get(gaid,'LoVInDegrees',this%grid%proj%lov)
4040#ifdef DEBUG
4041 CALL l4f_category_log(this%category,l4f_debug, &
4042 "griddim_import_gribapi, central meridian "//trim(to_char(this%grid%proj%lov)))
4043#endif
4044
4045! latitude at which dx and dy are valid
4046 IF (editionnumber == 1) THEN
4047! ECMWF (gribex/grib_api) says: Grid lengths are in metres, at the
4048! 60-degree parallel nearest to the pole on the projection plane.
4049! IF (IAND(this%projection_center_flag, 128) == 0) THEN
4050! this%grid%proj%polar%lad = 60.D0
4051! ELSE
4052! this%grid%proj%polar%lad = -60.D0
4053! ENDIF
4054! WMO says: Grid lengths are in units of metres, at the secant cone
4055! intersection parallel nearest to the pole on the projection plane.
4056 this%grid%proj%polar%lad = this%grid%proj%polar%latin1
4057 ELSE IF (editionnumber == 2) THEN
4058 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
4059 ENDIF
4060#ifdef DEBUG
4061 CALL l4f_category_log(this%category,l4f_debug, &
4062 "griddim_import_gribapi, lad "//trim(to_char(this%grid%proj%polar%lad)))
4063#endif
4064
4065! compute projected extremes from lon and lat of first point
4066 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
4067 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
4068 CALL long_reset_0_360(lofirst)
4069 CALL long_reset_to_cart_closest(this%grid%proj%lov, lofirst)
4070#ifdef DEBUG
4071 CALL l4f_category_log(this%category,l4f_debug, &
4072 "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
4073 CALL l4f_category_log(this%category,l4f_debug, &
4074 "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
4075#endif
4076
4077 CALL proj(this, lofirst, lafirst, x1, y1)
4078 IF (iscansnegatively == 0) THEN
4079 this%grid%grid%xmin = x1
4080 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
4081 ELSE
4082 this%grid%grid%xmax = x1
4083 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
4084 ENDIF
4085 IF (jscanspositively == 0) THEN
4086 this%grid%grid%ymax = y1
4087 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
4088 ELSE
4089 this%grid%grid%ymin = y1
4090 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
4091 ENDIF
4092! keep these values for personal pleasure
4093 this%grid%proj%polar%lon1 = lofirst
4094 this%grid%proj%polar%lat1 = lafirst
4095
4096! Keys for equatorial projections
4097CASE ('mercator')
4098
4099 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
4100 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
4101 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
4102 this%grid%proj%lov = 0.0d0 ! not in grib
4103
4104 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
4105 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
4106
4107 CALL proj(this, lofirst, lafirst, x1, y1)
4108 IF (iscansnegatively == 0) THEN
4109 this%grid%grid%xmin = x1
4110 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
4111 ELSE
4112 this%grid%grid%xmax = x1
4113 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
4114 ENDIF
4115 IF (jscanspositively == 0) THEN
4116 this%grid%grid%ymax = y1
4117 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
4118 ELSE
4119 this%grid%grid%ymin = y1
4120 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
4121 ENDIF
4122
4123 IF (editionnumber == 2) THEN
4124 CALL grib_get(gaid,'orientationOfTheGridInDegrees',orient)
4125 IF (orient /= 0.0d0) THEN
4126 CALL l4f_category_log(this%category,l4f_error, &
4127 "griddim_import_gribapi, Mercator grid orientation != 0 not supported")
4128 CALL raise_error()
4129 ENDIF
4130 ENDIF
4131
4132#ifdef DEBUG
4133 CALL unproj(this, x1, y1, lofirst, lafirst)
4134 CALL l4f_category_log(this%category,l4f_debug, &
4135 "griddim_import_gribapi, unprojected first point "//t2c(lofirst)//" "// &
4136 t2c(lafirst))
4137
4138 CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
4139 CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
4140 CALL proj(this, lolast, lalast, x1, y1)
4141 CALL l4f_category_log(this%category,l4f_debug, &
4142 "griddim_import_gribapi, extremes from grib "//t2c(x1)//" "//t2c(y1))
4143 CALL l4f_category_log(this%category,l4f_debug, &
4144 "griddim_import_gribapi, extremes from proj "//t2c(this%grid%grid%xmin)// &
4145 " "//t2c(this%grid%grid%ymin)//" "//t2c(this%grid%grid%xmax)//" "// &
4146 t2c(this%grid%grid%ymax))
4147#endif
4148
4149CASE ('UTM')
4150
4151 CALL grib_get(gaid,'zone',zone)
4152
4153 CALL grib_get(gaid,'datum',datum)
4154 IF (datum == 0) THEN
4155 CALL grib_get(gaid,'referenceLongitude',reflon)
4156 CALL grib_get(gaid,'falseEasting',this%grid%proj%xoff)
4157 CALL grib_get(gaid,'falseNorthing',this%grid%proj%yoff)
4158 CALL set_val(this%grid%proj, zone=zone, lov=reflon/1.0d6)
4159 ELSE
4160 CALL l4f_category_log(this%category,l4f_error,'only datum 0 supported')
4161 CALL raise_fatal_error()
4162 ENDIF
4163
4164 CALL grib_get(gaid,'eastingOfFirstGridPoint',lofirst)
4165 CALL grib_get(gaid,'eastingOfLastGridPoint',lolast)
4166 CALL grib_get(gaid,'northingOfFirstGridPoint',lafirst)
4167 CALL grib_get(gaid,'northingOfLastGridPoint',lalast)
4168
4169 IF (iscansnegatively == 0) THEN
4170 this%grid%grid%xmin = lofirst
4171 this%grid%grid%xmax = lolast
4172 ELSE
4173 this%grid%grid%xmax = lofirst
4174 this%grid%grid%xmin = lolast
4175 ENDIF
4176 IF (jscanspositively == 0) THEN
4177 this%grid%grid%ymax = lafirst
4178 this%grid%grid%ymin = lalast
4179 ELSE
4180 this%grid%grid%ymin = lafirst
4181 this%grid%grid%ymax = lalast
4182 ENDIF
4183
4184! compute dx and dy (should we get them from grib?)
4185 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
4186
4187CASE default
4188 CALL l4f_category_log(this%category,l4f_error, &
4189 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
4190 CALL raise_error()
4191
4192END SELECT
4193
4194CONTAINS
4195! utilities routines for grib_api, is there a better place?
4196SUBROUTINE grib_get_dmiss(gaid, key, value)
4197INTEGER,INTENT(in) :: gaid
4198CHARACTER(len=*),INTENT(in) :: key
4199DOUBLE PRECISION,INTENT(out) :: value
4200
4201INTEGER :: ierr
4202
4203CALL grib_get(gaid, key, value, ierr)
4204IF (ierr /= grib_success) value = dmiss
4205
4206END SUBROUTINE grib_get_dmiss
4207
4208SUBROUTINE grib_get_imiss(gaid, key, value)
4209INTEGER,INTENT(in) :: gaid
4210CHARACTER(len=*),INTENT(in) :: key
4211INTEGER,INTENT(out) :: value
4212
4213INTEGER :: ierr
4214
4215CALL grib_get(gaid, key, value, ierr)
4216IF (ierr /= grib_success) value = imiss
4217
4218END SUBROUTINE grib_get_imiss
4219
4220
4221SUBROUTINE griddim_import_ellipsoid(this, gaid)
4222TYPE(griddim_def),INTENT(inout) :: this
4223INTEGER,INTENT(in) :: gaid
4224
4225INTEGER :: shapeofearth, iv, is
4226DOUBLE PRECISION :: r1, r2
4227
4228IF (editionnumber == 2) THEN
4229 CALL grib_get(gaid, 'shapeOfTheEarth', shapeofearth)
4230 SELECT CASE(shapeofearth)
4231 CASE(0) ! spherical
4232 CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
4233 CASE(1) ! spherical generic
4234 CALL grib_get(gaid, 'scaleFactorOfRadiusOfSphericalEarth', is)
4235 CALL grib_get(gaid, 'scaledValueOfRadiusOfSphericalEarth', iv)
4236 r1 = dble(iv) / 10**is
4237 CALL set_val(this, ellips_smaj_axis=r1, ellips_flatt=0.0d0)
4238 CASE(2) ! iau65
4239 CALL set_val(this, ellips_smaj_axis=6378160.0d0, ellips_flatt=1.0d0/297.0d0)
4240 CASE(3,7) ! ellipsoidal generic
4241 CALL grib_get(gaid, 'scaleFactorOfEarthMajorAxis', is)
4242 CALL grib_get(gaid, 'scaledValueOfEarthMajorAxis', iv)
4243 r1 = dble(iv) / 10**is
4244 CALL grib_get(gaid, 'scaleFactorOfEarthMinorAxis', is)
4245 CALL grib_get(gaid, 'scaledValueOfEarthMinorAxis', iv)
4246 r2 = dble(iv) / 10**is
4247 IF (shapeofearth == 3) THEN ! km->m
4248 r1 = r1*1000.0d0
4249 r2 = r2*1000.0d0
4250 ENDIF
4251 IF (abs(r1) < 1.0d-6) THEN ! suspicious data read from grib
4252 CALL l4f_category_log(this%category,l4f_warn,'zero Earth major axis '// &
4253 'read from grib, going on with spherical Earth but the results may be wrong')
4254 CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
4255 ELSE
4256 CALL set_val(this, ellips_smaj_axis=r1, ellips_flatt=(r1-r2)/r1)
4257 ENDIF
4258 CASE(4) ! iag-grs80
4259 CALL set_val(this, ellips_type=ellips_grs80)
4260 CASE(5) ! wgs84
4261 CALL set_val(this, ellips_type=ellips_wgs84)
4262 CASE(6) ! spherical
4263 CALL set_val(this, ellips_smaj_axis=6371229.0d0, ellips_flatt=0.0d0)
4264! CASE(7) ! google earth-like?
4265 CASE default
4266 CALL l4f_category_log(this%category,l4f_error,'shapeOfTheEarth '// &
4267 t2c(shapeofearth)//' not supported in grib2')
4268 CALL raise_fatal_error()
4269
4270 END SELECT
4271
4272ELSE
4273
4274 CALL grib_get(gaid, 'earthIsOblate', shapeofearth)
4275 IF (shapeofearth == 0) THEN ! spherical
4276 CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
4277 ELSE ! iau65
4278 CALL set_val(this, ellips_smaj_axis=6378160.0d0, ellips_flatt=1.0d0/297.0d0)
4279 ENDIF
4280
4281ENDIF
4282
4283END SUBROUTINE griddim_import_ellipsoid
4284
4285
4286END SUBROUTINE griddim_import_gribapi
4287
4288
4289! grib_api driver
4290SUBROUTINE griddim_export_gribapi(this, gaid)
4291USE grib_api
4292TYPE(griddim_def),INTENT(in) :: this ! griddim object
4293INTEGER, INTENT(inout) :: gaid ! grib_api id of the grib loaded in memory to export
4294
4295INTEGER :: EditionNumber, iScansNegatively, jScansPositively, nv, pvl, zone, ierr
4296DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, reflon
4297DOUBLE PRECISION :: sdx, sdy, ratio, tol, ex, ey
4298
4299
4300! Generic keys
4301CALL grib_get(gaid,'GRIBEditionNumber', editionnumber)
4302! the following required since eccodes
4303IF (editionnumber == 2) CALL grib_set(gaid,'shapeOfTheEarth', 6)
4304CALL grib_set(gaid,'typeOfGrid', this%grid%proj%proj_type)
4305#ifdef DEBUG
4306CALL l4f_category_log(this%category,l4f_debug, &
4307 "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
4308#endif
4309
4310IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
4311 CALL grib_set(gaid,'numberOfDataPoints', this%dim%nx)
4312! reenable when possible
4313! CALL griddim_export_ellipsoid(this, gaid)
4314 RETURN
4315ENDIF
4316
4317
4318! Edition dependent setup
4319IF (editionnumber == 1) THEN
4320 ratio = 1.d3
4321ELSE IF (editionnumber == 2) THEN
4322 ratio = 1.d6
4323ELSE
4324 ratio = 0.0d0 ! signal error?!
4325ENDIF
4326
4327! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
4328CALL griddim_export_ellipsoid(this, gaid)
4329CALL grib_set(gaid, 'Ni', this%dim%nx)
4330CALL grib_set(gaid, 'Nj', this%dim%ny)
4331CALL grib_set(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
4332
4333CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
4334CALL grib_get(gaid,'jScansPositively',jscanspositively)
4335
4336! Keys for rotated grids (checked through missing values and/or error code)
4337!SELECT CASE (this%grid%proj%proj_type)
4338!CASE ('rotated_ll', 'stretched_rotated_ll', 'polar_stereographic', 'lambert', 'albers')
4339CALL grib_set_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
4340 this%grid%proj%rotated%longitude_south_pole, 0.0d0)
4341CALL grib_set_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
4342 this%grid%proj%rotated%latitude_south_pole, -90.0d0)
4343IF (editionnumber == 1) THEN
4344 CALL grib_set_dmiss(gaid,'angleOfRotationInDegrees', &
4345 this%grid%proj%rotated%angle_rotation, 0.0d0)
4346ELSE IF (editionnumber == 2)THEN
4347 CALL grib_set_dmiss(gaid,'angleOfRotationOfProjectionInDegrees', &
4348 this%grid%proj%rotated%angle_rotation, 0.0d0)
4349ENDIF
4350
4351! Keys for stretched grids (checked through missing values and/or error code)
4352! units must be verified, still experimental in grib_api
4353! # TODO: Is it a float? Is it signed?
4354IF (editionnumber == 1) THEN
4355 CALL grib_set_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
4356 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
4357 CALL grib_set_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
4358 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
4359 CALL grib_set_dmiss(gaid,'stretchingFactor', &
4360 this%grid%proj%stretched%stretch_factor, 1.0d0)
4361ELSE IF (editionnumber == 2) THEN
4362 CALL grib_set_dmiss(gaid,'longitudeOfThePoleOfStretching', &
4363 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
4364 CALL grib_set_dmiss(gaid,'latitudeOfThePoleOfStretching', &
4365 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
4366 CALL grib_set_dmiss(gaid,'stretchingFactor', &
4367 this%grid%proj%stretched%stretch_factor, 1.0d6, 1.0d6)
4368ENDIF
4369
4370! Projection-dependent keys
4371SELECT CASE (this%grid%proj%proj_type)
4372
4373! Keys for sphaerical coordinate systems
4374CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
4375
4376 IF (iscansnegatively == 0) THEN
4377 lofirst = this%grid%grid%xmin
4378 lolast = this%grid%grid%xmax
4379 ELSE
4380 lofirst = this%grid%grid%xmax
4381 lolast = this%grid%grid%xmin
4382 ENDIF
4383 IF (jscanspositively == 0) THEN
4384 lafirst = this%grid%grid%ymax
4385 lalast = this%grid%grid%ymin
4386 ELSE
4387 lafirst = this%grid%grid%ymin
4388 lalast = this%grid%grid%ymax
4389 ENDIF
4390
4391! reset lon in standard grib 2 definition [0,360]
4392 IF (editionnumber == 1) THEN
4393 CALL long_reset_m180_360(lofirst)
4394 CALL long_reset_m180_360(lolast)
4395 ELSE IF (editionnumber == 2) THEN
4396 CALL long_reset_0_360(lofirst)
4397 CALL long_reset_0_360(lolast)
4398 ENDIF
4399
4400 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
4401 CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
4402 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
4403 CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
4404
4405! test relative coordinate truncation error with respect to tol
4406! tol should be tuned
4407 sdx = this%grid%grid%dx*ratio
4408 sdy = this%grid%grid%dy*ratio
4409! error is computed relatively to the whole grid extension
4410 ex = abs(((sdx - nint(sdx))/ratio*this%dim%nx)/ &
4411 (this%grid%grid%xmax-this%grid%grid%xmin))
4412 ey = abs(((sdy - nint(sdy))/ratio*this%dim%ny)/ &
4413 (this%grid%grid%ymax-this%grid%grid%ymin))
4414 tol = 1.0d-3
4415 IF (ex > tol .OR. ey > tol) THEN
4416#ifdef DEBUG
4417 CALL l4f_category_log(this%category,l4f_debug, &
4418 "griddim_export_gribapi, error on x "//&
4419 trim(to_char(ex))//"/"//t2c(tol))
4420 CALL l4f_category_log(this%category,l4f_debug, &
4421 "griddim_export_gribapi, error on y "//&
4422 trim(to_char(ey))//"/"//t2c(tol))
4423#endif
4424! previous frmula relative to a single grid step
4425! tol = 1.0d0/ratio
4426! IF (ABS(NINT(sdx)/sdx - 1.0d0) > tol .OR. ABS(NINT(sdy)/sdy - 1.0d0) > tol) THEN
4427!#ifdef DEBUG
4428! CALL l4f_category_log(this%category,L4F_DEBUG, &
4429! "griddim_export_gribapi, dlon relative error: "//&
4430! TRIM(to_char(ABS(NINT(sdx)/sdx - 1.0d0)))//">"//TRIM(to_char(tol)))
4431! CALL l4f_category_log(this%category,L4F_DEBUG, &
4432! "griddim_export_gribapi, dlat relative error: "//&
4433! TRIM(to_char(ABS(NINT(sdy)/sdy - 1.0d0)))//">"//TRIM(to_char(tol)))
4434!#endif
4435 CALL l4f_category_log(this%category,l4f_info, &
4436 "griddim_export_gribapi, increments not given: inaccurate!")
4437 CALL grib_set_missing(gaid,'Di')
4438 CALL grib_set_missing(gaid,'Dj')
4439 CALL grib_set(gaid,'ijDirectionIncrementGiven',0)
4440 ELSE
4441#ifdef DEBUG
4442 CALL l4f_category_log(this%category,l4f_debug, &
4443 "griddim_export_gribapi, setting increments: "// &
4444 trim(to_char(this%grid%grid%dx))//' '//trim(to_char(this%grid%grid%dy)))
4445#endif
4446 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
4447 CALL grib_set(gaid,'iDirectionIncrement',nint(this%grid%grid%dx*ratio))
4448 CALL grib_set(gaid,'jDirectionIncrement',nint(this%grid%grid%dy*ratio))
4449! this does not work in grib_set
4450! CALL grib_set(gaid,'iDirectionIncrementInDegrees',this%grid%grid%dx)
4451! CALL grib_set(gaid,'jDirectionIncrementInDegrees',this%grid%grid%dy)
4452 ENDIF
4453
4454! Keys for polar projections
4455CASE ('polar_stereographic', 'lambert', 'albers')
4456! increments are required
4457 CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
4458 CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
4459 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
4460! latin1/latin2 may be missing (e.g. stereographic)
4461 CALL grib_set_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
4462 CALL grib_set_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
4463! projection center flag, aka hemisphere
4464 CALL grib_set(gaid,'projectionCenterFlag',&
4465 this%grid%proj%polar%projection_center_flag, ierr)
4466 IF (ierr /= grib_success) THEN ! try center/centre
4467 CALL grib_set(gaid,'projectionCentreFlag',&
4468 this%grid%proj%polar%projection_center_flag)
4469 ENDIF
4470
4471
4472! line of view, aka central meridian
4473 CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
4474! latitude at which dx and dy are valid
4475 IF (editionnumber == 2) THEN
4476 CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
4477 ENDIF
4478
4479! compute lon and lat of first point from projected extremes
4480 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
4481 lofirst, lafirst)
4482! reset lon in standard grib 2 definition [0,360]
4483 IF (editionnumber == 1) THEN
4484 CALL long_reset_m180_360(lofirst)
4485 ELSE IF (editionnumber == 2) THEN
4486 CALL long_reset_0_360(lofirst)
4487 ENDIF
4488 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
4489 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
4490
4491! Keys for equatorial projections
4492CASE ('mercator')
4493
4494! increments are required
4495 CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
4496 CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
4497 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
4498
4499! line of view, aka central meridian, not in grib
4500! CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
4501! latitude at which dx and dy are valid
4502 CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
4503
4504! compute lon and lat of first and last points from projected extremes
4505 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
4506 lofirst, lafirst)
4507 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
4508 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
4509 CALL get_unproj_last(this, iscansnegatively, jscanspositively, &
4510 lolast, lalast)
4511 CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
4512 CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
4513
4514 IF (editionnumber == 2) THEN
4515 CALL grib_set(gaid,'orientationOfTheGridInDegrees',0.)
4516 ENDIF
4517
4518CASE ('UTM')
4519
4520 CALL grib_set(gaid,'datum',0)
4521 CALL get_val(this, zone=zone, lov=reflon)
4522 CALL grib_set(gaid,'referenceLongitude',nint(reflon/1.0d6))
4523 CALL grib_set(gaid,'falseEasting',this%grid%proj%xoff)
4524 CALL grib_set(gaid,'falseNorthing',this%grid%proj%yoff)
4525
4526 CALL grib_set(gaid,'iDirectionIncrement',this%grid%grid%dx)
4527 CALL grib_set(gaid,'jDirectionIncrement',this%grid%grid%dy)
4528
4529!res/scann ??
4530
4531 CALL grib_set(gaid,'zone',zone)
4532
4533 IF (iscansnegatively == 0) THEN
4534 lofirst = this%grid%grid%xmin
4535 lolast = this%grid%grid%xmax
4536 ELSE
4537 lofirst = this%grid%grid%xmax
4538 lolast = this%grid%grid%xmin
4539 ENDIF
4540 IF (jscanspositively == 0) THEN
4541 lafirst = this%grid%grid%ymax
4542 lalast = this%grid%grid%ymin
4543 ELSE
4544 lafirst = this%grid%grid%ymin
4545 lalast = this%grid%grid%ymax
4546 ENDIF
4547
4548 CALL grib_set(gaid,'eastingOfFirstGridPoint',lofirst)
4549 CALL grib_set(gaid,'eastingOfLastGridPoint',lolast)
4550 CALL grib_set(gaid,'northingOfFirstGridPoint',lafirst)
4551 CALL grib_set(gaid,'northingOfLastGridPoint',lalast)
4552
4553CASE default
4554 CALL l4f_category_log(this%category,l4f_error, &
4555 "griddim_export_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
4556 CALL raise_error()
4557
4558END SELECT
4559
4560! hack for position of vertical coordinate parameters
4561! buggy in grib_api
4562IF (editionnumber == 1) THEN
4563! CALL grib_get(gaid,"PVPresent",pvp) ! alias, probably useless
4564 CALL grib_get(gaid,"NV",nv)
4565#ifdef DEBUG
4566 CALL l4f_category_log(this%category,l4f_debug,"griddim_export_gribapi, coding "// &
4567 trim(to_char(nv))//" vertical coordinate parameters")
4568#endif
4569
4570 IF (nv == 0) THEN
4571 pvl = 255
4572 ELSE
4573 SELECT CASE (this%grid%proj%proj_type)
4574 CASE ('regular_ll') ! check whether "29-32 Set to zero (reserved)" are required
4575 pvl = 33
4576 CASE ('polar_stereographic')
4577 pvl = 33
4578 CASE ('rotated_ll', 'stretched_ll', 'lambert', 'albers')
4579 pvl = 43
4580 CASE ('stretched_rotated_ll')
4581 pvl = 43
4582 CASE DEFAULT
4583 pvl = 43 !?
4584 END SELECT
4585 ENDIF
4586
4587 CALL grib_set(gaid,"pvlLocation",pvl)
4588#ifdef DEBUG
4589 CALL l4f_category_log(this%category,l4f_debug,"griddim_export_gribapi, coding "// &
4590 trim(to_char(pvl))//" as vertical coordinate parameter location")
4591#endif
4592
4593ENDIF
4594
4595
4596CONTAINS
4597! utilities routines for grib_api, is there a better place?
4598SUBROUTINE grib_set_dmiss(gaid, key, val, default, factor)
4599INTEGER,INTENT(in) :: gaid
4600CHARACTER(len=*),INTENT(in) :: key
4601DOUBLE PRECISION,INTENT(in) :: val
4602DOUBLE PRECISION,INTENT(in),OPTIONAL :: default
4603DOUBLE PRECISION,INTENT(in),OPTIONAL :: factor
4604
4605INTEGER :: ierr
4606
4607IF (c_e(val)) THEN
4608 IF (PRESENT(factor)) THEN
4609 CALL grib_set(gaid, key, val*factor, ierr)
4610 ELSE
4611 CALL grib_set(gaid, key, val, ierr)
4612 ENDIF
4613ELSE IF (PRESENT(default)) THEN
4614 CALL grib_set(gaid, key, default, ierr)
4615ENDIF
4616
4617END SUBROUTINE grib_set_dmiss
4618
4619SUBROUTINE grib_set_imiss(gaid, key, value, default)
4620INTEGER,INTENT(in) :: gaid
4621CHARACTER(len=*),INTENT(in) :: key
4622INTEGER,INTENT(in) :: value
4623INTEGER,INTENT(in),OPTIONAL :: default
4624
4625INTEGER :: ierr
4626
4627IF (c_e(value)) THEN
4628 CALL grib_set(gaid, key, value, ierr)
4629ELSE IF (PRESENT(default)) THEN
4630 CALL grib_set(gaid, key, default, ierr)
4631ENDIF
4632
4633END SUBROUTINE grib_set_imiss
4634
4635SUBROUTINE griddim_export_ellipsoid(this, gaid)
4636TYPE(griddim_def),INTENT(in) :: this
4637INTEGER,INTENT(in) :: gaid
4638
4639INTEGER :: ellips_type, ierr
4640DOUBLE PRECISION :: r1, r2, f
4641
4642CALL get_val(this, ellips_smaj_axis=r1, ellips_flatt=f, ellips_type=ellips_type)
4643
4644IF (editionnumber == 2) THEN
4645
4646! why does it produce a message even with ierr?
4647 CALL grib_set_missing(gaid, 'scaleFactorOfRadiusOfSphericalEarth', ierr)
4648 CALL grib_set_missing(gaid, 'scaledValueOfRadiusOfSphericalEarth', ierr)
4649 CALL grib_set_missing(gaid, 'scaleFactorOfEarthMajorAxis', ierr)
4650 CALL grib_set_missing(gaid, 'scaledValueOfEarthMajorAxis', ierr)
4651 CALL grib_set_missing(gaid, 'scaleFactorOfEarthMinorAxis', ierr)
4652 CALL grib_set_missing(gaid, 'scaledValueOfEarthMinorAxis', ierr)
4653
4654 SELECT CASE(ellips_type)
4655 CASE(ellips_grs80) ! iag-grs80
4656 CALL grib_set(gaid, 'shapeOfTheEarth', 4)
4657 CASE(ellips_wgs84) ! wgs84
4658 CALL grib_set(gaid, 'shapeOfTheEarth', 5)
4659 CASE default
4660 IF (f == 0.0d0) THEN ! spherical Earth
4661 IF (r1 == 6367470.0d0) THEN ! spherical
4662 CALL grib_set(gaid, 'shapeOfTheEarth', 0)
4663 ELSE IF (r1 == 6371229.0d0) THEN ! spherical
4664 CALL grib_set(gaid, 'shapeOfTheEarth', 6)
4665 ELSE ! spherical generic
4666 CALL grib_set(gaid, 'shapeOfTheEarth', 1)
4667 CALL grib_set(gaid, 'scaleFactorOfRadiusOfSphericalEarth', 2)
4668 CALL grib_set(gaid, 'scaledValueOfRadiusOfSphericalEarth', int(r1*100.0d0))
4669 ENDIF
4670 ELSE ! ellipsoidal
4671 IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
4672 CALL grib_set(gaid, 'shapeOfTheEarth', 2)
4673 ELSE ! ellipsoidal generic
4674 CALL grib_set(gaid, 'shapeOfTheEarth', 3)
4675 r2 = r1*(1.0d0 - f)
4676 CALL grib_set(gaid, 'scaleFactorOfEarthMajorAxis', 5)
4677 CALL grib_set(gaid, 'scaledValueOfEarthMajorAxis', &
4678 int(r1*100.0d0))
4679 CALL grib_set(gaid, 'scaleFactorOfEarthMinorAxis', 5)
4680 CALL grib_set(gaid, 'scaledValueOfEarthMinorAxis', &
4681 int(r2*100.0d0))
4682 ENDIF
4683 ENDIF
4684 END SELECT
4685
4686ELSE
4687
4688 IF (r1 == 6367470.0d0 .AND. f == 0.0d0) THEN ! spherical
4689 CALL grib_set(gaid, 'earthIsOblate', 0)
4690 ELSE IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
4691 CALL grib_set(gaid, 'earthIsOblate', 1)
4692 ELSE IF (f == 0.0d0) THEN ! generic spherical
4693 CALL grib_set(gaid, 'earthIsOblate', 0)
4694 CALL l4f_category_log(this%category,l4f_warn,'desired spherical Earth radius &
4695 &not supported in grib 1, coding with default radius of 6367470 m')
4696 ELSE ! generic ellipsoidal
4697 CALL grib_set(gaid, 'earthIsOblate', 1)
4698 CALL l4f_category_log(this%category,l4f_warn,'desired Earth ellipsoid &
4699 &not supported in grib 1, coding with default iau65 ellipsoid')
4700 ENDIF
4701
4702ENDIF
4703
4704END SUBROUTINE griddim_export_ellipsoid
4705
4706SUBROUTINE get_unproj_first(this, iScansNegatively, jScansPositively, &
4707 loFirst, laFirst)
4708TYPE(griddim_def),INTENT(in) :: this ! griddim object
4709INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
4710DOUBLE PRECISION,INTENT(out) :: loFirst, laFirst
4711
4712! compute lon and lat of first point from projected extremes
4713IF (iscansnegatively == 0) THEN
4714 IF (jscanspositively == 0) THEN
4715 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymax, lofirst, lafirst)
4716 ELSE
4717 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymin, lofirst, lafirst)
4718 ENDIF
4719ELSE
4720 IF (jscanspositively == 0) THEN
4721 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymax, lofirst, lafirst)
4722 ELSE
4723 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymin, lofirst, lafirst)
4724 ENDIF
4725ENDIF
4726! use the values kept for personal pleasure ?
4727! loFirst = this%grid%proj%polar%lon1
4728! laFirst = this%grid%proj%polar%lat1
4729END SUBROUTINE get_unproj_first
4730
4731SUBROUTINE get_unproj_last(this, iScansNegatively, jScansPositively, &
4732 loLast, laLast)
4733TYPE(griddim_def),INTENT(in) :: this ! griddim object
4734INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
4735DOUBLE PRECISION,INTENT(out) :: loLast, laLast
4736
4737! compute lon and lat of last point from projected extremes
4738IF (iscansnegatively == 0) THEN
4739 IF (jscanspositively == 0) THEN
4740 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymin, lolast, lalast)
4741 ELSE
4742 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymax, lolast, lalast)
4743 ENDIF
4744ELSE
4745 IF (jscanspositively == 0) THEN
4746 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymin, lolast, lalast)
4747 ELSE
4748 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymax, lolast, lalast)
4749 ENDIF
4750ENDIF
4751! use the values kept for personal pleasure ?
4752! loLast = this%grid%proj%polar%lon?
4753! laLast = this%grid%proj%polar%lat?
4754END SUBROUTINE get_unproj_last
4755
4756END SUBROUTINE griddim_export_gribapi
4757#endif
4758
4759
4760#ifdef HAVE_LIBGDAL
4761! gdal driver
4762SUBROUTINE griddim_import_gdal(this, gdalid, gdal_options)
4763USE gdal
4764TYPE(griddim_def),INTENT(inout) :: this ! griddim object
4765TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal rasterband pointer
4766TYPE(gdal_file_id_options),INTENT(in) :: gdal_options
4767
4768TYPE(gdaldataseth) :: hds
4769REAL(kind=c_double) :: geotrans(6) !, invgeotrans(6), x1, y1, x2, y2
4770INTEGER(kind=c_int) :: offsetx, offsety
4771INTEGER :: ier
4772
4773hds = gdalgetbanddataset(gdalid) ! go back to dataset
4774ier = gdalgetgeotransform(hds, geotrans)
4775
4776IF (ier /= 0) THEN
4777 CALL l4f_category_log(this%category, l4f_error, &
4778 'griddim_import_gdal: error in accessing gdal dataset')
4779 CALL raise_error()
4780 RETURN
4781ENDIF
4782IF (geotrans(3) /= 0.0_c_double .OR. geotrans(5) /= 0.0_c_double) THEN ! dataset has not a diagonal transformation
4783 CALL l4f_category_log(this%category, l4f_error, &
4784 'griddim_import_gdal: dataset has a non-diagonal transformation matrix, unsupported')
4785 CALL raise_error()
4786 RETURN
4787ENDIF
4788
4789CALL gdaldatasetbbsize_f(hds, gdal_options%xmin, gdal_options%ymin, &
4790 gdal_options%xmax, gdal_options%ymax, &
4791 this%dim%nx, this%dim%ny, offsetx, offsety, &
4792 this%grid%grid%xmin, this%grid%grid%ymin, this%grid%grid%xmax, this%grid%grid%ymax)
4793
4794IF (this%dim%nx == 0 .OR. this%dim%ny == 0) THEN
4795 CALL l4f_category_log(this%category, l4f_warn, &
4796 'griddim_import_gdal: requested bounding box '//t2c(gdal_options%xmin)//','// &
4797 t2c(gdal_options%ymin)//','//t2c(gdal_options%xmax)//','//&
4798 t2c(gdal_options%ymax))
4799 CALL l4f_category_log(this%category, l4f_warn, &
4800 'determines an empty dataset '//t2c(this%dim%nx)//'x'//t2c(this%dim%ny))
4801ENDIF
4802
4803! get grid corners
4804!CALL gdalapplygeotransform(geotrans, 0.5_c_double, 0.5_c_double, x1, y1)
4805!CALL gdalapplygeotransform(geotrans, gdalgetrasterbandxsize(gdalid)-0.5_c_double, &
4806! gdalgetrasterbandysize(gdalid)-0.5_c_double, x2, y2)
4807
4808!IF (geotrans(3) == 0.0_c_double .AND. geotrans(5) == 0.0_c_double) THEN ! transformation is diagonal, no transposing
4809! this%dim%nx = gdalgetrasterbandxsize(gdalid)
4810! this%dim%ny = gdalgetrasterbandysize(gdalid)
4811! this%grid%grid%xmin = MIN(x1, x2)
4812! this%grid%grid%xmax = MAX(x1, x2)
4813! this%grid%grid%ymin = MIN(y1, y2)
4814! this%grid%grid%ymax = MAX(y1, y2)
4815!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
4816!
4817! this%dim%nx = gdalgetrasterbandysize(gdalid)
4818! this%dim%ny = gdalgetrasterbandxsize(gdalid)
4819! this%grid%grid%xmin = MIN(y1, y2)
4820! this%grid%grid%xmax = MAX(y1, y2)
4821! this%grid%grid%ymin = MIN(x1, x2)
4822! this%grid%grid%ymax = MAX(x1, x2)
4823!
4824!ELSE ! transformation is a rotation, not supported
4825!ENDIF
4826
4827this%grid%proj%proj_type = 'regular_ll' ! forced, only one supported (check?)
4828
4829CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
4830this%grid%grid%component_flag = 0
4831
4832END SUBROUTINE griddim_import_gdal
4833#endif
4834
4835
4836!> Display on the screen a brief content of the object.
4837SUBROUTINE griddim_display(this)
4838TYPE(griddim_def),INTENT(in) :: this !< griddim object to display
4839
4840print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
4841
4842CALL display(this%grid%proj)
4843CALL display(this%grid%grid)
4844CALL display(this%dim)
4845
4846print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
4847
4848END SUBROUTINE griddim_display
4849
4850
4851#define VOL7D_POLY_TYPE TYPE(grid_def)
4852#define VOL7D_POLY_TYPES _grid
4853#include "array_utilities_inc.F90"
4854#undef VOL7D_POLY_TYPE
4855#undef VOL7D_POLY_TYPES
4856
4857#define VOL7D_POLY_TYPE TYPE(griddim_def)
4858#define VOL7D_POLY_TYPES _griddim
4859#include "array_utilities_inc.F90"
4860#undef VOL7D_POLY_TYPE
4861#undef VOL7D_POLY_TYPES
4862
4863
4864!> Compute rotation matrix for wind unrotation. It allocates and
4865!! computes a matrix suitable for recomputing wind components in the
4866!! geographical system from the components in the projected
4867!! system. The rotation matrix \a rot_mat has to be passed as a
4868!! pointer and successively deallocated by the caller; it is a
4869!! 3-dimensional array where the first two dimensions are lon and lat
4870!! and the third, with extension 4, contains the packed rotation
4871!! matrix for the given grid point. It should work for every
4872!! projection. In order for the method to work, the \a griddim_unproj
4873!! method must have already been called for the \a griddim_def object.
4874!! \todo Check the algorithm and add some orthogonality tests.
4875SUBROUTINE griddim_wind_unrot(this, rot_mat)
4877TYPE(griddim_def), INTENT(in) :: this !< object describing the grid
4878DOUBLE PRECISION, POINTER :: rot_mat(:,:,:) !< rotation matrix for every grid point, to be deallocated by the caller, if \c .NOT. \c ASSOCIATED() an error occurred
4879
4880DOUBLE PRECISION :: dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j
4881DOUBLE PRECISION :: lat_factor
4882INTEGER :: i, j, ip1, im1, jp1, jm1;
4883
4884IF (this%dim%nx <= 0 .OR. this%dim%ny <= 0 .OR. &
4885 .NOT. ASSOCIATED(this%dim%lon) .OR. .NOT. ASSOCIATED(this%dim%lat)) THEN
4886 CALL l4f_category_log(this%category, l4f_error, 'rotate_uv must be called after correct init of griddim object')
4887 NULLIFY(rot_mat)
4888 RETURN
4889ENDIF
4890
4891ALLOCATE(rot_mat(this%dim%nx, this%dim%ny, 4))
4892
4893DO j = 1, this%dim%ny
4894 jp1 = min(j+1, this%dim%ny)
4895 jm1 = max(j-1, 1)
4896 DO i = 1, this%dim%nx
4897 ip1 = min(i+1, this%dim%nx)
4898 im1 = max(i-1, 1)
4899
4900 dlat_i = this%dim%lat(ip1,j) - this%dim%lat(im1,j)
4901 dlat_j = this%dim%lat(i,jp1) - this%dim%lat(i,jm1)
4902
4903 dlon_i = this%dim%lon(ip1,j) - this%dim%lon(im1,j)
4904! if ( dlon_i > pi ) dlon_i -= 2*pi;
4905! if ( dlon_i < (-pi) ) dlon_i += 2*pi;
4906 dlon_j = this%dim%lon(i,jp1) - this%dim%lon(i,jm1)
4907! if ( dlon_j > pi ) dlon_j -= 2*pi;
4908! if ( dlon_j < (-pi) ) dlon_j += 2*pi;
4909
4910! check whether this is really necessary !!!!
4911 lat_factor = cos(degrad*this%dim%lat(i,j))
4912 dlon_i = dlon_i * lat_factor
4913 dlon_j = dlon_j * lat_factor
4914
4915 dist_i = sqrt(dlon_i*dlon_i + dlat_i*dlat_i)
4916 dist_j = sqrt(dlon_j*dlon_j + dlat_j*dlat_j)
4917
4918 IF (dist_i > 0.d0) THEN
4919 rot_mat(i,j,1) = dlon_i/dist_i
4920 rot_mat(i,j,2) = dlat_i/dist_i
4921 ELSE
4922 rot_mat(i,j,1) = 0.0d0
4923 rot_mat(i,j,2) = 0.0d0
4924 ENDIF
4925 IF (dist_j > 0.d0) THEN
4926 rot_mat(i,j,3) = dlon_j/dist_j
4927 rot_mat(i,j,4) = dlat_j/dist_j
4928 ELSE
4929 rot_mat(i,j,3) = 0.0d0
4930 rot_mat(i,j,4) = 0.0d0
4931 ENDIF
4932
4933 ENDDO
4934ENDDO
4935
4936END SUBROUTINE griddim_wind_unrot
4937
4938
4939! compute zoom indices from geographical zoom coordinates
4940SUBROUTINE griddim_zoom_coord(this, ilon, ilat, flon, flat, ix, iy, fx, fy)
4941TYPE(griddim_def),INTENT(in) :: this
4942DOUBLE PRECISION,INTENT(in) :: ilon, ilat, flon, flat
4943INTEGER,INTENT(out) :: ix, iy, fx, fy
4944
4945DOUBLE PRECISION :: ix1, iy1, fx1, fy1
4946
4947! compute projected coordinates of vertices of desired lonlat rectangle
4948CALL proj(this, ilon, ilat, ix1, iy1)
4949CALL proj(this, flon, flat, fx1, fy1)
4950
4951CALL griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
4952
4953END SUBROUTINE griddim_zoom_coord
4954
4955
4956! compute zoom indices from projected zoom coordinates
4957SUBROUTINE griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
4958TYPE(griddim_def),INTENT(in) :: this
4959DOUBLE PRECISION,INTENT(in) :: ix1, iy1, fx1, fy1
4960INTEGER,INTENT(out) :: ix, iy, fx, fy
4961
4962INTEGER :: lix, liy, lfx, lfy
4963
4964! compute projected indices
4965lix = nint((ix1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
4966liy = nint((iy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
4967lfx = nint((fx1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
4968lfy = nint((fy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
4969! swap projected indices, in case grid is "strongly rotated"
4970ix = min(lix, lfx)
4971fx = max(lix, lfx)
4972iy = min(liy, lfy)
4973fy = max(liy, lfy)
4974
4975END SUBROUTINE griddim_zoom_projcoord
4976
4977
4978!> Reset a longitude value in the interval [0-360[.
4979!! The value is reset in place. This is usually useful in connection
4980!! with grib2 coding/decoding.
4981SUBROUTINE long_reset_0_360(lon)
4982DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
4983
4984IF (.NOT.c_e(lon)) RETURN
4985DO WHILE(lon < 0.0d0)
4986 lon = lon + 360.0d0
4987END DO
4988DO WHILE(lon >= 360.0d0)
4989 lon = lon - 360.0d0
4990END DO
4991
4992END SUBROUTINE long_reset_0_360
4993
4994
4995!> Reset a longitude value in the interval [-180-360[.
4996!! The value is reset in place. This is usually useful in connection
4997!! with grib1 coding/decoding.
4998SUBROUTINE long_reset_m180_360(lon)
4999DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
5000
5001IF (.NOT.c_e(lon)) RETURN
5002DO WHILE(lon < -180.0d0)
5003 lon = lon + 360.0d0
5004END DO
5005DO WHILE(lon >= 360.0d0)
5006 lon = lon - 360.0d0
5007END DO
5008
5009END SUBROUTINE long_reset_m180_360
5010
5011
5012!> Reset a longitude value in the interval [-90-270[.
5013!! The value is reset in place. This is usually useful in connection
5014!! with grib2 coding/decoding.
5015!SUBROUTINE long_reset_m90_270(lon)
5016!DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
5017!
5018!IF (.NOT.c_e(lon)) RETURN
5019!DO WHILE(lon < -90.0D0)
5020! lon = lon + 360.0D0
5021!END DO
5022!DO WHILE(lon >= 270.0D0)
5023! lon = lon - 360.0D0
5024!END DO
5025!
5026!END SUBROUTINE long_reset_m90_270
5027
5028
5029!> Reset a longitude value in the interval [-180-180[.
5030!! The value is reset in place. This is usually useful in connection
5031!! with grib2 coding/decoding.
5032SUBROUTINE long_reset_m180_180(lon)
5033DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
5034
5035IF (.NOT.c_e(lon)) RETURN
5036DO WHILE(lon < -180.0d0)
5037 lon = lon + 360.0d0
5038END DO
5039DO WHILE(lon >= 180.0d0)
5040 lon = lon - 360.0d0
5041END DO
5042
5043END SUBROUTINE long_reset_m180_180
5044
5045
5046SUBROUTINE long_reset_to_cart_closest(lon, lonref)
5047DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
5048DOUBLE PRECISION,INTENT(in) :: lonref !< the longitude to compare
5049
5050IF (.NOT.c_e(lon) .OR. .NOT.c_e(lonref)) RETURN
5051IF (abs(lon-lonref) < 180.0d0) RETURN ! nothing to do
5052lon = lon - nint((lon-lonref)/360.0d0)*360.0d0 ! se non e` vera e` ben trovata
5053
5054END SUBROUTINE long_reset_to_cart_closest
5055
5056
5057END MODULE grid_class
5058
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.

Generated with Doxygen.