libsim Versione 7.2.6

◆ map_distinct_grid()

integer function, dimension(size(vect)) map_distinct_grid ( type(grid_def), dimension(:), intent(in) vect,
logical, dimension(:), intent(in), optional mask,
logical, intent(in), optional back )
private

map distinct

Definizione alla linea 2207 del file grid_class.F90.

2208! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2209! authors:
2210! Davide Cesari <dcesari@arpa.emr.it>
2211! Paolo Patruno <ppatruno@arpa.emr.it>
2212
2213! This program is free software; you can redistribute it and/or
2214! modify it under the terms of the GNU General Public License as
2215! published by the Free Software Foundation; either version 2 of
2216! the License, or (at your option) any later version.
2217
2218! This program is distributed in the hope that it will be useful,
2219! but WITHOUT ANY WARRANTY; without even the implied warranty of
2220! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2221! GNU General Public License for more details.
2222
2223! You should have received a copy of the GNU General Public License
2224! along with this program. If not, see <http://www.gnu.org/licenses/>.
2225#include "config.h"
2226!> Module for describing geographically referenced regular grids.
2227!! This module defines classes and methods describing rectangular
2228!! georeferenced grids in different geographical projections. The grid
2229!! and projection definition can be specified explicitely by the
2230!! caller or they can be entirely imported from a grib file (through
2231!! grib_api) or from a file format supported by gdal (through gdal
2232!! fortran interface).
2233!!
2234!! The projection is internally stored following the WMO grib
2235!! conventions (gridType in grib_api). The projections currently
2236!! supported or for which support is planned are:
2237!! - For grib edition 1 and 2:
2238!! - regular_ll (works)
2239!! - mercator (to be done)
2240!! - lambert (works, to be completed)
2241!! - polar_stereographic (to be tested)
2242!! - albers (to be done)
2243!! - rotated_ll (works)
2244!! - stretched_ll (to be completed and tested)
2245!! - stretched_rotated_ll (to be completed and tested)
2246!! - For grib edition 2 only:
2247!! - UTM (ARPA-SIM extension)
2248!! - equatorial_azimuthal_equidistant (to be done)
2249!! - For gdal-supported formats:
2250!! - regular_ll
2251!!
2252!!
2253!! See the example program \include example_vg6d_1.f90
2254!!
2255!!\ingroup volgrid6d
2256MODULE grid_class
2257use geo_proj_class
2259use grid_rect_class
2260use grid_id_class
2261use err_handling
2264use log4fortran
2265implicit none
2266
2267
2268character (len=255),parameter:: subcategory="grid_class"
2269
2270
2271!> This object, mainly for internal use, describes a grid on
2272!! a geographical projection, except the grid dimensions.
2273!! The object is opaque, thus all its members have to be set and
2274!! accessed through the constructor and the ::get_val and ::set_val
2275!! methods.
2276type grid_def
2277 private
2278 type(geo_proj) :: proj
2279 type(grid_rect) :: grid
2280 integer :: category = 0
2281end type grid_def
2282
2283
2284!> This object completely describes a grid on a geographic projection.
2285!! It is the main public object of this module. The grid definition \a
2286!! grid is separated from the definition of the grid dimensions \a dim
2287!! in order to make members of \a grid \a PRIVATE while maintaining
2288!! free access to the members of \a dim.
2289type griddim_def
2290 type(grid_def) :: grid !< grid and projection definition
2291 type(grid_dim) :: dim !< grid dimensions definition
2292 integer :: category = 0 !< category for log4fortran
2293end type griddim_def
2294
2295
2296!> Logical equality operators for objects of the classes \a grid_def,
2297!! and \a griddim_def. They are all defined as \c
2298!! ELEMENTAL thus work also on arrays of any shape.
2299INTERFACE OPERATOR (==)
2300 MODULE PROCEDURE grid_eq, griddim_eq
2301END INTERFACE
2302
2303!> Logical inequality operators for objects of the classes \a grid_def,
2304!! and \a griddim_def. They are all defined as \c
2305!! ELEMENTAL thus work also on arrays of any shape.
2306INTERFACE OPERATOR (/=)
2307 MODULE PROCEDURE grid_ne, griddim_ne
2308END INTERFACE
2309
2310!> Constructors of the corresponding objects.
2311INTERFACE init
2312 MODULE PROCEDURE griddim_init
2313END INTERFACE
2314
2315!> Destructors of the corresponding objects.
2316INTERFACE delete
2317 MODULE PROCEDURE griddim_delete
2318END INTERFACE
2319
2320!> Copy an object, creating a fully new instance.
2321INTERFACE copy
2322 MODULE PROCEDURE griddim_copy
2323END INTERFACE
2324
2325!> Compute forward coordinate transformation from geographical system to
2326!! projected system.
2327INTERFACE proj
2328 MODULE PROCEDURE griddim_coord_proj!, griddim_proj
2329END INTERFACE
2330
2331!> Compute backward coordinate transformation from projected system
2332!! to geographical system.
2333INTERFACE unproj
2334 MODULE PROCEDURE griddim_coord_unproj, griddim_unproj
2335END INTERFACE
2336
2337!> Method for returning the contents of the object.
2338INTERFACE get_val
2339 MODULE PROCEDURE griddim_get_val
2340END INTERFACE
2341
2342!> Method for setting the contents of the object.
2343INTERFACE set_val
2344 MODULE PROCEDURE griddim_set_val
2345END INTERFACE
2346
2347!> Write the object on a formatted or unformatted file.
2348INTERFACE write_unit
2349 MODULE PROCEDURE griddim_write_unit
2350END INTERFACE
2351
2352!> Read the object from a formatted or unformatted file.
2353INTERFACE read_unit
2354 MODULE PROCEDURE griddim_read_unit
2355END INTERFACE
2356
2357!> Import griddim object from grid_id.
2358INTERFACE import
2359 MODULE PROCEDURE griddim_import_grid_id
2360END INTERFACE
2361
2362!> Export griddim object to grid_id.
2363INTERFACE export
2364 MODULE PROCEDURE griddim_export_grid_id
2365END INTERFACE
2366
2367!> Print a brief description on stdout.
2368INTERFACE display
2369 MODULE PROCEDURE griddim_display
2370END INTERFACE
2371
2372#define VOL7D_POLY_TYPE TYPE(grid_def)
2373#define VOL7D_POLY_TYPES _grid
2374#include "array_utilities_pre.F90"
2375#undef VOL7D_POLY_TYPE
2376#undef VOL7D_POLY_TYPES
2377
2378#define VOL7D_POLY_TYPE TYPE(griddim_def)
2379#define VOL7D_POLY_TYPES _griddim
2380#include "array_utilities_pre.F90"
2381#undef VOL7D_POLY_TYPE
2382#undef VOL7D_POLY_TYPES
2383
2384INTERFACE wind_unrot
2385 MODULE PROCEDURE griddim_wind_unrot
2386END INTERFACE
2387
2388
2389PRIVATE
2390
2391PUBLIC proj, unproj, griddim_unproj, griddim_gen_coord, &
2392 griddim_zoom_coord, griddim_zoom_projcoord, &
2393 griddim_setsteps, griddim_def, grid_def, grid_dim
2394PUBLIC init, delete, copy
2396PUBLIC OPERATOR(==),OPERATOR(/=)
2397PUBLIC count_distinct, pack_distinct, count_and_pack_distinct, &
2398 map_distinct, map_inv_distinct,index
2399PUBLIC wind_unrot, import, export
2400PUBLIC griddim_central_lon, griddim_set_central_lon
2401CONTAINS
2402
2403!> Constructor for a \a griddim_def object.
2404SUBROUTINE griddim_init(this, nx, ny, &
2405 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
2406 proj_type, lov, zone, xoff, yoff, &
2407 longitude_south_pole, latitude_south_pole, angle_rotation, &
2408 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
2409 latin1, latin2, lad, projection_center_flag, &
2410 ellips_smaj_axis, ellips_flatt, ellips_type, &
2411 categoryappend)
2412TYPE(griddim_def),INTENT(inout) :: this !< object to be created
2413INTEGER,INTENT(in),OPTIONAL :: nx !< number of points along the x axis
2414INTEGER,INTENT(in),OPTIONAL :: ny !< number of points along the y axis
2415DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin !< lower bound for x coordinate on grid in projection units (degrees or meters depending on the projection type)
2416DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmax !< upper bound for x coordinate
2417DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymin !< lower bound for y coordinate
2418DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymax !< upper bound for y coordinate
2419DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx !< grid step in x direction
2420DOUBLE PRECISION,INTENT(in),OPTIONAL :: dy !< grid step in y direction
2421!> Resolved u- and v- components of vector quantities relative to 0=the easterly and northerly directions
2422!! 1=the defined grid in the direction of increasing x and y (or i and j) coordinates respectively (0=north, 128=south)
2423INTEGER,INTENT(in),OPTIONAL :: component_flag
2424CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type !< type of projection
2425DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov !< line of view, also known as reference longitude or orientation of the grid (polar projections)
2426INTEGER,INTENT(in),OPTIONAL :: zone !< Earth zone (mainly for UTM), sets lov to the correct zone central meridian
2427DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff !< offset on x axis (false easting)
2428DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff !< offset on y axis (false northing)
2429DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole !< longitude of the southern pole of projection
2430DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole !< latitude of the southern pole of projection
2431DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation !< angle of rotation of projection
2432DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole !< longitude of the pole of stretching
2433DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole !< latitude of the pole of stretching
2434DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor !< stretching factor
2435DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1 !< first standard latitude from main pole (Lambert)
2436DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2 !< second standard latitude from main pole (Lambert)
2437DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad !< latitude at which dx and dy (in m) are specified (Lambert, grib2 only)
2438INTEGER,INTENT(in),OPTIONAL :: projection_center_flag !< flag indicating which pole is represented
2439DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis !< Earth semi-major axis
2440DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt !< Earth flattening
2441INTEGER,INTENT(in),OPTIONAL :: ellips_type !< number in the interval [1,nellips] indicating a predefined ellipsoid, alternative to the previous arguments
2442CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
2443
2444CHARACTER(len=512) :: a_name
2445
2446IF (PRESENT(categoryappend)) THEN
2447 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
2448 trim(categoryappend))
2449ELSE
2450 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2451ENDIF
2452this%category=l4f_category_get(a_name)
2453
2454! geographical projection
2455this%grid%proj = geo_proj_new( &
2456 proj_type=proj_type, lov=lov, zone=zone, xoff=xoff, yoff=yoff, &
2457 longitude_south_pole=longitude_south_pole, &
2458 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
2459 longitude_stretch_pole=longitude_stretch_pole, &
2460 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
2461 latin1=latin1, latin2=latin2, lad=lad, projection_center_flag=projection_center_flag, &
2462 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, ellips_type=ellips_type)
2463! grid extension
2464this%grid%grid = grid_rect_new( &
2465 xmin, xmax, ymin, ymax, dx, dy, component_flag)
2466! grid size
2467this%dim = grid_dim_new(nx, ny)
2468
2469#ifdef DEBUG
2470call l4f_category_log(this%category,l4f_debug,"init gtype: "//this%grid%proj%proj_type )
2471#endif
2472
2473END SUBROUTINE griddim_init
2474
2475
2476!> Destroy a \a griddim_def object.
2477SUBROUTINE griddim_delete(this)
2478TYPE(griddim_def),INTENT(inout) :: this !< object to be destroyed
2479
2480CALL delete(this%dim)
2481CALL delete(this%grid%proj)
2482CALL delete(this%grid%grid)
2483
2484call l4f_category_delete(this%category)
2485
2486END SUBROUTINE griddim_delete
2487
2488
2489!> Create an independent copy of a \a griddim_def object.
2490SUBROUTINE griddim_copy(this, that, categoryappend)
2491TYPE(griddim_def),INTENT(in) :: this !< object to be copied
2492TYPE(griddim_def),INTENT(out) :: that !< copied object
2493CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
2494
2495CHARACTER(len=512) :: a_name
2496
2497CALL init(that)
2498
2499CALL copy(this%grid%proj, that%grid%proj)
2500CALL copy(this%grid%grid, that%grid%grid)
2501CALL copy(this%dim, that%dim)
2502
2503! new category
2504IF (PRESENT(categoryappend)) THEN
2505 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
2506 trim(categoryappend))
2507ELSE
2508 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2509ENDIF
2510that%category=l4f_category_get(a_name)
2511
2512END SUBROUTINE griddim_copy
2513
2514
2515!> Computes and returns coordinates in the projected system given the
2516!! geographical coordinates.
2517ELEMENTAL SUBROUTINE griddim_coord_proj(this, lon, lat, x, y)
2518TYPE(griddim_def),INTENT(in) :: this !< definition of projection
2519!> geographical coordinates
2520DOUBLE PRECISION,INTENT(in) :: lon, lat
2521!> projected coordinates
2522DOUBLE PRECISION,INTENT(out) :: x, y
2523
2524CALL proj(this%grid%proj, lon, lat, x, y)
2525
2526END SUBROUTINE griddim_coord_proj
2527
2528
2529!> Computes and returns geographical coordinates given the coordinates
2530!! in the projected system.
2531ELEMENTAL SUBROUTINE griddim_coord_unproj(this, x, y, lon, lat)
2532TYPE(griddim_def),INTENT(in) :: this !< definition of projection
2533!> projected coordinates
2534DOUBLE PRECISION,INTENT(in) :: x, y
2535!> geographical coordinates
2536DOUBLE PRECISION,INTENT(out) :: lon, lat
2537
2538CALL unproj(this%grid%proj, x, y, lon, lat)
2539
2540END SUBROUTINE griddim_coord_unproj
2541
2542
2543! Computes and sets the grid parameters required to compute
2544! coordinates of grid points in the projected system.
2545! probably meaningless
2546!SUBROUTINE griddim_proj(this)
2547!TYPE(griddim_def),INTENT(inout) :: this !< definition of projection and grid
2548!
2549!CALL proj(this, this%dim%lon(1,1), this%dim%lat(1,1), &
2550! this%grid%grid%xmin, this%grid%grid%ymin)
2551!
2552!CALL proj(this, this%dim%lon(this%dim%nx,this%dim%ny), &
2553! this%dim%lat(this%dim%nx,this%dim%ny), &
2554! this%grid%grid%xmax, this%grid%grid%ymax)
2555!
2556!END SUBROUTINE griddim_proj
2557
2558!> Computes the geographical coordinates of all the grid points in the
2559!! \a griddim_def object and stores them in the object itself. The \a
2560!! this::dim::lon and \a this::dim:lat arrays are allocated, and upon
2561!! return they will contain the coordinates' values. These arrays can
2562!! be directly accessed by the user, and their association status
2563!! should be checked with the \a ASSOCIATED() intrinsic before using
2564!! them.
2565SUBROUTINE griddim_unproj(this)
2566TYPE(griddim_def),INTENT(inout) :: this !< definition of projection and grid
2567
2568IF (.NOT.c_e(this%dim%nx) .OR. .NOT.c_e(this%dim%ny)) RETURN
2569CALL alloc(this%dim)
2570CALL griddim_unproj_internal(this)
2571
2572END SUBROUTINE griddim_unproj
2573
2574! internal subroutine needed for allocating automatic arrays
2575SUBROUTINE griddim_unproj_internal(this)
2576TYPE(griddim_def),INTENT(inout) ::this ! definition of projection and grid
2577
2578DOUBLE PRECISION :: x(this%dim%nx,this%dim%ny), y(this%dim%nx,this%dim%ny)
2579
2580CALL grid_rect_coordinates(this%grid%grid, x, y)
2581CALL unproj(this, x, y, this%dim%lon, this%dim%lat)
2582
2583END SUBROUTINE griddim_unproj_internal
2584
2585
2586!> Query the object content.
2587SUBROUTINE griddim_get_val(this, nx, ny, &
2588 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
2589 proj, proj_type, lov, zone, xoff, yoff, &
2590 longitude_south_pole, latitude_south_pole, angle_rotation, &
2591 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
2592 latin1, latin2, lad, projection_center_flag, &
2593 ellips_smaj_axis, ellips_flatt, ellips_type)
2594TYPE(griddim_def),INTENT(in) :: this !< object to be queried
2595INTEGER,INTENT(out),OPTIONAL :: nx !< number of points along the x axis
2596INTEGER,INTENT(out),OPTIONAL :: ny !< number of points along the y axis
2597!> longitudini e latitudini minime e massime
2598DOUBLE PRECISION,INTENT(out),OPTIONAL :: xmin, xmax, ymin, ymax !< grid extremes in projection units (degrees or meters depending on the projection type)
2599!> grid steps in x and y directions
2600DOUBLE PRECISION,INTENT(out),OPTIONAL :: dx, dy
2601!> Resolved u- and v- components of vector quantities relative to 0=the easterly and northerly directions
2602!! 1=the defined grid in the direction of increasing x and y (or i and j) coordinates respectively (0=north, 128=south)
2603INTEGER,INTENT(out),OPTIONAL :: component_flag
2604TYPE(geo_proj),INTENT(out),OPTIONAL :: proj !< the complete projection object associated
2605CHARACTER(len=*),INTENT(out),OPTIONAL :: proj_type !< type of projection
2606DOUBLE PRECISION,INTENT(out),OPTIONAL :: lov !< line of view, also known as reference longitude or orientation of the grid (polar projections)
2607INTEGER,INTENT(out),OPTIONAL :: zone !< Earth zone (mainly for UTM), sets lov to the correct zone central meridian
2608DOUBLE PRECISION,INTENT(out),OPTIONAL :: xoff !< offset on x axis (false easting)
2609DOUBLE PRECISION,INTENT(out),OPTIONAL :: yoff !< offset on y axis (false northing)
2610DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_south_pole !< longitude of the southern pole of projection
2611DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_south_pole !< latitude of the southern pole of projection
2612DOUBLE PRECISION,INTENT(out),OPTIONAL :: angle_rotation !< angle of rotation of projection
2613DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_stretch_pole !< longitude of the pole of stretching
2614DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_stretch_pole !< latitude of the pole of stretching
2615DOUBLE PRECISION,INTENT(out),OPTIONAL :: stretch_factor !< stretching factor
2616DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin1 !< first standard latitude from main pole (Lambert)
2617DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin2 !< second standard latitude from main pole (Lambert)
2618DOUBLE PRECISION,INTENT(out),OPTIONAL :: lad !< latitude at which dx and dy (in m) are specified (Lambert, grib2 only)
2619INTEGER,INTENT(out),OPTIONAL :: projection_center_flag !< flag indicating which pole is represented
2620DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_smaj_axis !< Earth semi-major axis
2621DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_flatt !< Earth flattening
2622INTEGER,INTENT(out),OPTIONAL :: ellips_type !< number in the interval [1,nellips] indicating a predefined ellipsoid, alternative to the previous arguments
2623
2624IF (PRESENT(nx)) nx = this%dim%nx
2625IF (PRESENT(ny)) ny = this%dim%ny
2626
2627IF (PRESENT(proj)) proj = this%grid%proj
2628
2629CALL get_val(this%grid%proj, proj_type=proj_type, lov=lov, zone=zone, &
2630 xoff=xoff, yoff=yoff, &
2631 longitude_south_pole=longitude_south_pole, &
2632 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
2633 longitude_stretch_pole=longitude_stretch_pole, &
2634 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
2635 latin1=latin1, latin2=latin2, lad=lad, &
2636 projection_center_flag=projection_center_flag, &
2637 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
2638 ellips_type=ellips_type)
2639
2640CALL get_val(this%grid%grid, &
2641 xmin, xmax, ymin, ymax, dx, dy, component_flag)
2642
2643END SUBROUTINE griddim_get_val
2644
2645
2646!> Set the object content.
2647SUBROUTINE griddim_set_val(this, nx, ny, &
2648 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
2649 proj_type, lov, zone, xoff, yoff, &
2650 longitude_south_pole, latitude_south_pole, angle_rotation, &
2651 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
2652 latin1, latin2, lad, projection_center_flag, &
2653 ellips_smaj_axis, ellips_flatt, ellips_type)
2654TYPE(griddim_def),INTENT(inout) :: this !< object to be queried
2655INTEGER,INTENT(in),OPTIONAL :: nx !< number of points along the x axis
2656INTEGER,INTENT(in),OPTIONAL :: ny !< number of points along the y axis
2657!> longitudini e latitudini minime e massime
2658DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin, xmax, ymin, ymax !< grid extremes in projection units (degrees or meters depending on the projection type)
2659!> grid steps in x and y directions
2660DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx, dy
2661!> Resolved u- and v- components of vector quantities relative to 0=the easterly and northerly directions
2662!! 1=the defined grid in the direction of increasing x and y (or i and j) coordinates respectively (0=north, 128=south)
2663INTEGER,INTENT(in),OPTIONAL :: component_flag
2664CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type !< type of projection
2665DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov !< line of view, also known as reference longitude or orientation of the grid (polar projections)
2666INTEGER,INTENT(in),OPTIONAL :: zone !< Earth zone (mainly for UTM), sets lov to the correct zone central meridian
2667DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff !< offset on x axis (false easting)
2668DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff !< offset on y axis (false northing)
2669DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole !< longitude of the southern pole of projection
2670DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole !< latitude of the southern pole of projection
2671DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation !< angle of rotation of projection
2672DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole !< longitude of the pole of stretching
2673DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole !< latitude of the pole of stretching
2674DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor !< stretching factor
2675DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1 !< first standard latitude from main pole (Lambert)
2676DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2 !< second standard latitude from main pole (Lambert)
2677DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad !< latitude at which dx and dy (in m) are specified (Lambert, grib2 only)
2678INTEGER,INTENT(in),OPTIONAL :: projection_center_flag !< flag indicating which pole is represented
2679DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis !< Earth semi-major axis
2680DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt !< Earth flattening
2681INTEGER,INTENT(in),OPTIONAL :: ellips_type !< number in the interval [1,nellips] indicating a predefined ellipsoid, alternative to the previous arguments
2682
2683IF (PRESENT(nx)) this%dim%nx = nx
2684IF (PRESENT(ny)) this%dim%ny = ny
2685
2686CALL set_val(this%grid%proj, proj_type=proj_type, lov=lov, zone=zone, &
2687 xoff=xoff, yoff=yoff, longitude_south_pole=longitude_south_pole, &
2688 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
2689 longitude_stretch_pole=longitude_stretch_pole, &
2690 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
2691 latin1=latin1, latin2=latin2, lad=lad, &
2692 projection_center_flag=projection_center_flag, &
2693 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
2694 ellips_type=ellips_type)
2695
2696CALL set_val(this%grid%grid, &
2697 xmin, xmax, ymin, ymax, dx, dy, component_flag)
2698
2699END SUBROUTINE griddim_set_val
2700
2701
2702!> This method reads from a Fortran file unit the contents of the
2703!! object \a this. The record to be read must have been written with
2704!! the ::write_unit method. The method works both on formatted and
2705!! unformatted files.
2706SUBROUTINE griddim_read_unit(this, unit)
2707TYPE(griddim_def),INTENT(out) :: this !< object to be read
2708INTEGER, INTENT(in) :: unit !< unit from which to read, it must be an opened Fortran file unit
2709
2710
2711CALL read_unit(this%dim, unit)
2712CALL read_unit(this%grid%proj, unit)
2713CALL read_unit(this%grid%grid, unit)
2714
2715END SUBROUTINE griddim_read_unit
2716
2717
2718!> This method writes on a Fortran file unit the contents of the
2719!! object \a this. The record can successively be read by the
2720!! ::read_unit method. The method works both on formatted and
2721!! unformatted files.
2722SUBROUTINE griddim_write_unit(this, unit)
2723TYPE(griddim_def),INTENT(in) :: this !< object to be written
2724INTEGER, INTENT(in) :: unit !< unit where to write, it must be an opened Fortran file unit
2725
2726
2727CALL write_unit(this%dim, unit)
2728CALL write_unit(this%grid%proj, unit)
2729CALL write_unit(this%grid%grid, unit)
2730
2731END SUBROUTINE griddim_write_unit
2732
2733
2734!> Euristically determine the approximate central longitude of the
2735!! grid in degrees.
2736!! The method depends on the projection used.
2737FUNCTION griddim_central_lon(this) RESULT(lon)
2738TYPE(griddim_def),INTENT(inout) :: this !< grid descriptor
2739
2740DOUBLE PRECISION :: lon
2741
2742CALL griddim_pistola_central_lon(this, lon)
2743
2744END FUNCTION griddim_central_lon
2745
2746
2747!> Euristically reset the approximate central longitude of the
2748!! grid to a value compatible to the provided longitude \a lonref.
2749!! The method depends on the projection used.
2750SUBROUTINE griddim_set_central_lon(this, lonref)
2751TYPE(griddim_def),INTENT(inout) :: this !< grid descriptor
2752DOUBLE PRECISION,INTENT(in) :: lonref !< reference longitude
2753
2754DOUBLE PRECISION :: lon
2755
2756CALL griddim_pistola_central_lon(this, lon, lonref)
2757
2758END SUBROUTINE griddim_set_central_lon
2759
2760
2761! internal subroutine for performing tasks common to the prevous two
2762SUBROUTINE griddim_pistola_central_lon(this, lon, lonref)
2763TYPE(griddim_def),INTENT(inout) :: this ! grid descriptor
2764DOUBLE PRECISION,INTENT(inout) :: lon ! central longitude
2765DOUBLE PRECISION,INTENT(in),OPTIONAL :: lonref ! reference longitude
2766
2767INTEGER :: unit
2768DOUBLE PRECISION :: lonsp, latsp, londelta, lov, lonrot
2769CHARACTER(len=80) :: ptype
2770
2771lon = dmiss
2772CALL get_val(this%grid%proj, unit=unit)
2773IF (unit == geo_proj_unit_meter) THEN ! it is a plane projection
2774 CALL get_val(this%grid%proj, lov=lon)
2775 IF (PRESENT(lonref)) THEN
2776 CALL long_reset_to_cart_closest(lov, lonref)
2777 CALL set_val(this%grid%proj, lov=lon)
2778 ENDIF
2779
2780ELSE IF (unit == geo_proj_unit_degree) THEN ! it is a spheric projection
2781 CALL get_val(this%grid%proj, proj_type=ptype, &
2782 longitude_south_pole=lonsp, latitude_south_pole=latsp)
2783 SELECT CASE(ptype)
2784 CASE('rotated_ll','stretched_rotated_ll') ! use origin of rotated system
2785 IF (latsp < 0.0d0) THEN
2786 lon = lonsp
2787 IF (PRESENT(lonref)) THEN
2788 CALL long_reset_to_cart_closest(lon, lonref)
2789 CALL set_val(this%grid%proj, longitude_south_pole=lonref)
2790! now reset rotated coordinates around zero
2791 IF (c_e(this%grid%grid%xmin) .AND. c_e(this%grid%grid%xmax)) THEN
2792 lonrot = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
2793 ENDIF
2794 londelta = lonrot
2795 CALL long_reset_to_cart_closest(londelta, 0.0d0)
2796 londelta = londelta - lonrot
2797 this%grid%grid%xmin = this%grid%grid%xmin + londelta
2798 this%grid%grid%xmax = this%grid%grid%xmax + londelta
2799 ENDIF
2800 ELSE ! this part to be checked
2801 lon = modulo(lonsp + 180.0d0, 360.0d0)
2802! IF (PRESENT(lonref)) THEN
2803! CALL long_reset_to_cart_closest(lov, lonref)
2804! CALL set_val(this%grid%proj, longitude_south_pole=lonref)
2805! ENDIF
2806 ENDIF
2807 CASE default ! use real grid limits
2808 IF (c_e(this%grid%grid%xmin) .AND. c_e(this%grid%grid%xmax)) THEN
2809 lon = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
2810 ENDIF
2811 IF (PRESENT(lonref)) THEN
2812 londelta = lon
2813 CALL long_reset_to_cart_closest(londelta, lonref)
2814 londelta = londelta - lon
2815 this%grid%grid%xmin = this%grid%grid%xmin + londelta
2816 this%grid%grid%xmax = this%grid%grid%xmax + londelta
2817 ENDIF
2818 END SELECT
2819ENDIF
2820
2821END SUBROUTINE griddim_pistola_central_lon
2822
2823
2824!> Generates coordinates of every point of a generic grid from the
2825!! grid description. The number of grid points along both direction is
2826!! guessed from the shape of x and y arrays, which must be conformal.
2827SUBROUTINE griddim_gen_coord(this, x, y)
2828TYPE(griddim_def),INTENT(in) :: this !< generic grid descriptor
2829DOUBLE PRECISION,INTENT(out) :: x(:,:) !< x coordinate of every point, linearly computed between grid extremes
2830DOUBLE PRECISION,INTENT(out) :: y(:,:) !< y coordinate of every point, linearly computed between grid extremes, it should have the same shape as x(:,:)
2831
2832
2833CALL grid_rect_coordinates(this%grid%grid, x, y)
2834
2835END SUBROUTINE griddim_gen_coord
2836
2837
2838!> Compute and return grid steps.
2839SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
2840TYPE(griddim_def), INTENT(in) :: this !< generic grid descriptor
2841INTEGER,INTENT(in) :: nx !< number of points along x direction
2842INTEGER,INTENT(in) :: ny !< number of points along y direction
2843DOUBLE PRECISION,INTENT(out) :: dx !< grid step along x direction
2844DOUBLE PRECISION,INTENT(out) :: dy !< grid step along y direction
2845
2846CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
2847
2848END SUBROUTINE griddim_steps
2849
2850
2851!> Compute and set grid steps.
2852SUBROUTINE griddim_setsteps(this)
2853TYPE(griddim_def), INTENT(inout) :: this !< generic grid descriptor
2854
2855CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
2856
2857END SUBROUTINE griddim_setsteps
2858
2859
2860! TODO
2861! bisogna sviluppare gli altri operatori
2862ELEMENTAL FUNCTION grid_eq(this, that) RESULT(res)
2863TYPE(grid_def),INTENT(IN) :: this, that
2864LOGICAL :: res
2865
2866res = this%proj == that%proj .AND. &
2867 this%grid == that%grid
2868
2869END FUNCTION grid_eq
2870
2871
2872ELEMENTAL FUNCTION griddim_eq(this, that) RESULT(res)
2873TYPE(griddim_def),INTENT(IN) :: this, that
2874LOGICAL :: res
2875
2876res = this%grid == that%grid .AND. &
2877 this%dim == that%dim
2878
2879END FUNCTION griddim_eq
2880
2881
2882ELEMENTAL FUNCTION grid_ne(this, that) RESULT(res)
2883TYPE(grid_def),INTENT(IN) :: this, that
2884LOGICAL :: res
2885
2886res = .NOT.(this == that)
2887
2888END FUNCTION grid_ne
2889
2890
2891ELEMENTAL FUNCTION griddim_ne(this, that) RESULT(res)
2892TYPE(griddim_def),INTENT(IN) :: this, that
2893LOGICAL :: res
2894
2895res = .NOT.(this == that)
2896
2897END FUNCTION griddim_ne
2898
2899
2900!> Import a griddim object from a grid_id object associated to a
2901!! supported gridded dataset driver (typically a grib message from
2902!! grib_api or a raster band from gdal). The griddim object is
2903!! populated with all the grid information (size, projection, etc.)
2904!! carried by the grid_id object provided.
2905SUBROUTINE griddim_import_grid_id(this, ingrid_id)
2906#ifdef HAVE_LIBGDAL
2907USE gdal
2908#endif
2909TYPE(griddim_def),INTENT(inout) :: this !< griddim object
2910TYPE(grid_id),INTENT(in) :: ingrid_id !< grid_id object with information about the grid
2911
2912#ifdef HAVE_LIBGRIBAPI
2913INTEGER :: gaid
2914#endif
2915#ifdef HAVE_LIBGDAL
2916TYPE(gdalrasterbandh) :: gdalid
2917#endif
2918CALL init(this)
2919
2920#ifdef HAVE_LIBGRIBAPI
2921gaid = grid_id_get_gaid(ingrid_id)
2922IF (c_e(gaid)) CALL griddim_import_gribapi(this, gaid)
2923#endif
2924#ifdef HAVE_LIBGDAL
2925gdalid = grid_id_get_gdalid(ingrid_id)
2926IF (gdalassociated(gdalid)) CALL griddim_import_gdal(this, gdalid, &
2927 grid_id_get_gdal_options(ingrid_id))
2928#endif
2929
2930END SUBROUTINE griddim_import_grid_id
2931
2932
2933!> Export a griddim object to a grid_id object associated to a
2934!! supported gridded dataset driver (typically a grib message from
2935!! grib_api). All the grid information (size, projection, etc.)
2936!! contained in the griddim object is exported to the grid_id object.
2937SUBROUTINE griddim_export_grid_id(this, outgrid_id)
2938#ifdef HAVE_LIBGDAL
2939USE gdal
2940#endif
2941TYPE(griddim_def),INTENT(in) :: this !< griddim object
2942TYPE(grid_id),INTENT(inout) :: outgrid_id !< grid_id object which will contain information about the grid
2943
2944#ifdef HAVE_LIBGRIBAPI
2945INTEGER :: gaid
2946#endif
2947#ifdef HAVE_LIBGDAL
2948TYPE(gdalrasterbandh) :: gdalid
2949#endif
2950
2951#ifdef HAVE_LIBGRIBAPI
2952gaid = grid_id_get_gaid(outgrid_id)
2953IF (c_e(gaid)) CALL griddim_export_gribapi(this, gaid)
2954#endif
2955#ifdef HAVE_LIBGDAL
2956gdalid = grid_id_get_gdalid(outgrid_id)
2957!IF (gdalassociated(gdalid)
2958! export for gdal not implemented, log?
2959#endif
2960
2961END SUBROUTINE griddim_export_grid_id
2962
2963
2964#ifdef HAVE_LIBGRIBAPI
2965! grib_api driver
2966SUBROUTINE griddim_import_gribapi(this, gaid)
2967USE grib_api
2968TYPE(griddim_def),INTENT(inout) :: this ! griddim object
2969INTEGER, INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
2970
2971DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, x1, y1, orient
2972INTEGER :: EditionNumber, iScansNegatively, jScansPositively, zone, datum, &
2973 reflon, ierr
2974
2975! Generic keys
2976CALL grib_get(gaid, 'typeOfGrid', this%grid%proj%proj_type)
2977#ifdef DEBUG
2978call l4f_category_log(this%category,l4f_debug, &
2979 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type))
2980#endif
2981CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
2982
2983IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
2984 CALL grib_get(gaid,'numberOfDataPoints', this%dim%nx)
2985 this%dim%ny = 1
2986 this%grid%grid%component_flag = 0
2987 CALL griddim_import_ellipsoid(this, gaid)
2988 RETURN
2989ENDIF
2990
2991! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
2992CALL grib_get(gaid, 'Ni', this%dim%nx)
2993CALL grib_get(gaid, 'Nj', this%dim%ny)
2994CALL griddim_import_ellipsoid(this, gaid) ! placed here, not valid for utm datum /= 1
2995
2996CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
2997CALL grib_get(gaid,'jScansPositively',jscanspositively)
2998CALL grib_get(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
2999
3000! Keys for rotated grids (checked through missing values)
3001CALL grib_get_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
3002 this%grid%proj%rotated%longitude_south_pole)
3003CALL grib_get_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
3004 this%grid%proj%rotated%latitude_south_pole)
3005CALL grib_get_dmiss(gaid,'angleOfRotationInDegrees', &
3006 this%grid%proj%rotated%angle_rotation)
3007
3008! Keys for stretched grids (checked through missing values)
3009! units must be verified, still experimental in grib_api
3010! # TODO: Is it a float? Is it signed?
3011IF (editionnumber == 1) THEN
3012 CALL grib_get_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
3013 this%grid%proj%stretched%longitude_stretch_pole)
3014 CALL grib_get_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
3015 this%grid%proj%stretched%latitude_stretch_pole)
3016 CALL grib_get_dmiss(gaid,'stretchingFactor', &
3017 this%grid%proj%stretched%stretch_factor)
3018ELSE IF (editionnumber == 2) THEN
3019 CALL grib_get_dmiss(gaid,'longitudeOfThePoleOfStretching', &
3020 this%grid%proj%stretched%longitude_stretch_pole)
3021 CALL grib_get_dmiss(gaid,'latitudeOfThePoleOfStretching', &
3022 this%grid%proj%stretched%latitude_stretch_pole)
3023 CALL grib_get_dmiss(gaid,'stretchingFactor', &
3024 this%grid%proj%stretched%stretch_factor)
3025 IF (c_e(this%grid%proj%stretched%stretch_factor)) &
3026 this%grid%proj%stretched%stretch_factor = this%grid%proj%stretched%stretch_factor*1.0d-6
3027ENDIF
3028
3029! Projection-dependent keys
3030SELECT CASE (this%grid%proj%proj_type)
3031
3032! Keys for spherical coordinate systems
3033CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
3034
3035 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3036 CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3037 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3038 CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3039
3040! longitudes are sometimes wrongly coded even in grib2 and even by the
3041! Metoffice!
3042! longitudeOfFirstGridPointInDegrees = 354.911;
3043! longitudeOfLastGridPointInDegrees = 363.311;
3044 CALL long_reset_0_360(lofirst)
3045 CALL long_reset_0_360(lolast)
3046
3047 IF (iscansnegatively == 0) THEN
3048 this%grid%grid%xmin = lofirst
3049 this%grid%grid%xmax = lolast
3050 ELSE
3051 this%grid%grid%xmax = lofirst
3052 this%grid%grid%xmin = lolast
3053 ENDIF
3054 IF (jscanspositively == 0) THEN
3055 this%grid%grid%ymax = lafirst
3056 this%grid%grid%ymin = lalast
3057 ELSE
3058 this%grid%grid%ymin = lafirst
3059 this%grid%grid%ymax = lalast
3060 ENDIF
3061
3062! reset longitudes in order to have a Cartesian plane
3063 IF (this%grid%grid%xmax < this%grid%grid%xmin) &
3064 this%grid%grid%xmin = this%grid%grid%xmin - 360.d0
3065
3066! compute dx and dy (should we get them from grib?)
3067 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3068
3069! Keys for polar projections
3070CASE ('polar_stereographic', 'lambert', 'albers')
3071
3072 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
3073 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
3074! latin1/latin2 may be missing (e.g. stereographic)
3075 CALL grib_get_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
3076 CALL grib_get_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
3077#ifdef DEBUG
3078 CALL l4f_category_log(this%category,l4f_debug, &
3079 "griddim_import_gribapi, latin1/2 "// &
3080 trim(to_char(this%grid%proj%polar%latin1))//" "// &
3081 trim(to_char(this%grid%proj%polar%latin2)))
3082#endif
3083! projection center flag, aka hemisphere
3084 CALL grib_get(gaid,'projectionCenterFlag',&
3085 this%grid%proj%polar%projection_center_flag, ierr)
3086 IF (ierr /= grib_success) THEN ! try center/centre
3087 CALL grib_get(gaid,'projectionCentreFlag',&
3088 this%grid%proj%polar%projection_center_flag)
3089 ENDIF
3090
3091 IF (iand(this%grid%proj%polar%projection_center_flag,64) == 1) THEN
3092 CALL l4f_category_log(this%category,l4f_error, &
3093 "griddim_import_gribapi, bi-polar projections not supported")
3094 CALL raise_error()
3095 ENDIF
3096! line of view, aka central meridian
3097 CALL grib_get(gaid,'LoVInDegrees',this%grid%proj%lov)
3098#ifdef DEBUG
3099 CALL l4f_category_log(this%category,l4f_debug, &
3100 "griddim_import_gribapi, central meridian "//trim(to_char(this%grid%proj%lov)))
3101#endif
3102
3103! latitude at which dx and dy are valid
3104 IF (editionnumber == 1) THEN
3105! ECMWF (gribex/grib_api) says: Grid lengths are in metres, at the
3106! 60-degree parallel nearest to the pole on the projection plane.
3107! IF (IAND(this%projection_center_flag, 128) == 0) THEN
3108! this%grid%proj%polar%lad = 60.D0
3109! ELSE
3110! this%grid%proj%polar%lad = -60.D0
3111! ENDIF
3112! WMO says: Grid lengths are in units of metres, at the secant cone
3113! intersection parallel nearest to the pole on the projection plane.
3114 this%grid%proj%polar%lad = this%grid%proj%polar%latin1
3115 ELSE IF (editionnumber == 2) THEN
3116 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3117 ENDIF
3118#ifdef DEBUG
3119 CALL l4f_category_log(this%category,l4f_debug, &
3120 "griddim_import_gribapi, lad "//trim(to_char(this%grid%proj%polar%lad)))
3121#endif
3122
3123! compute projected extremes from lon and lat of first point
3124 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3125 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3126 CALL long_reset_0_360(lofirst)
3127 CALL long_reset_to_cart_closest(this%grid%proj%lov, lofirst)
3128#ifdef DEBUG
3129 CALL l4f_category_log(this%category,l4f_debug, &
3130 "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
3131 CALL l4f_category_log(this%category,l4f_debug, &
3132 "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
3133#endif
3134
3135 CALL proj(this, lofirst, lafirst, x1, y1)
3136 IF (iscansnegatively == 0) THEN
3137 this%grid%grid%xmin = x1
3138 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
3139 ELSE
3140 this%grid%grid%xmax = x1
3141 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
3142 ENDIF
3143 IF (jscanspositively == 0) THEN
3144 this%grid%grid%ymax = y1
3145 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
3146 ELSE
3147 this%grid%grid%ymin = y1
3148 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
3149 ENDIF
3150! keep these values for personal pleasure
3151 this%grid%proj%polar%lon1 = lofirst
3152 this%grid%proj%polar%lat1 = lafirst
3153
3154! Keys for equatorial projections
3155CASE ('mercator')
3156
3157 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
3158 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
3159 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3160 this%grid%proj%lov = 0.0d0 ! not in grib
3161
3162 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3163 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3164
3165 CALL proj(this, lofirst, lafirst, x1, y1)
3166 IF (iscansnegatively == 0) THEN
3167 this%grid%grid%xmin = x1
3168 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
3169 ELSE
3170 this%grid%grid%xmax = x1
3171 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
3172 ENDIF
3173 IF (jscanspositively == 0) THEN
3174 this%grid%grid%ymax = y1
3175 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
3176 ELSE
3177 this%grid%grid%ymin = y1
3178 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
3179 ENDIF
3180
3181 IF (editionnumber == 2) THEN
3182 CALL grib_get(gaid,'orientationOfTheGridInDegrees',orient)
3183 IF (orient /= 0.0d0) THEN
3184 CALL l4f_category_log(this%category,l4f_error, &
3185 "griddim_import_gribapi, Mercator grid orientation != 0 not supported")
3186 CALL raise_error()
3187 ENDIF
3188 ENDIF
3189
3190#ifdef DEBUG
3191 CALL unproj(this, x1, y1, lofirst, lafirst)
3192 CALL l4f_category_log(this%category,l4f_debug, &
3193 "griddim_import_gribapi, unprojected first point "//t2c(lofirst)//" "// &
3194 t2c(lafirst))
3195
3196 CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3197 CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3198 CALL proj(this, lolast, lalast, x1, y1)
3199 CALL l4f_category_log(this%category,l4f_debug, &
3200 "griddim_import_gribapi, extremes from grib "//t2c(x1)//" "//t2c(y1))
3201 CALL l4f_category_log(this%category,l4f_debug, &
3202 "griddim_import_gribapi, extremes from proj "//t2c(this%grid%grid%xmin)// &
3203 " "//t2c(this%grid%grid%ymin)//" "//t2c(this%grid%grid%xmax)//" "// &
3204 t2c(this%grid%grid%ymax))
3205#endif
3206
3207CASE ('UTM')
3208
3209 CALL grib_get(gaid,'zone',zone)
3210
3211 CALL grib_get(gaid,'datum',datum)
3212 IF (datum == 0) THEN
3213 CALL grib_get(gaid,'referenceLongitude',reflon)
3214 CALL grib_get(gaid,'falseEasting',this%grid%proj%xoff)
3215 CALL grib_get(gaid,'falseNorthing',this%grid%proj%yoff)
3216 CALL set_val(this%grid%proj, zone=zone, lov=reflon/1.0d6)
3217 ELSE
3218 CALL l4f_category_log(this%category,l4f_error,'only datum 0 supported')
3219 CALL raise_fatal_error()
3220 ENDIF
3221
3222 CALL grib_get(gaid,'eastingOfFirstGridPoint',lofirst)
3223 CALL grib_get(gaid,'eastingOfLastGridPoint',lolast)
3224 CALL grib_get(gaid,'northingOfFirstGridPoint',lafirst)
3225 CALL grib_get(gaid,'northingOfLastGridPoint',lalast)
3226
3227 IF (iscansnegatively == 0) THEN
3228 this%grid%grid%xmin = lofirst
3229 this%grid%grid%xmax = lolast
3230 ELSE
3231 this%grid%grid%xmax = lofirst
3232 this%grid%grid%xmin = lolast
3233 ENDIF
3234 IF (jscanspositively == 0) THEN
3235 this%grid%grid%ymax = lafirst
3236 this%grid%grid%ymin = lalast
3237 ELSE
3238 this%grid%grid%ymin = lafirst
3239 this%grid%grid%ymax = lalast
3240 ENDIF
3241
3242! compute dx and dy (should we get them from grib?)
3243 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3244
3245CASE default
3246 CALL l4f_category_log(this%category,l4f_error, &
3247 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
3248 CALL raise_error()
3249
3250END SELECT
3251
3252CONTAINS
3253! utilities routines for grib_api, is there a better place?
3254SUBROUTINE grib_get_dmiss(gaid, key, value)
3255INTEGER,INTENT(in) :: gaid
3256CHARACTER(len=*),INTENT(in) :: key
3257DOUBLE PRECISION,INTENT(out) :: value
3258
3259INTEGER :: ierr
3260
3261CALL grib_get(gaid, key, value, ierr)
3262IF (ierr /= grib_success) value = dmiss
3263
3264END SUBROUTINE grib_get_dmiss
3265
3266SUBROUTINE grib_get_imiss(gaid, key, value)
3267INTEGER,INTENT(in) :: gaid
3268CHARACTER(len=*),INTENT(in) :: key
3269INTEGER,INTENT(out) :: value
3270
3271INTEGER :: ierr
3272
3273CALL grib_get(gaid, key, value, ierr)
3274IF (ierr /= grib_success) value = imiss
3275
3276END SUBROUTINE grib_get_imiss
3277
3278
3279SUBROUTINE griddim_import_ellipsoid(this, gaid)
3280TYPE(griddim_def),INTENT(inout) :: this
3281INTEGER,INTENT(in) :: gaid
3282
3283INTEGER :: shapeofearth, iv, is
3284DOUBLE PRECISION :: r1, r2
3285
3286IF (editionnumber == 2) THEN
3287 CALL grib_get(gaid, 'shapeOfTheEarth', shapeofearth)
3288 SELECT CASE(shapeofearth)
3289 CASE(0) ! spherical
3290 CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
3291 CASE(1) ! spherical generic
3292 CALL grib_get(gaid, 'scaleFactorOfRadiusOfSphericalEarth', is)
3293 CALL grib_get(gaid, 'scaledValueOfRadiusOfSphericalEarth', iv)
3294 r1 = dble(iv) / 10**is
3295 CALL set_val(this, ellips_smaj_axis=r1, ellips_flatt=0.0d0)
3296 CASE(2) ! iau65
3297 CALL set_val(this, ellips_smaj_axis=6378160.0d0, ellips_flatt=1.0d0/297.0d0)
3298 CASE(3,7) ! ellipsoidal generic
3299 CALL grib_get(gaid, 'scaleFactorOfEarthMajorAxis', is)
3300 CALL grib_get(gaid, 'scaledValueOfEarthMajorAxis', iv)
3301 r1 = dble(iv) / 10**is
3302 CALL grib_get(gaid, 'scaleFactorOfEarthMinorAxis', is)
3303 CALL grib_get(gaid, 'scaledValueOfEarthMinorAxis', iv)
3304 r2 = dble(iv) / 10**is
3305 IF (shapeofearth == 3) THEN ! km->m
3306 r1 = r1*1000.0d0
3307 r2 = r2*1000.0d0
3308 ENDIF
3309 IF (abs(r1) < 1.0d-6) THEN ! suspicious data read from grib
3310 CALL l4f_category_log(this%category,l4f_warn,'zero Earth major axis '// &
3311 'read from grib, going on with spherical Earth but the results may be wrong')
3312 CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
3313 ELSE
3314 CALL set_val(this, ellips_smaj_axis=r1, ellips_flatt=(r1-r2)/r1)
3315 ENDIF
3316 CASE(4) ! iag-grs80
3317 CALL set_val(this, ellips_type=ellips_grs80)
3318 CASE(5) ! wgs84
3319 CALL set_val(this, ellips_type=ellips_wgs84)
3320 CASE(6) ! spherical
3321 CALL set_val(this, ellips_smaj_axis=6371229.0d0, ellips_flatt=0.0d0)
3322! CASE(7) ! google earth-like?
3323 CASE default
3324 CALL l4f_category_log(this%category,l4f_error,'shapeOfTheEarth '// &
3325 t2c(shapeofearth)//' not supported in grib2')
3326 CALL raise_fatal_error()
3327
3328 END SELECT
3329
3330ELSE
3331
3332 CALL grib_get(gaid, 'earthIsOblate', shapeofearth)
3333 IF (shapeofearth == 0) THEN ! spherical
3334 CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
3335 ELSE ! iau65
3336 CALL set_val(this, ellips_smaj_axis=6378160.0d0, ellips_flatt=1.0d0/297.0d0)
3337 ENDIF
3338
3339ENDIF
3340
3341END SUBROUTINE griddim_import_ellipsoid
3342
3343
3344END SUBROUTINE griddim_import_gribapi
3345
3346
3347! grib_api driver
3348SUBROUTINE griddim_export_gribapi(this, gaid)
3349USE grib_api
3350TYPE(griddim_def),INTENT(in) :: this ! griddim object
3351INTEGER, INTENT(inout) :: gaid ! grib_api id of the grib loaded in memory to export
3352
3353INTEGER :: EditionNumber, iScansNegatively, jScansPositively, nv, pvl, zone, ierr
3354DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, reflon
3355DOUBLE PRECISION :: sdx, sdy, ratio, tol, ex, ey
3356
3357
3358! Generic keys
3359CALL grib_get(gaid,'GRIBEditionNumber', editionnumber)
3360! the following required since eccodes
3361IF (editionnumber == 2) CALL grib_set(gaid,'shapeOfTheEarth', 6)
3362CALL grib_set(gaid,'typeOfGrid', this%grid%proj%proj_type)
3363#ifdef DEBUG
3364CALL l4f_category_log(this%category,l4f_debug, &
3365 "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
3366#endif
3367
3368IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
3369 CALL grib_set(gaid,'numberOfDataPoints', this%dim%nx)
3370! reenable when possible
3371! CALL griddim_export_ellipsoid(this, gaid)
3372 RETURN
3373ENDIF
3374
3375
3376! Edition dependent setup
3377IF (editionnumber == 1) THEN
3378 ratio = 1.d3
3379ELSE IF (editionnumber == 2) THEN
3380 ratio = 1.d6
3381ELSE
3382 ratio = 0.0d0 ! signal error?!
3383ENDIF
3384
3385! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
3386CALL griddim_export_ellipsoid(this, gaid)
3387CALL grib_set(gaid, 'Ni', this%dim%nx)
3388CALL grib_set(gaid, 'Nj', this%dim%ny)
3389CALL grib_set(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
3390
3391CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
3392CALL grib_get(gaid,'jScansPositively',jscanspositively)
3393
3394! Keys for rotated grids (checked through missing values and/or error code)
3395!SELECT CASE (this%grid%proj%proj_type)
3396!CASE ('rotated_ll', 'stretched_rotated_ll', 'polar_stereographic', 'lambert', 'albers')
3397CALL grib_set_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
3398 this%grid%proj%rotated%longitude_south_pole, 0.0d0)
3399CALL grib_set_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
3400 this%grid%proj%rotated%latitude_south_pole, -90.0d0)
3401IF (editionnumber == 1) THEN
3402 CALL grib_set_dmiss(gaid,'angleOfRotationInDegrees', &
3403 this%grid%proj%rotated%angle_rotation, 0.0d0)
3404ELSE IF (editionnumber == 2)THEN
3405 CALL grib_set_dmiss(gaid,'angleOfRotationOfProjectionInDegrees', &
3406 this%grid%proj%rotated%angle_rotation, 0.0d0)
3407ENDIF
3408
3409! Keys for stretched grids (checked through missing values and/or error code)
3410! units must be verified, still experimental in grib_api
3411! # TODO: Is it a float? Is it signed?
3412IF (editionnumber == 1) THEN
3413 CALL grib_set_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
3414 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
3415 CALL grib_set_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
3416 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
3417 CALL grib_set_dmiss(gaid,'stretchingFactor', &
3418 this%grid%proj%stretched%stretch_factor, 1.0d0)
3419ELSE IF (editionnumber == 2) THEN
3420 CALL grib_set_dmiss(gaid,'longitudeOfThePoleOfStretching', &
3421 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
3422 CALL grib_set_dmiss(gaid,'latitudeOfThePoleOfStretching', &
3423 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
3424 CALL grib_set_dmiss(gaid,'stretchingFactor', &
3425 this%grid%proj%stretched%stretch_factor, 1.0d6, 1.0d6)
3426ENDIF
3427
3428! Projection-dependent keys
3429SELECT CASE (this%grid%proj%proj_type)
3430
3431! Keys for sphaerical coordinate systems
3432CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
3433
3434 IF (iscansnegatively == 0) THEN
3435 lofirst = this%grid%grid%xmin
3436 lolast = this%grid%grid%xmax
3437 ELSE
3438 lofirst = this%grid%grid%xmax
3439 lolast = this%grid%grid%xmin
3440 ENDIF
3441 IF (jscanspositively == 0) THEN
3442 lafirst = this%grid%grid%ymax
3443 lalast = this%grid%grid%ymin
3444 ELSE
3445 lafirst = this%grid%grid%ymin
3446 lalast = this%grid%grid%ymax
3447 ENDIF
3448
3449! reset lon in standard grib 2 definition [0,360]
3450 IF (editionnumber == 1) THEN
3451 CALL long_reset_m180_360(lofirst)
3452 CALL long_reset_m180_360(lolast)
3453 ELSE IF (editionnumber == 2) THEN
3454 CALL long_reset_0_360(lofirst)
3455 CALL long_reset_0_360(lolast)
3456 ENDIF
3457
3458 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3459 CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3460 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3461 CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3462
3463! test relative coordinate truncation error with respect to tol
3464! tol should be tuned
3465 sdx = this%grid%grid%dx*ratio
3466 sdy = this%grid%grid%dy*ratio
3467! error is computed relatively to the whole grid extension
3468 ex = abs(((sdx - nint(sdx))/ratio*this%dim%nx)/ &
3469 (this%grid%grid%xmax-this%grid%grid%xmin))
3470 ey = abs(((sdy - nint(sdy))/ratio*this%dim%ny)/ &
3471 (this%grid%grid%ymax-this%grid%grid%ymin))
3472 tol = 1.0d-3
3473 IF (ex > tol .OR. ey > tol) THEN
3474#ifdef DEBUG
3475 CALL l4f_category_log(this%category,l4f_debug, &
3476 "griddim_export_gribapi, error on x "//&
3477 trim(to_char(ex))//"/"//t2c(tol))
3478 CALL l4f_category_log(this%category,l4f_debug, &
3479 "griddim_export_gribapi, error on y "//&
3480 trim(to_char(ey))//"/"//t2c(tol))
3481#endif
3482! previous frmula relative to a single grid step
3483! tol = 1.0d0/ratio
3484! IF (ABS(NINT(sdx)/sdx - 1.0d0) > tol .OR. ABS(NINT(sdy)/sdy - 1.0d0) > tol) THEN
3485!#ifdef DEBUG
3486! CALL l4f_category_log(this%category,L4F_DEBUG, &
3487! "griddim_export_gribapi, dlon relative error: "//&
3488! TRIM(to_char(ABS(NINT(sdx)/sdx - 1.0d0)))//">"//TRIM(to_char(tol)))
3489! CALL l4f_category_log(this%category,L4F_DEBUG, &
3490! "griddim_export_gribapi, dlat relative error: "//&
3491! TRIM(to_char(ABS(NINT(sdy)/sdy - 1.0d0)))//">"//TRIM(to_char(tol)))
3492!#endif
3493 CALL l4f_category_log(this%category,l4f_info, &
3494 "griddim_export_gribapi, increments not given: inaccurate!")
3495 CALL grib_set_missing(gaid,'Di')
3496 CALL grib_set_missing(gaid,'Dj')
3497 CALL grib_set(gaid,'ijDirectionIncrementGiven',0)
3498 ELSE
3499#ifdef DEBUG
3500 CALL l4f_category_log(this%category,l4f_debug, &
3501 "griddim_export_gribapi, setting increments: "// &
3502 trim(to_char(this%grid%grid%dx))//' '//trim(to_char(this%grid%grid%dy)))
3503#endif
3504 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
3505 CALL grib_set(gaid,'iDirectionIncrement',nint(this%grid%grid%dx*ratio))
3506 CALL grib_set(gaid,'jDirectionIncrement',nint(this%grid%grid%dy*ratio))
3507! this does not work in grib_set
3508! CALL grib_set(gaid,'iDirectionIncrementInDegrees',this%grid%grid%dx)
3509! CALL grib_set(gaid,'jDirectionIncrementInDegrees',this%grid%grid%dy)
3510 ENDIF
3511
3512! Keys for polar projections
3513CASE ('polar_stereographic', 'lambert', 'albers')
3514! increments are required
3515 CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
3516 CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
3517 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
3518! latin1/latin2 may be missing (e.g. stereographic)
3519 CALL grib_set_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
3520 CALL grib_set_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
3521! projection center flag, aka hemisphere
3522 CALL grib_set(gaid,'projectionCenterFlag',&
3523 this%grid%proj%polar%projection_center_flag, ierr)
3524 IF (ierr /= grib_success) THEN ! try center/centre
3525 CALL grib_set(gaid,'projectionCentreFlag',&
3526 this%grid%proj%polar%projection_center_flag)
3527 ENDIF
3528
3529
3530! line of view, aka central meridian
3531 CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
3532! latitude at which dx and dy are valid
3533 IF (editionnumber == 2) THEN
3534 CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3535 ENDIF
3536
3537! compute lon and lat of first point from projected extremes
3538 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
3539 lofirst, lafirst)
3540! reset lon in standard grib 2 definition [0,360]
3541 IF (editionnumber == 1) THEN
3542 CALL long_reset_m180_360(lofirst)
3543 ELSE IF (editionnumber == 2) THEN
3544 CALL long_reset_0_360(lofirst)
3545 ENDIF
3546 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3547 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3548
3549! Keys for equatorial projections
3550CASE ('mercator')
3551
3552! increments are required
3553 CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
3554 CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
3555 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
3556
3557! line of view, aka central meridian, not in grib
3558! CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
3559! latitude at which dx and dy are valid
3560 CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3561
3562! compute lon and lat of first and last points from projected extremes
3563 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
3564 lofirst, lafirst)
3565 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3566 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3567 CALL get_unproj_last(this, iscansnegatively, jscanspositively, &
3568 lolast, lalast)
3569 CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3570 CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3571
3572 IF (editionnumber == 2) THEN
3573 CALL grib_set(gaid,'orientationOfTheGridInDegrees',0.)
3574 ENDIF
3575
3576CASE ('UTM')
3577
3578 CALL grib_set(gaid,'datum',0)
3579 CALL get_val(this, zone=zone, lov=reflon)
3580 CALL grib_set(gaid,'referenceLongitude',nint(reflon/1.0d6))
3581 CALL grib_set(gaid,'falseEasting',this%grid%proj%xoff)
3582 CALL grib_set(gaid,'falseNorthing',this%grid%proj%yoff)
3583
3584 CALL grib_set(gaid,'iDirectionIncrement',this%grid%grid%dx)
3585 CALL grib_set(gaid,'jDirectionIncrement',this%grid%grid%dy)
3586
3587!res/scann ??
3588
3589 CALL grib_set(gaid,'zone',zone)
3590
3591 IF (iscansnegatively == 0) THEN
3592 lofirst = this%grid%grid%xmin
3593 lolast = this%grid%grid%xmax
3594 ELSE
3595 lofirst = this%grid%grid%xmax
3596 lolast = this%grid%grid%xmin
3597 ENDIF
3598 IF (jscanspositively == 0) THEN
3599 lafirst = this%grid%grid%ymax
3600 lalast = this%grid%grid%ymin
3601 ELSE
3602 lafirst = this%grid%grid%ymin
3603 lalast = this%grid%grid%ymax
3604 ENDIF
3605
3606 CALL grib_set(gaid,'eastingOfFirstGridPoint',lofirst)
3607 CALL grib_set(gaid,'eastingOfLastGridPoint',lolast)
3608 CALL grib_set(gaid,'northingOfFirstGridPoint',lafirst)
3609 CALL grib_set(gaid,'northingOfLastGridPoint',lalast)
3610
3611CASE default
3612 CALL l4f_category_log(this%category,l4f_error, &
3613 "griddim_export_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
3614 CALL raise_error()
3615
3616END SELECT
3617
3618! hack for position of vertical coordinate parameters
3619! buggy in grib_api
3620IF (editionnumber == 1) THEN
3621! CALL grib_get(gaid,"PVPresent",pvp) ! alias, probably useless
3622 CALL grib_get(gaid,"NV",nv)
3623#ifdef DEBUG
3624 CALL l4f_category_log(this%category,l4f_debug,"griddim_export_gribapi, coding "// &
3625 trim(to_char(nv))//" vertical coordinate parameters")
3626#endif
3627
3628 IF (nv == 0) THEN
3629 pvl = 255
3630 ELSE
3631 SELECT CASE (this%grid%proj%proj_type)
3632 CASE ('regular_ll') ! check whether "29-32 Set to zero (reserved)" are required
3633 pvl = 33
3634 CASE ('polar_stereographic')
3635 pvl = 33
3636 CASE ('rotated_ll', 'stretched_ll', 'lambert', 'albers')
3637 pvl = 43
3638 CASE ('stretched_rotated_ll')
3639 pvl = 43
3640 CASE DEFAULT
3641 pvl = 43 !?
3642 END SELECT
3643 ENDIF
3644
3645 CALL grib_set(gaid,"pvlLocation",pvl)
3646#ifdef DEBUG
3647 CALL l4f_category_log(this%category,l4f_debug,"griddim_export_gribapi, coding "// &
3648 trim(to_char(pvl))//" as vertical coordinate parameter location")
3649#endif
3650
3651ENDIF
3652
3653
3654CONTAINS
3655! utilities routines for grib_api, is there a better place?
3656SUBROUTINE grib_set_dmiss(gaid, key, val, default, factor)
3657INTEGER,INTENT(in) :: gaid
3658CHARACTER(len=*),INTENT(in) :: key
3659DOUBLE PRECISION,INTENT(in) :: val
3660DOUBLE PRECISION,INTENT(in),OPTIONAL :: default
3661DOUBLE PRECISION,INTENT(in),OPTIONAL :: factor
3662
3663INTEGER :: ierr
3664
3665IF (c_e(val)) THEN
3666 IF (PRESENT(factor)) THEN
3667 CALL grib_set(gaid, key, val*factor, ierr)
3668 ELSE
3669 CALL grib_set(gaid, key, val, ierr)
3670 ENDIF
3671ELSE IF (PRESENT(default)) THEN
3672 CALL grib_set(gaid, key, default, ierr)
3673ENDIF
3674
3675END SUBROUTINE grib_set_dmiss
3676
3677SUBROUTINE grib_set_imiss(gaid, key, value, default)
3678INTEGER,INTENT(in) :: gaid
3679CHARACTER(len=*),INTENT(in) :: key
3680INTEGER,INTENT(in) :: value
3681INTEGER,INTENT(in),OPTIONAL :: default
3682
3683INTEGER :: ierr
3684
3685IF (c_e(value)) THEN
3686 CALL grib_set(gaid, key, value, ierr)
3687ELSE IF (PRESENT(default)) THEN
3688 CALL grib_set(gaid, key, default, ierr)
3689ENDIF
3690
3691END SUBROUTINE grib_set_imiss
3692
3693SUBROUTINE griddim_export_ellipsoid(this, gaid)
3694TYPE(griddim_def),INTENT(in) :: this
3695INTEGER,INTENT(in) :: gaid
3696
3697INTEGER :: ellips_type, ierr
3698DOUBLE PRECISION :: r1, r2, f
3699
3700CALL get_val(this, ellips_smaj_axis=r1, ellips_flatt=f, ellips_type=ellips_type)
3701
3702IF (editionnumber == 2) THEN
3703
3704! why does it produce a message even with ierr?
3705 CALL grib_set_missing(gaid, 'scaleFactorOfRadiusOfSphericalEarth', ierr)
3706 CALL grib_set_missing(gaid, 'scaledValueOfRadiusOfSphericalEarth', ierr)
3707 CALL grib_set_missing(gaid, 'scaleFactorOfEarthMajorAxis', ierr)
3708 CALL grib_set_missing(gaid, 'scaledValueOfEarthMajorAxis', ierr)
3709 CALL grib_set_missing(gaid, 'scaleFactorOfEarthMinorAxis', ierr)
3710 CALL grib_set_missing(gaid, 'scaledValueOfEarthMinorAxis', ierr)
3711
3712 SELECT CASE(ellips_type)
3713 CASE(ellips_grs80) ! iag-grs80
3714 CALL grib_set(gaid, 'shapeOfTheEarth', 4)
3715 CASE(ellips_wgs84) ! wgs84
3716 CALL grib_set(gaid, 'shapeOfTheEarth', 5)
3717 CASE default
3718 IF (f == 0.0d0) THEN ! spherical Earth
3719 IF (r1 == 6367470.0d0) THEN ! spherical
3720 CALL grib_set(gaid, 'shapeOfTheEarth', 0)
3721 ELSE IF (r1 == 6371229.0d0) THEN ! spherical
3722 CALL grib_set(gaid, 'shapeOfTheEarth', 6)
3723 ELSE ! spherical generic
3724 CALL grib_set(gaid, 'shapeOfTheEarth', 1)
3725 CALL grib_set(gaid, 'scaleFactorOfRadiusOfSphericalEarth', 2)
3726 CALL grib_set(gaid, 'scaledValueOfRadiusOfSphericalEarth', int(r1*100.0d0))
3727 ENDIF
3728 ELSE ! ellipsoidal
3729 IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
3730 CALL grib_set(gaid, 'shapeOfTheEarth', 2)
3731 ELSE ! ellipsoidal generic
3732 CALL grib_set(gaid, 'shapeOfTheEarth', 3)
3733 r2 = r1*(1.0d0 - f)
3734 CALL grib_set(gaid, 'scaleFactorOfEarthMajorAxis', 5)
3735 CALL grib_set(gaid, 'scaledValueOfEarthMajorAxis', &
3736 int(r1*100.0d0))
3737 CALL grib_set(gaid, 'scaleFactorOfEarthMinorAxis', 5)
3738 CALL grib_set(gaid, 'scaledValueOfEarthMinorAxis', &
3739 int(r2*100.0d0))
3740 ENDIF
3741 ENDIF
3742 END SELECT
3743
3744ELSE
3745
3746 IF (r1 == 6367470.0d0 .AND. f == 0.0d0) THEN ! spherical
3747 CALL grib_set(gaid, 'earthIsOblate', 0)
3748 ELSE IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
3749 CALL grib_set(gaid, 'earthIsOblate', 1)
3750 ELSE IF (f == 0.0d0) THEN ! generic spherical
3751 CALL grib_set(gaid, 'earthIsOblate', 0)
3752 CALL l4f_category_log(this%category,l4f_warn,'desired spherical Earth radius &
3753 &not supported in grib 1, coding with default radius of 6367470 m')
3754 ELSE ! generic ellipsoidal
3755 CALL grib_set(gaid, 'earthIsOblate', 1)
3756 CALL l4f_category_log(this%category,l4f_warn,'desired Earth ellipsoid &
3757 &not supported in grib 1, coding with default iau65 ellipsoid')
3758 ENDIF
3759
3760ENDIF
3761
3762END SUBROUTINE griddim_export_ellipsoid
3763
3764SUBROUTINE get_unproj_first(this, iScansNegatively, jScansPositively, &
3765 loFirst, laFirst)
3766TYPE(griddim_def),INTENT(in) :: this ! griddim object
3767INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
3768DOUBLE PRECISION,INTENT(out) :: loFirst, laFirst
3769
3770! compute lon and lat of first point from projected extremes
3771IF (iscansnegatively == 0) THEN
3772 IF (jscanspositively == 0) THEN
3773 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymax, lofirst, lafirst)
3774 ELSE
3775 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymin, lofirst, lafirst)
3776 ENDIF
3777ELSE
3778 IF (jscanspositively == 0) THEN
3779 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymax, lofirst, lafirst)
3780 ELSE
3781 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymin, lofirst, lafirst)
3782 ENDIF
3783ENDIF
3784! use the values kept for personal pleasure ?
3785! loFirst = this%grid%proj%polar%lon1
3786! laFirst = this%grid%proj%polar%lat1
3787END SUBROUTINE get_unproj_first
3788
3789SUBROUTINE get_unproj_last(this, iScansNegatively, jScansPositively, &
3790 loLast, laLast)
3791TYPE(griddim_def),INTENT(in) :: this ! griddim object
3792INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
3793DOUBLE PRECISION,INTENT(out) :: loLast, laLast
3794
3795! compute lon and lat of last point from projected extremes
3796IF (iscansnegatively == 0) THEN
3797 IF (jscanspositively == 0) THEN
3798 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymin, lolast, lalast)
3799 ELSE
3800 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymax, lolast, lalast)
3801 ENDIF
3802ELSE
3803 IF (jscanspositively == 0) THEN
3804 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymin, lolast, lalast)
3805 ELSE
3806 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymax, lolast, lalast)
3807 ENDIF
3808ENDIF
3809! use the values kept for personal pleasure ?
3810! loLast = this%grid%proj%polar%lon?
3811! laLast = this%grid%proj%polar%lat?
3812END SUBROUTINE get_unproj_last
3813
3814END SUBROUTINE griddim_export_gribapi
3815#endif
3816
3817
3818#ifdef HAVE_LIBGDAL
3819! gdal driver
3820SUBROUTINE griddim_import_gdal(this, gdalid, gdal_options)
3821USE gdal
3822TYPE(griddim_def),INTENT(inout) :: this ! griddim object
3823TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal rasterband pointer
3824TYPE(gdal_file_id_options),INTENT(in) :: gdal_options
3825
3826TYPE(gdaldataseth) :: hds
3827REAL(kind=c_double) :: geotrans(6) !, invgeotrans(6), x1, y1, x2, y2
3828INTEGER(kind=c_int) :: offsetx, offsety
3829INTEGER :: ier
3830
3831hds = gdalgetbanddataset(gdalid) ! go back to dataset
3832ier = gdalgetgeotransform(hds, geotrans)
3833
3834IF (ier /= 0) THEN
3835 CALL l4f_category_log(this%category, l4f_error, &
3836 'griddim_import_gdal: error in accessing gdal dataset')
3837 CALL raise_error()
3838 RETURN
3839ENDIF
3840IF (geotrans(3) /= 0.0_c_double .OR. geotrans(5) /= 0.0_c_double) THEN ! dataset has not a diagonal transformation
3841 CALL l4f_category_log(this%category, l4f_error, &
3842 'griddim_import_gdal: dataset has a non-diagonal transformation matrix, unsupported')
3843 CALL raise_error()
3844 RETURN
3845ENDIF
3846
3847CALL gdaldatasetbbsize_f(hds, gdal_options%xmin, gdal_options%ymin, &
3848 gdal_options%xmax, gdal_options%ymax, &
3849 this%dim%nx, this%dim%ny, offsetx, offsety, &
3850 this%grid%grid%xmin, this%grid%grid%ymin, this%grid%grid%xmax, this%grid%grid%ymax)
3851
3852IF (this%dim%nx == 0 .OR. this%dim%ny == 0) THEN
3853 CALL l4f_category_log(this%category, l4f_warn, &
3854 'griddim_import_gdal: requested bounding box '//t2c(gdal_options%xmin)//','// &
3855 t2c(gdal_options%ymin)//','//t2c(gdal_options%xmax)//','//&
3856 t2c(gdal_options%ymax))
3857 CALL l4f_category_log(this%category, l4f_warn, &
3858 'determines an empty dataset '//t2c(this%dim%nx)//'x'//t2c(this%dim%ny))
3859ENDIF
3860
3861! get grid corners
3862!CALL gdalapplygeotransform(geotrans, 0.5_c_double, 0.5_c_double, x1, y1)
3863!CALL gdalapplygeotransform(geotrans, gdalgetrasterbandxsize(gdalid)-0.5_c_double, &
3864! gdalgetrasterbandysize(gdalid)-0.5_c_double, x2, y2)
3865
3866!IF (geotrans(3) == 0.0_c_double .AND. geotrans(5) == 0.0_c_double) THEN ! transformation is diagonal, no transposing
3867! this%dim%nx = gdalgetrasterbandxsize(gdalid)
3868! this%dim%ny = gdalgetrasterbandysize(gdalid)
3869! this%grid%grid%xmin = MIN(x1, x2)
3870! this%grid%grid%xmax = MAX(x1, x2)
3871! this%grid%grid%ymin = MIN(y1, y2)
3872! this%grid%grid%ymax = MAX(y1, y2)
3873!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
3874!
3875! this%dim%nx = gdalgetrasterbandysize(gdalid)
3876! this%dim%ny = gdalgetrasterbandxsize(gdalid)
3877! this%grid%grid%xmin = MIN(y1, y2)
3878! this%grid%grid%xmax = MAX(y1, y2)
3879! this%grid%grid%ymin = MIN(x1, x2)
3880! this%grid%grid%ymax = MAX(x1, x2)
3881!
3882!ELSE ! transformation is a rotation, not supported
3883!ENDIF
3884
3885this%grid%proj%proj_type = 'regular_ll' ! forced, only one supported (check?)
3886
3887CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3888this%grid%grid%component_flag = 0
3889
3890END SUBROUTINE griddim_import_gdal
3891#endif
3892
3893
3894!> Display on the screen a brief content of the object.
3895SUBROUTINE griddim_display(this)
3896TYPE(griddim_def),INTENT(in) :: this !< griddim object to display
3897
3898print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
3899
3900CALL display(this%grid%proj)
3901CALL display(this%grid%grid)
3902CALL display(this%dim)
3903
3904print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
3905
3906END SUBROUTINE griddim_display
3907
3908
3909#define VOL7D_POLY_TYPE TYPE(grid_def)
3910#define VOL7D_POLY_TYPES _grid
3911#include "array_utilities_inc.F90"
3912#undef VOL7D_POLY_TYPE
3913#undef VOL7D_POLY_TYPES
3914
3915#define VOL7D_POLY_TYPE TYPE(griddim_def)
3916#define VOL7D_POLY_TYPES _griddim
3917#include "array_utilities_inc.F90"
3918#undef VOL7D_POLY_TYPE
3919#undef VOL7D_POLY_TYPES
3920
3921
3922!> Compute rotation matrix for wind unrotation. It allocates and
3923!! computes a matrix suitable for recomputing wind components in the
3924!! geographical system from the components in the projected
3925!! system. The rotation matrix \a rot_mat has to be passed as a
3926!! pointer and successively deallocated by the caller; it is a
3927!! 3-dimensional array where the first two dimensions are lon and lat
3928!! and the third, with extension 4, contains the packed rotation
3929!! matrix for the given grid point. It should work for every
3930!! projection. In order for the method to work, the \a griddim_unproj
3931!! method must have already been called for the \a griddim_def object.
3932!! \todo Check the algorithm and add some orthogonality tests.
3933SUBROUTINE griddim_wind_unrot(this, rot_mat)
3935TYPE(griddim_def), INTENT(in) :: this !< object describing the grid
3936DOUBLE PRECISION, POINTER :: rot_mat(:,:,:) !< rotation matrix for every grid point, to be deallocated by the caller, if \c .NOT. \c ASSOCIATED() an error occurred
3937
3938DOUBLE PRECISION :: dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j
3939DOUBLE PRECISION :: lat_factor
3940INTEGER :: i, j, ip1, im1, jp1, jm1;
3941
3942IF (this%dim%nx <= 0 .OR. this%dim%ny <= 0 .OR. &
3943 .NOT. ASSOCIATED(this%dim%lon) .OR. .NOT. ASSOCIATED(this%dim%lat)) THEN
3944 CALL l4f_category_log(this%category, l4f_error, 'rotate_uv must be called after correct init of griddim object')
3945 NULLIFY(rot_mat)
3946 RETURN
3947ENDIF
3948
3949ALLOCATE(rot_mat(this%dim%nx, this%dim%ny, 4))
3950
3951DO j = 1, this%dim%ny
3952 jp1 = min(j+1, this%dim%ny)
3953 jm1 = max(j-1, 1)
3954 DO i = 1, this%dim%nx
3955 ip1 = min(i+1, this%dim%nx)
3956 im1 = max(i-1, 1)
3957
3958 dlat_i = this%dim%lat(ip1,j) - this%dim%lat(im1,j)
3959 dlat_j = this%dim%lat(i,jp1) - this%dim%lat(i,jm1)
3960
3961 dlon_i = this%dim%lon(ip1,j) - this%dim%lon(im1,j)
3962! if ( dlon_i > pi ) dlon_i -= 2*pi;
3963! if ( dlon_i < (-pi) ) dlon_i += 2*pi;
3964 dlon_j = this%dim%lon(i,jp1) - this%dim%lon(i,jm1)
3965! if ( dlon_j > pi ) dlon_j -= 2*pi;
3966! if ( dlon_j < (-pi) ) dlon_j += 2*pi;
3967
3968! check whether this is really necessary !!!!
3969 lat_factor = cos(degrad*this%dim%lat(i,j))
3970 dlon_i = dlon_i * lat_factor
3971 dlon_j = dlon_j * lat_factor
3972
3973 dist_i = sqrt(dlon_i*dlon_i + dlat_i*dlat_i)
3974 dist_j = sqrt(dlon_j*dlon_j + dlat_j*dlat_j)
3975
3976 IF (dist_i > 0.d0) THEN
3977 rot_mat(i,j,1) = dlon_i/dist_i
3978 rot_mat(i,j,2) = dlat_i/dist_i
3979 ELSE
3980 rot_mat(i,j,1) = 0.0d0
3981 rot_mat(i,j,2) = 0.0d0
3982 ENDIF
3983 IF (dist_j > 0.d0) THEN
3984 rot_mat(i,j,3) = dlon_j/dist_j
3985 rot_mat(i,j,4) = dlat_j/dist_j
3986 ELSE
3987 rot_mat(i,j,3) = 0.0d0
3988 rot_mat(i,j,4) = 0.0d0
3989 ENDIF
3990
3991 ENDDO
3992ENDDO
3993
3994END SUBROUTINE griddim_wind_unrot
3995
3996
3997! compute zoom indices from geographical zoom coordinates
3998SUBROUTINE griddim_zoom_coord(this, ilon, ilat, flon, flat, ix, iy, fx, fy)
3999TYPE(griddim_def),INTENT(in) :: this
4000DOUBLE PRECISION,INTENT(in) :: ilon, ilat, flon, flat
4001INTEGER,INTENT(out) :: ix, iy, fx, fy
4002
4003DOUBLE PRECISION :: ix1, iy1, fx1, fy1
4004
4005! compute projected coordinates of vertices of desired lonlat rectangle
4006CALL proj(this, ilon, ilat, ix1, iy1)
4007CALL proj(this, flon, flat, fx1, fy1)
4008
4009CALL griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
4010
4011END SUBROUTINE griddim_zoom_coord
4012
4013
4014! compute zoom indices from projected zoom coordinates
4015SUBROUTINE griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
4016TYPE(griddim_def),INTENT(in) :: this
4017DOUBLE PRECISION,INTENT(in) :: ix1, iy1, fx1, fy1
4018INTEGER,INTENT(out) :: ix, iy, fx, fy
4019
4020INTEGER :: lix, liy, lfx, lfy
4021
4022! compute projected indices
4023lix = nint((ix1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
4024liy = nint((iy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
4025lfx = nint((fx1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
4026lfy = nint((fy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
4027! swap projected indices, in case grid is "strongly rotated"
4028ix = min(lix, lfx)
4029fx = max(lix, lfx)
4030iy = min(liy, lfy)
4031fy = max(liy, lfy)
4032
4033END SUBROUTINE griddim_zoom_projcoord
4034
4035
4036!> Reset a longitude value in the interval [0-360[.
4037!! The value is reset in place. This is usually useful in connection
4038!! with grib2 coding/decoding.
4039SUBROUTINE long_reset_0_360(lon)
4040DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
4041
4042IF (.NOT.c_e(lon)) RETURN
4043DO WHILE(lon < 0.0d0)
4044 lon = lon + 360.0d0
4045END DO
4046DO WHILE(lon >= 360.0d0)
4047 lon = lon - 360.0d0
4048END DO
4049
4050END SUBROUTINE long_reset_0_360
4051
4052
4053!> Reset a longitude value in the interval [-180-360[.
4054!! The value is reset in place. This is usually useful in connection
4055!! with grib1 coding/decoding.
4056SUBROUTINE long_reset_m180_360(lon)
4057DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
4058
4059IF (.NOT.c_e(lon)) RETURN
4060DO WHILE(lon < -180.0d0)
4061 lon = lon + 360.0d0
4062END DO
4063DO WHILE(lon >= 360.0d0)
4064 lon = lon - 360.0d0
4065END DO
4066
4067END SUBROUTINE long_reset_m180_360
4068
4069
4070!> Reset a longitude value in the interval [-90-270[.
4071!! The value is reset in place. This is usually useful in connection
4072!! with grib2 coding/decoding.
4073!SUBROUTINE long_reset_m90_270(lon)
4074!DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
4075!
4076!IF (.NOT.c_e(lon)) RETURN
4077!DO WHILE(lon < -90.0D0)
4078! lon = lon + 360.0D0
4079!END DO
4080!DO WHILE(lon >= 270.0D0)
4081! lon = lon - 360.0D0
4082!END DO
4083!
4084!END SUBROUTINE long_reset_m90_270
4085
4086
4087!> Reset a longitude value in the interval [-180-180[.
4088!! The value is reset in place. This is usually useful in connection
4089!! with grib2 coding/decoding.
4090SUBROUTINE long_reset_m180_180(lon)
4091DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
4092
4093IF (.NOT.c_e(lon)) RETURN
4094DO WHILE(lon < -180.0d0)
4095 lon = lon + 360.0d0
4096END DO
4097DO WHILE(lon >= 180.0d0)
4098 lon = lon - 360.0d0
4099END DO
4100
4101END SUBROUTINE long_reset_m180_180
4102
4103
4104SUBROUTINE long_reset_to_cart_closest(lon, lonref)
4105DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
4106DOUBLE PRECISION,INTENT(in) :: lonref !< the longitude to compare
4107
4108IF (.NOT.c_e(lon) .OR. .NOT.c_e(lonref)) RETURN
4109IF (abs(lon-lonref) < 180.0d0) RETURN ! nothing to do
4110lon = lon - nint((lon-lonref)/360.0d0)*360.0d0 ! se non e` vera e` ben trovata
4111
4112END SUBROUTINE long_reset_to_cart_closest
4113
4114
4115END MODULE grid_class
4116
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.