libsim Versione 7.2.4

◆ map_inv_distinct_grid()

integer function, dimension(dim) map_inv_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

map inv distinct

Definizione alla linea 2303 del file grid_class.F90.

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