libsim Versione 7.1.11

◆ map_inv_distinct_griddim()

integer function, dimension(dim) map_inv_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

map inv distinct

Definizione alla linea 2823 del file grid_class.F90.

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

Generated with Doxygen.