libsim Versione 7.2.0

◆ count_distinct_grid()

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

conta gli elementi distinti in vect

Definizione alla linea 1980 del file grid_class.F90.

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