libsim Versione 7.1.11

◆ pack_distinct_griddim()

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

compatta gli elementi distinti di vect in un array

Definizione alla linea 2578 del file grid_class.F90.

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