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