libsim Versione 7.2.4

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

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

Generated with Doxygen.