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