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