libsim Versione 7.2.1

◆ pack_distinct_griddim()

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

compatta gli elementi distinti di vect in un array

Definizione alla linea 2572 del file grid_class.F90.

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

Generated with Doxygen.