libsim Versione 7.1.11

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

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