libsim Versione 7.1.11

◆ index_griddim()

integer function index_griddim ( type(griddim_def), dimension(:), intent(in)  vect,
type(griddim_def), intent(in)  search,
logical, dimension(:), intent(in), optional  mask,
logical, intent(in), optional  back,
integer, intent(in), optional  cache 
)

Cerca l'indice del primo o ultimo elemento di vect uguale a search.

Definizione alla linea 2909 del file grid_class.F90.

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