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