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