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