libsim Versione 7.1.11

◆ 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 1944 del file grid_class.F90.

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

Generated with Doxygen.