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