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