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