libsim Versione 7.2.1

◆ index_grid()

integer function index_grid ( type(grid_def), dimension(:), intent(in) vect,
type(grid_def), intent(in) search,
logical, dimension(:), intent(in), optional mask,
logical, intent(in), optional back,
integer, intent(in), optional cache )

Cerca l'indice del primo o ultimo elemento di vect uguale a search.

Definizione alla linea 2389 del file grid_class.F90.

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

Generated with Doxygen.