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