libsim Versione 7.2.1

◆ griddim_display()

subroutine griddim_display ( type(griddim_def), intent(in)  this)
private

Display on the screen a brief content of the object.

Parametri
[in]thisgriddim object to display

Definizione alla linea 1938 del file grid_class.F90.

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

Generated with Doxygen.