libsim Versione 7.1.11

◆ long_reset_m180_180()

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

Reset a longitude value in the interval [-90-270[.

The value is reset in place. This is usually useful in connection with grib2 coding/decoding. Reset a longitude value in the interval [-180-180[. The value is reset in place. This is usually useful in connection with grib2 coding/decoding.

Parametri
[in,out]lonthe longitude to reset

Definizione alla linea 3155 del file grid_class.F90.

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