libsim Versione 7.1.11

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

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