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