libsim Versione 7.2.0

◆ 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 3115 del file grid_class.F90.

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

Generated with Doxygen.