libsim Versione 7.1.11

◆ long_reset_0_360()

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

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

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

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