libsim Versione 7.2.6

◆ map_inv_distinct_griddim()

integer function, dimension(dim) map_inv_distinct_griddim ( type(griddim_def), dimension(:), intent(in) vect,
integer, intent(in) dim,
logical, dimension(:), intent(in), optional mask,
logical, intent(in), optional back )
private

map inv distinct

Definizione alla linea 2817 del file grid_class.F90.

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