libsim Versione 7.2.1

◆ griddim_wind_unrot()

subroutine griddim_wind_unrot ( type(griddim_def), intent(in) this,
double precision, dimension(:,:,:), pointer rot_mat )
private

Compute rotation matrix for wind unrotation.

It allocates and computes a matrix suitable for recomputing wind components in the geographical system from the components in the projected system. The rotation matrix rot_mat has to be passed as a pointer and successively deallocated by the caller; it is a 3-dimensional array where the first two dimensions are lon and lat and the third, with extension 4, contains the packed rotation matrix for the given grid point. It should work for every projection. In order for the method to work, the griddim_unproj method must have already been called for the griddim_def object.

Da fare
Check the algorithm and add some orthogonality tests.
Parametri
[in]thisobject describing the grid
rot_matrotation matrix for every grid point, to be deallocated by the caller, if .NOT. ASSOCIATED() an error occurred

Definizione alla linea 2992 del file grid_class.F90.

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

Generated with Doxygen.