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