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