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