libsim Versione 7.1.11

◆ map_distinct_griddim()

integer function, dimension(size(vect)) map_distinct_griddim ( type(griddim_def), dimension(:), intent(in)  vect,
logical, dimension(:), intent(in), optional  mask,
logical, intent(in), optional  back 
)
private

map distinct

Definizione alla linea 2727 del file grid_class.F90.

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

Generated with Doxygen.