libsim Versione 7.2.4

◆ count_distinct_griddim()

integer function count_distinct_griddim ( type(griddim_def), dimension(:), intent(in) vect,
logical, dimension(:), intent(in), optional mask,
logical, intent(in), optional back )
private

conta gli elementi distinti in vect

Definizione alla linea 2494 del file grid_class.F90.

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