libsim Versione 7.2.6

◆ index_griddim()

integer function index_griddim ( type(griddim_def), dimension(:), intent(in) vect,
type(griddim_def), intent(in) search,
logical, dimension(:), intent(in), optional mask,
logical, intent(in), optional back,
integer, intent(in), optional cache )
private

Cerca l'indice del primo o ultimo elemento di vect uguale a search.

Definizione alla linea 2903 del file grid_class.F90.

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

Generated with Doxygen.