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