libsim Versione 7.2.0
|
◆ index_grid()
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
2440use geo_proj_class
2442use grid_rect_class
2448implicit none
2449
2450
2451character (len=255),parameter:: subcategory="grid_class"
2452
2453
2460 private
2461 type(geo_proj) :: proj
2462 type(grid_rect) :: grid
2463 integer :: category = 0
2465
2466
2473 type(grid_def) :: grid
2474 type(grid_dim) :: dim
2475 integer :: category = 0
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
2495 MODULE PROCEDURE griddim_init
2496END INTERFACE
2497
2500 MODULE PROCEDURE griddim_delete
2501END INTERFACE
2502
2505 MODULE PROCEDURE griddim_copy
2506END INTERFACE
2507
2511 MODULE PROCEDURE griddim_coord_proj!, griddim_proj
2512END INTERFACE
2513
2517 MODULE PROCEDURE griddim_coord_unproj, griddim_unproj
2518END INTERFACE
2519
2522 MODULE PROCEDURE griddim_get_val
2523END INTERFACE
2524
2527 MODULE PROCEDURE griddim_set_val
2528END INTERFACE
2529
2532 MODULE PROCEDURE griddim_write_unit
2533END INTERFACE
2534
2537 MODULE PROCEDURE griddim_read_unit
2538END INTERFACE
2539
2541INTERFACE import
2542 MODULE PROCEDURE griddim_import_grid_id
2543END INTERFACE
2544
2547 MODULE PROCEDURE griddim_export_grid_id
2548END INTERFACE
2549
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
2575 griddim_zoom_coord, griddim_zoom_projcoord, &
2579PUBLIC OPERATOR(==),OPERATOR(/=)
2580PUBLIC count_distinct, pack_distinct, count_and_pack_distinct, &
2581 map_distinct, map_inv_distinct,index
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
2654#endif
2655
2656END SUBROUTINE griddim_init
2657
2658
2660SUBROUTINE griddim_delete(this)
2661TYPE(griddim_def),INTENT(inout) :: this
2662
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
2681
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
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
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
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)
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
2811
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
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
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
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
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
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
2956IF (unit == geo_proj_unit_meter) THEN ! it is a plane projection
2958 IF (PRESENT(lonref)) THEN
2959 CALL long_reset_to_cart_closest(lov, lonref)
2961 ENDIF
2962
2963ELSE IF (unit == geo_proj_unit_degree) THEN ! it is a spheric projection
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)
2973! now reset rotated coordinates around zero
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
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
3102
3103#ifdef HAVE_LIBGRIBAPI
3104gaid = grid_id_get_gaid(ingrid_id)
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)
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
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)
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
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
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
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
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
3313 "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
3315 "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
3316#endif
3317
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
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
3368 "griddim_import_gribapi, Mercator grid orientation != 0 not supported")
3369 CALL raise_error()
3370 ENDIF
3371 ENDIF
3372
3373#ifdef 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)
3383 "griddim_import_gribapi, extremes from grib "//t2c(x1)//" "//t2c(y1))
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)
3400 ELSE
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
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
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
3479 CASE(2) ! iau65
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
3494 'read from grib, going on with spherical Earth but the results may be wrong')
3496 ELSE
3498 ENDIF
3499 CASE(4) ! iag-grs80
3501 CASE(5) ! wgs84
3503 CASE(6) ! spherical
3505! CASE(7) ! google earth-like?
3506 CASE default
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
3518 ELSE ! iau65
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
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
3659 "griddim_export_gribapi, error on x "//&
3660 trim(to_char(ex))//"/"//t2c(tol))
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
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
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)
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
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
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
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
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
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
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)
3936 ¬ supported in grib 1, coding with default radius of 6367470 m')
3937 ELSE ! generic ellipsoidal
3938 CALL grib_set(gaid, 'earthIsOblate', 1)
3940 ¬ 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
3957 ELSE
3959 ENDIF
3960ELSE
3961 IF (jscanspositively == 0) THEN
3963 ELSE
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
3982 ELSE
3984 ENDIF
3985ELSE
3986 IF (jscanspositively == 0) THEN
3988 ELSE
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
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
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
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))
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
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
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
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
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
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
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
4299
Method for testing the existence of the object. Definition: geo_proj_class.F90:268 Compute forward coordinate transformation from geographical system to projected system. Definition: grid_class.F90:308 Read the object from a formatted or unformatted file. Definition: grid_class.F90:334 Compute backward coordinate transformation from projected system to geographical system. Definition: grid_class.F90:314 Write the object on a formatted or unformatted file. Definition: grid_class.F90:329 Emit log message for a category with specific priority. Definition: log4fortran.F90:457 Module for describing geographically referenced regular grids. Definition: grid_class.F90:237 Module for defining the extension and coordinates of a rectangular georeferenced grid. Definition: grid_dim_class.F90:211 This module defines an abstract interface to different drivers for access to files containing gridded... Definition: grid_id_class.F90:249 Definitions of constants and functions for working with missing values. Definition: missing_values.f90:50 Module for quickly interpreting the OPTIONAL parameters passed to a subprogram. Definition: optional_values.f90:28 This object, mainly for internal use, describes a grid on a geographical projection,... Definition: grid_class.F90:257 This object completely describes a grid on a geographic projection. Definition: grid_class.F90:270 Derived type describing the extension of a grid and the geographical coordinates of each point. Definition: grid_dim_class.F90:221 |