libsim Versione 7.1.11

◆ long_reset_to_cart_closest()

subroutine long_reset_to_cart_closest ( double precision, intent(inout)  lon,
double precision, intent(in)  lonref 
)
private
Parametri
[in,out]lonthe longitude to reset
[in]lonrefthe longitude to compare

Definizione alla linea 3169 del file grid_class.F90.

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