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