libsim Versione 7.2.4

◆ pack_distinct_grid()

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

compatta gli elementi distinti di vect in un array

Definizione alla linea 2058 del file grid_class.F90.

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

Generated with Doxygen.