libsim Versione 7.2.1

◆ map_distinct_griddim()

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

map distinct

Definizione alla linea 2721 del file grid_class.F90.

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

Generated with Doxygen.