libsim Versione 7.1.11

◆ griddim_wind_unrot()

subroutine griddim_wind_unrot ( type(griddim_def), intent(in)  this,
double precision, dimension(:,:,:), pointer  rot_mat 
)
private

Compute rotation matrix for wind unrotation.

It allocates and computes a matrix suitable for recomputing wind components in the geographical system from the components in the projected system. The rotation matrix rot_mat has to be passed as a pointer and successively deallocated by the caller; it is a 3-dimensional array where the first two dimensions are lon and lat and the third, with extension 4, contains the packed rotation matrix for the given grid point. It should work for every projection. In order for the method to work, the griddim_unproj method must have already been called for the griddim_def object.

Da fare:
Check the algorithm and add some orthogonality tests.
Parametri
[in]thisobject describing the grid
rot_matrotation matrix for every grid point, to be deallocated by the caller, if .NOT. ASSOCIATED() an error occurred

Definizione alla linea 2998 del file grid_class.F90.

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

Generated with Doxygen.