libsim Versione 7.1.11

◆ map_distinct_grid()

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

map distinct

Definizione alla linea 2213 del file grid_class.F90.

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