libsim Versione 7.1.11

◆ 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 
)
private

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

Definizione alla linea 2395 del file grid_class.F90.

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

Generated with Doxygen.