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