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