libsim Versione 7.1.11

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

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