libsim Versione 7.1.11

◆ long_reset_m180_360()

subroutine long_reset_m180_360 ( double precision, intent(inout)  lon)
private

Reset a longitude value in the interval [-180-360[.

The value is reset in place. This is usually useful in connection with grib1 coding/decoding.

Parametri
[in,out]lonthe longitude to reset

Definizione alla linea 3121 del file grid_class.F90.

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