libsim Versione 7.1.11
|
◆ map_distinct_grid()
map distinct Definizione alla linea 2213 del file grid_class.F90. 2214! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2215! authors:
2216! Davide Cesari <dcesari@arpa.emr.it>
2217! Paolo Patruno <ppatruno@arpa.emr.it>
2218
2219! This program is free software; you can redistribute it and/or
2220! modify it under the terms of the GNU General Public License as
2221! published by the Free Software Foundation; either version 2 of
2222! the License, or (at your option) any later version.
2223
2224! This program is distributed in the hope that it will be useful,
2225! but WITHOUT ANY WARRANTY; without even the implied warranty of
2226! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2227! GNU General Public License for more details.
2228
2229! You should have received a copy of the GNU General Public License
2230! along with this program. If not, see <http://www.gnu.org/licenses/>.
2231#include "config.h"
2232
2263use geo_proj_class
2265use grid_rect_class
2271implicit none
2272
2273
2274character (len=255),parameter:: subcategory="grid_class"
2275
2276
2283 private
2284 type(geo_proj) :: proj
2285 type(grid_rect) :: grid
2286 integer :: category = 0
2288
2289
2296 type(grid_def) :: grid
2297 type(grid_dim) :: dim
2298 integer :: category = 0
2300
2301
2305INTERFACE OPERATOR (==)
2306 MODULE PROCEDURE grid_eq, griddim_eq
2307END INTERFACE
2308
2312INTERFACE OPERATOR (/=)
2313 MODULE PROCEDURE grid_ne, griddim_ne
2314END INTERFACE
2315
2318 MODULE PROCEDURE griddim_init
2319END INTERFACE
2320
2323 MODULE PROCEDURE griddim_delete
2324END INTERFACE
2325
2328 MODULE PROCEDURE griddim_copy
2329END INTERFACE
2330
2334 MODULE PROCEDURE griddim_coord_proj!, griddim_proj
2335END INTERFACE
2336
2340 MODULE PROCEDURE griddim_coord_unproj, griddim_unproj
2341END INTERFACE
2342
2345 MODULE PROCEDURE griddim_get_val
2346END INTERFACE
2347
2350 MODULE PROCEDURE griddim_set_val
2351END INTERFACE
2352
2355 MODULE PROCEDURE griddim_write_unit
2356END INTERFACE
2357
2360 MODULE PROCEDURE griddim_read_unit
2361END INTERFACE
2362
2364INTERFACE import
2365 MODULE PROCEDURE griddim_import_grid_id
2366END INTERFACE
2367
2370 MODULE PROCEDURE griddim_export_grid_id
2371END INTERFACE
2372
2375 MODULE PROCEDURE griddim_display
2376END INTERFACE
2377
2378#define VOL7D_POLY_TYPE TYPE(grid_def)
2379#define VOL7D_POLY_TYPES _grid
2380#include "array_utilities_pre.F90"
2381#undef VOL7D_POLY_TYPE
2382#undef VOL7D_POLY_TYPES
2383
2384#define VOL7D_POLY_TYPE TYPE(griddim_def)
2385#define VOL7D_POLY_TYPES _griddim
2386#include "array_utilities_pre.F90"
2387#undef VOL7D_POLY_TYPE
2388#undef VOL7D_POLY_TYPES
2389
2390INTERFACE wind_unrot
2391 MODULE PROCEDURE griddim_wind_unrot
2392END INTERFACE
2393
2394
2395PRIVATE
2396
2398 griddim_zoom_coord, griddim_zoom_projcoord, &
2402PUBLIC OPERATOR(==),OPERATOR(/=)
2403PUBLIC count_distinct, pack_distinct, count_and_pack_distinct, &
2404 map_distinct, map_inv_distinct,index
2406PUBLIC griddim_central_lon, griddim_set_central_lon
2407CONTAINS
2408
2410SUBROUTINE griddim_init(this, nx, ny, &
2411 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
2412 proj_type, lov, zone, xoff, yoff, &
2413 longitude_south_pole, latitude_south_pole, angle_rotation, &
2414 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
2415 latin1, latin2, lad, projection_center_flag, &
2416 ellips_smaj_axis, ellips_flatt, ellips_type, &
2417 categoryappend)
2418TYPE(griddim_def),INTENT(inout) :: this
2419INTEGER,INTENT(in),OPTIONAL :: nx
2420INTEGER,INTENT(in),OPTIONAL :: ny
2421DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin
2422DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmax
2423DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymin
2424DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymax
2425DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx
2426DOUBLE PRECISION,INTENT(in),OPTIONAL :: dy
2429INTEGER,INTENT(in),OPTIONAL :: component_flag
2430CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
2431DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
2432INTEGER,INTENT(in),OPTIONAL :: zone
2433DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
2434DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
2435DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
2436DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
2437DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
2438DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
2439DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
2440DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
2441DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
2442DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
2443DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
2444INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
2445DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
2446DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
2447INTEGER,INTENT(in),OPTIONAL :: ellips_type
2448CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2449
2450CHARACTER(len=512) :: a_name
2451
2452IF (PRESENT(categoryappend)) THEN
2453 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
2454 trim(categoryappend))
2455ELSE
2456 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2457ENDIF
2458this%category=l4f_category_get(a_name)
2459
2460! geographical projection
2461this%grid%proj = geo_proj_new( &
2462 proj_type=proj_type, lov=lov, zone=zone, xoff=xoff, yoff=yoff, &
2463 longitude_south_pole=longitude_south_pole, &
2464 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
2465 longitude_stretch_pole=longitude_stretch_pole, &
2466 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
2467 latin1=latin1, latin2=latin2, lad=lad, projection_center_flag=projection_center_flag, &
2468 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, ellips_type=ellips_type)
2469! grid extension
2470this%grid%grid = grid_rect_new( &
2471 xmin, xmax, ymin, ymax, dx, dy, component_flag)
2472! grid size
2473this%dim = grid_dim_new(nx, ny)
2474
2475#ifdef DEBUG
2477#endif
2478
2479END SUBROUTINE griddim_init
2480
2481
2483SUBROUTINE griddim_delete(this)
2484TYPE(griddim_def),INTENT(inout) :: this
2485
2489
2490call l4f_category_delete(this%category)
2491
2492END SUBROUTINE griddim_delete
2493
2494
2496SUBROUTINE griddim_copy(this, that, categoryappend)
2497TYPE(griddim_def),INTENT(in) :: this
2498TYPE(griddim_def),INTENT(out) :: that
2499CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2500
2501CHARACTER(len=512) :: a_name
2502
2504
2508
2509! new category
2510IF (PRESENT(categoryappend)) THEN
2511 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
2512 trim(categoryappend))
2513ELSE
2514 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2515ENDIF
2516that%category=l4f_category_get(a_name)
2517
2518END SUBROUTINE griddim_copy
2519
2520
2523ELEMENTAL SUBROUTINE griddim_coord_proj(this, lon, lat, x, y)
2524TYPE(griddim_def),INTENT(in) :: this
2526DOUBLE PRECISION,INTENT(in) :: lon, lat
2528DOUBLE PRECISION,INTENT(out) :: x, y
2529
2531
2532END SUBROUTINE griddim_coord_proj
2533
2534
2537ELEMENTAL SUBROUTINE griddim_coord_unproj(this, x, y, lon, lat)
2538TYPE(griddim_def),INTENT(in) :: this
2540DOUBLE PRECISION,INTENT(in) :: x, y
2542DOUBLE PRECISION,INTENT(out) :: lon, lat
2543
2545
2546END SUBROUTINE griddim_coord_unproj
2547
2548
2549! Computes and sets the grid parameters required to compute
2550! coordinates of grid points in the projected system.
2551! probably meaningless
2552!SUBROUTINE griddim_proj(this)
2553!TYPE(griddim_def),INTENT(inout) :: this !< definition of projection and grid
2554!
2555!CALL proj(this, this%dim%lon(1,1), this%dim%lat(1,1), &
2556! this%grid%grid%xmin, this%grid%grid%ymin)
2557!
2558!CALL proj(this, this%dim%lon(this%dim%nx,this%dim%ny), &
2559! this%dim%lat(this%dim%nx,this%dim%ny), &
2560! this%grid%grid%xmax, this%grid%grid%ymax)
2561!
2562!END SUBROUTINE griddim_proj
2563
2571SUBROUTINE griddim_unproj(this)
2572TYPE(griddim_def),INTENT(inout) :: this
2573
2575CALL alloc(this%dim)
2576CALL griddim_unproj_internal(this)
2577
2578END SUBROUTINE griddim_unproj
2579
2580! internal subroutine needed for allocating automatic arrays
2581SUBROUTINE griddim_unproj_internal(this)
2582TYPE(griddim_def),INTENT(inout) ::this ! definition of projection and grid
2583
2584DOUBLE PRECISION :: x(this%dim%nx,this%dim%ny), y(this%dim%nx,this%dim%ny)
2585
2586CALL grid_rect_coordinates(this%grid%grid, x, y)
2588
2589END SUBROUTINE griddim_unproj_internal
2590
2591
2593SUBROUTINE griddim_get_val(this, nx, ny, &
2594 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
2595 proj, proj_type, lov, zone, xoff, yoff, &
2596 longitude_south_pole, latitude_south_pole, angle_rotation, &
2597 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
2598 latin1, latin2, lad, projection_center_flag, &
2599 ellips_smaj_axis, ellips_flatt, ellips_type)
2600TYPE(griddim_def),INTENT(in) :: this
2601INTEGER,INTENT(out),OPTIONAL :: nx
2602INTEGER,INTENT(out),OPTIONAL :: ny
2604DOUBLE PRECISION,INTENT(out),OPTIONAL :: xmin, xmax, ymin, ymax
2606DOUBLE PRECISION,INTENT(out),OPTIONAL :: dx, dy
2609INTEGER,INTENT(out),OPTIONAL :: component_flag
2610TYPE(geo_proj),INTENT(out),OPTIONAL :: proj
2611CHARACTER(len=*),INTENT(out),OPTIONAL :: proj_type
2612DOUBLE PRECISION,INTENT(out),OPTIONAL :: lov
2613INTEGER,INTENT(out),OPTIONAL :: zone
2614DOUBLE PRECISION,INTENT(out),OPTIONAL :: xoff
2615DOUBLE PRECISION,INTENT(out),OPTIONAL :: yoff
2616DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_south_pole
2617DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_south_pole
2618DOUBLE PRECISION,INTENT(out),OPTIONAL :: angle_rotation
2619DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_stretch_pole
2620DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_stretch_pole
2621DOUBLE PRECISION,INTENT(out),OPTIONAL :: stretch_factor
2622DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin1
2623DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin2
2624DOUBLE PRECISION,INTENT(out),OPTIONAL :: lad
2625INTEGER,INTENT(out),OPTIONAL :: projection_center_flag
2626DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_smaj_axis
2627DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_flatt
2628INTEGER,INTENT(out),OPTIONAL :: ellips_type
2629
2630IF (PRESENT(nx)) nx = this%dim%nx
2631IF (PRESENT(ny)) ny = this%dim%ny
2632
2634
2636 xoff=xoff, yoff=yoff, &
2637 longitude_south_pole=longitude_south_pole, &
2638 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
2639 longitude_stretch_pole=longitude_stretch_pole, &
2640 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
2641 latin1=latin1, latin2=latin2, lad=lad, &
2642 projection_center_flag=projection_center_flag, &
2643 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
2644 ellips_type=ellips_type)
2645
2647 xmin, xmax, ymin, ymax, dx, dy, component_flag)
2648
2649END SUBROUTINE griddim_get_val
2650
2651
2653SUBROUTINE griddim_set_val(this, nx, ny, &
2654 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
2655 proj_type, lov, zone, xoff, yoff, &
2656 longitude_south_pole, latitude_south_pole, angle_rotation, &
2657 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
2658 latin1, latin2, lad, projection_center_flag, &
2659 ellips_smaj_axis, ellips_flatt, ellips_type)
2660TYPE(griddim_def),INTENT(inout) :: this
2661INTEGER,INTENT(in),OPTIONAL :: nx
2662INTEGER,INTENT(in),OPTIONAL :: ny
2664DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin, xmax, ymin, ymax
2666DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx, dy
2669INTEGER,INTENT(in),OPTIONAL :: component_flag
2670CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
2671DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
2672INTEGER,INTENT(in),OPTIONAL :: zone
2673DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
2674DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
2675DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
2676DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
2677DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
2678DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
2679DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
2680DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
2681DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
2682DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
2683DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
2684INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
2685DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
2686DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
2687INTEGER,INTENT(in),OPTIONAL :: ellips_type
2688
2689IF (PRESENT(nx)) this%dim%nx = nx
2690IF (PRESENT(ny)) this%dim%ny = ny
2691
2693 xoff=xoff, yoff=yoff, longitude_south_pole=longitude_south_pole, &
2694 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
2695 longitude_stretch_pole=longitude_stretch_pole, &
2696 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
2697 latin1=latin1, latin2=latin2, lad=lad, &
2698 projection_center_flag=projection_center_flag, &
2699 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
2700 ellips_type=ellips_type)
2701
2703 xmin, xmax, ymin, ymax, dx, dy, component_flag)
2704
2705END SUBROUTINE griddim_set_val
2706
2707
2712SUBROUTINE griddim_read_unit(this, unit)
2713TYPE(griddim_def),INTENT(out) :: this
2714INTEGER, INTENT(in) :: unit
2715
2716
2720
2721END SUBROUTINE griddim_read_unit
2722
2723
2728SUBROUTINE griddim_write_unit(this, unit)
2729TYPE(griddim_def),INTENT(in) :: this
2730INTEGER, INTENT(in) :: unit
2731
2732
2736
2737END SUBROUTINE griddim_write_unit
2738
2739
2743FUNCTION griddim_central_lon(this) RESULT(lon)
2744TYPE(griddim_def),INTENT(inout) :: this
2745
2746DOUBLE PRECISION :: lon
2747
2748CALL griddim_pistola_central_lon(this, lon)
2749
2750END FUNCTION griddim_central_lon
2751
2752
2756SUBROUTINE griddim_set_central_lon(this, lonref)
2757TYPE(griddim_def),INTENT(inout) :: this
2758DOUBLE PRECISION,INTENT(in) :: lonref
2759
2760DOUBLE PRECISION :: lon
2761
2762CALL griddim_pistola_central_lon(this, lon, lonref)
2763
2764END SUBROUTINE griddim_set_central_lon
2765
2766
2767! internal subroutine for performing tasks common to the prevous two
2768SUBROUTINE griddim_pistola_central_lon(this, lon, lonref)
2769TYPE(griddim_def),INTENT(inout) :: this ! grid descriptor
2770DOUBLE PRECISION,INTENT(inout) :: lon ! central longitude
2771DOUBLE PRECISION,INTENT(in),OPTIONAL :: lonref ! reference longitude
2772
2773INTEGER :: unit
2774DOUBLE PRECISION :: lonsp, latsp, londelta, lov, lonrot
2775CHARACTER(len=80) :: ptype
2776
2777lon = dmiss
2779IF (unit == geo_proj_unit_meter) THEN ! it is a plane projection
2781 IF (PRESENT(lonref)) THEN
2782 CALL long_reset_to_cart_closest(lov, lonref)
2784 ENDIF
2785
2786ELSE IF (unit == geo_proj_unit_degree) THEN ! it is a spheric projection
2788 longitude_south_pole=lonsp, latitude_south_pole=latsp)
2789 SELECT CASE(ptype)
2790 CASE('rotated_ll','stretched_rotated_ll') ! use origin of rotated system
2791 IF (latsp < 0.0d0) THEN
2792 lon = lonsp
2793 IF (PRESENT(lonref)) THEN
2794 CALL long_reset_to_cart_closest(lon, lonref)
2796! now reset rotated coordinates around zero
2798 lonrot = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
2799 ENDIF
2800 londelta = lonrot
2801 CALL long_reset_to_cart_closest(londelta, 0.0d0)
2802 londelta = londelta - lonrot
2803 this%grid%grid%xmin = this%grid%grid%xmin + londelta
2804 this%grid%grid%xmax = this%grid%grid%xmax + londelta
2805 ENDIF
2806 ELSE ! this part to be checked
2807 lon = modulo(lonsp + 180.0d0, 360.0d0)
2808! IF (PRESENT(lonref)) THEN
2809! CALL long_reset_to_cart_closest(lov, lonref)
2810! CALL set_val(this%grid%proj, longitude_south_pole=lonref)
2811! ENDIF
2812 ENDIF
2813 CASE default ! use real grid limits
2815 lon = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
2816 ENDIF
2817 IF (PRESENT(lonref)) THEN
2818 londelta = lon
2819 CALL long_reset_to_cart_closest(londelta, lonref)
2820 londelta = londelta - lon
2821 this%grid%grid%xmin = this%grid%grid%xmin + londelta
2822 this%grid%grid%xmax = this%grid%grid%xmax + londelta
2823 ENDIF
2824 END SELECT
2825ENDIF
2826
2827END SUBROUTINE griddim_pistola_central_lon
2828
2829
2833SUBROUTINE griddim_gen_coord(this, x, y)
2834TYPE(griddim_def),INTENT(in) :: this
2835DOUBLE PRECISION,INTENT(out) :: x(:,:)
2836DOUBLE PRECISION,INTENT(out) :: y(:,:)
2837
2838
2839CALL grid_rect_coordinates(this%grid%grid, x, y)
2840
2841END SUBROUTINE griddim_gen_coord
2842
2843
2845SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
2846TYPE(griddim_def), INTENT(in) :: this
2847INTEGER,INTENT(in) :: nx
2848INTEGER,INTENT(in) :: ny
2849DOUBLE PRECISION,INTENT(out) :: dx
2850DOUBLE PRECISION,INTENT(out) :: dy
2851
2852CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
2853
2854END SUBROUTINE griddim_steps
2855
2856
2858SUBROUTINE griddim_setsteps(this)
2859TYPE(griddim_def), INTENT(inout) :: this
2860
2861CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
2862
2863END SUBROUTINE griddim_setsteps
2864
2865
2866! TODO
2867! bisogna sviluppare gli altri operatori
2868ELEMENTAL FUNCTION grid_eq(this, that) RESULT(res)
2869TYPE(grid_def),INTENT(IN) :: this, that
2870LOGICAL :: res
2871
2872res = this%proj == that%proj .AND. &
2873 this%grid == that%grid
2874
2875END FUNCTION grid_eq
2876
2877
2878ELEMENTAL FUNCTION griddim_eq(this, that) RESULT(res)
2879TYPE(griddim_def),INTENT(IN) :: this, that
2880LOGICAL :: res
2881
2882res = this%grid == that%grid .AND. &
2883 this%dim == that%dim
2884
2885END FUNCTION griddim_eq
2886
2887
2888ELEMENTAL FUNCTION grid_ne(this, that) RESULT(res)
2889TYPE(grid_def),INTENT(IN) :: this, that
2890LOGICAL :: res
2891
2892res = .NOT.(this == that)
2893
2894END FUNCTION grid_ne
2895
2896
2897ELEMENTAL FUNCTION griddim_ne(this, that) RESULT(res)
2898TYPE(griddim_def),INTENT(IN) :: this, that
2899LOGICAL :: res
2900
2901res = .NOT.(this == that)
2902
2903END FUNCTION griddim_ne
2904
2905
2911SUBROUTINE griddim_import_grid_id(this, ingrid_id)
2912#ifdef HAVE_LIBGDAL
2913USE gdal
2914#endif
2915TYPE(griddim_def),INTENT(inout) :: this
2916TYPE(grid_id),INTENT(in) :: ingrid_id
2917
2918#ifdef HAVE_LIBGRIBAPI
2919INTEGER :: gaid
2920#endif
2921#ifdef HAVE_LIBGDAL
2922TYPE(gdalrasterbandh) :: gdalid
2923#endif
2925
2926#ifdef HAVE_LIBGRIBAPI
2927gaid = grid_id_get_gaid(ingrid_id)
2929#endif
2930#ifdef HAVE_LIBGDAL
2931gdalid = grid_id_get_gdalid(ingrid_id)
2932IF (gdalassociated(gdalid)) CALL griddim_import_gdal(this, gdalid, &
2933 grid_id_get_gdal_options(ingrid_id))
2934#endif
2935
2936END SUBROUTINE griddim_import_grid_id
2937
2938
2943SUBROUTINE griddim_export_grid_id(this, outgrid_id)
2944#ifdef HAVE_LIBGDAL
2945USE gdal
2946#endif
2947TYPE(griddim_def),INTENT(in) :: this
2948TYPE(grid_id),INTENT(inout) :: outgrid_id
2949
2950#ifdef HAVE_LIBGRIBAPI
2951INTEGER :: gaid
2952#endif
2953#ifdef HAVE_LIBGDAL
2954TYPE(gdalrasterbandh) :: gdalid
2955#endif
2956
2957#ifdef HAVE_LIBGRIBAPI
2958gaid = grid_id_get_gaid(outgrid_id)
2960#endif
2961#ifdef HAVE_LIBGDAL
2962gdalid = grid_id_get_gdalid(outgrid_id)
2963!IF (gdalassociated(gdalid)
2964! export for gdal not implemented, log?
2965#endif
2966
2967END SUBROUTINE griddim_export_grid_id
2968
2969
2970#ifdef HAVE_LIBGRIBAPI
2971! grib_api driver
2972SUBROUTINE griddim_import_gribapi(this, gaid)
2973USE grib_api
2974TYPE(griddim_def),INTENT(inout) :: this ! griddim object
2975INTEGER, INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
2976
2977DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, x1, y1, orient
2978INTEGER :: EditionNumber, iScansNegatively, jScansPositively, zone, datum, &
2979 reflon, ierr
2980
2981! Generic keys
2982CALL grib_get(gaid, 'typeOfGrid', this%grid%proj%proj_type)
2983#ifdef DEBUG
2985 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type))
2986#endif
2987CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
2988
2989IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
2990 CALL grib_get(gaid,'numberOfDataPoints', this%dim%nx)
2991 this%dim%ny = 1
2992 this%grid%grid%component_flag = 0
2993 CALL griddim_import_ellipsoid(this, gaid)
2994 RETURN
2995ENDIF
2996
2997! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
2998CALL grib_get(gaid, 'Ni', this%dim%nx)
2999CALL grib_get(gaid, 'Nj', this%dim%ny)
3000CALL griddim_import_ellipsoid(this, gaid) ! placed here, not valid for utm datum /= 1
3001
3002CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
3003CALL grib_get(gaid,'jScansPositively',jscanspositively)
3004CALL grib_get(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
3005
3006! Keys for rotated grids (checked through missing values)
3007CALL grib_get_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
3008 this%grid%proj%rotated%longitude_south_pole)
3009CALL grib_get_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
3010 this%grid%proj%rotated%latitude_south_pole)
3011CALL grib_get_dmiss(gaid,'angleOfRotationInDegrees', &
3012 this%grid%proj%rotated%angle_rotation)
3013
3014! Keys for stretched grids (checked through missing values)
3015! units must be verified, still experimental in grib_api
3016! # TODO: Is it a float? Is it signed?
3017IF (editionnumber == 1) THEN
3018 CALL grib_get_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
3019 this%grid%proj%stretched%longitude_stretch_pole)
3020 CALL grib_get_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
3021 this%grid%proj%stretched%latitude_stretch_pole)
3022 CALL grib_get_dmiss(gaid,'stretchingFactor', &
3023 this%grid%proj%stretched%stretch_factor)
3024ELSE IF (editionnumber == 2) THEN
3025 CALL grib_get_dmiss(gaid,'longitudeOfThePoleOfStretching', &
3026 this%grid%proj%stretched%longitude_stretch_pole)
3027 CALL grib_get_dmiss(gaid,'latitudeOfThePoleOfStretching', &
3028 this%grid%proj%stretched%latitude_stretch_pole)
3029 CALL grib_get_dmiss(gaid,'stretchingFactor', &
3030 this%grid%proj%stretched%stretch_factor)
3032 this%grid%proj%stretched%stretch_factor = this%grid%proj%stretched%stretch_factor*1.0d-6
3033ENDIF
3034
3035! Projection-dependent keys
3036SELECT CASE (this%grid%proj%proj_type)
3037
3038! Keys for spherical coordinate systems
3039CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
3040
3041 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3042 CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3043 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3044 CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3045
3046! longitudes are sometimes wrongly coded even in grib2 and even by the
3047! Metoffice!
3048! longitudeOfFirstGridPointInDegrees = 354.911;
3049! longitudeOfLastGridPointInDegrees = 363.311;
3050 CALL long_reset_0_360(lofirst)
3051 CALL long_reset_0_360(lolast)
3052
3053 IF (iscansnegatively == 0) THEN
3054 this%grid%grid%xmin = lofirst
3055 this%grid%grid%xmax = lolast
3056 ELSE
3057 this%grid%grid%xmax = lofirst
3058 this%grid%grid%xmin = lolast
3059 ENDIF
3060 IF (jscanspositively == 0) THEN
3061 this%grid%grid%ymax = lafirst
3062 this%grid%grid%ymin = lalast
3063 ELSE
3064 this%grid%grid%ymin = lafirst
3065 this%grid%grid%ymax = lalast
3066 ENDIF
3067
3068! reset longitudes in order to have a Cartesian plane
3069 IF (this%grid%grid%xmax < this%grid%grid%xmin) &
3070 this%grid%grid%xmin = this%grid%grid%xmin - 360.d0
3071
3072! compute dx and dy (should we get them from grib?)
3073 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3074
3075! Keys for polar projections
3076CASE ('polar_stereographic', 'lambert', 'albers')
3077
3078 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
3079 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
3080! latin1/latin2 may be missing (e.g. stereographic)
3081 CALL grib_get_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
3082 CALL grib_get_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
3083#ifdef DEBUG
3085 "griddim_import_gribapi, latin1/2 "// &
3086 trim(to_char(this%grid%proj%polar%latin1))//" "// &
3087 trim(to_char(this%grid%proj%polar%latin2)))
3088#endif
3089! projection center flag, aka hemisphere
3090 CALL grib_get(gaid,'projectionCenterFlag',&
3091 this%grid%proj%polar%projection_center_flag, ierr)
3092 IF (ierr /= grib_success) THEN ! try center/centre
3093 CALL grib_get(gaid,'projectionCentreFlag',&
3094 this%grid%proj%polar%projection_center_flag)
3095 ENDIF
3096
3097 IF (iand(this%grid%proj%polar%projection_center_flag,64) == 1) THEN
3099 "griddim_import_gribapi, bi-polar projections not supported")
3100 CALL raise_error()
3101 ENDIF
3102! line of view, aka central meridian
3103 CALL grib_get(gaid,'LoVInDegrees',this%grid%proj%lov)
3104#ifdef DEBUG
3106 "griddim_import_gribapi, central meridian "//trim(to_char(this%grid%proj%lov)))
3107#endif
3108
3109! latitude at which dx and dy are valid
3110 IF (editionnumber == 1) THEN
3111! ECMWF (gribex/grib_api) says: Grid lengths are in metres, at the
3112! 60-degree parallel nearest to the pole on the projection plane.
3113! IF (IAND(this%projection_center_flag, 128) == 0) THEN
3114! this%grid%proj%polar%lad = 60.D0
3115! ELSE
3116! this%grid%proj%polar%lad = -60.D0
3117! ENDIF
3118! WMO says: Grid lengths are in units of metres, at the secant cone
3119! intersection parallel nearest to the pole on the projection plane.
3120 this%grid%proj%polar%lad = this%grid%proj%polar%latin1
3121 ELSE IF (editionnumber == 2) THEN
3122 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3123 ENDIF
3124#ifdef DEBUG
3126 "griddim_import_gribapi, lad "//trim(to_char(this%grid%proj%polar%lad)))
3127#endif
3128
3129! compute projected extremes from lon and lat of first point
3130 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3131 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3132 CALL long_reset_0_360(lofirst)
3133 CALL long_reset_to_cart_closest(this%grid%proj%lov, lofirst)
3134#ifdef DEBUG
3136 "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
3138 "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
3139#endif
3140
3142 IF (iscansnegatively == 0) THEN
3143 this%grid%grid%xmin = x1
3144 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
3145 ELSE
3146 this%grid%grid%xmax = x1
3147 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
3148 ENDIF
3149 IF (jscanspositively == 0) THEN
3150 this%grid%grid%ymax = y1
3151 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
3152 ELSE
3153 this%grid%grid%ymin = y1
3154 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
3155 ENDIF
3156! keep these values for personal pleasure
3157 this%grid%proj%polar%lon1 = lofirst
3158 this%grid%proj%polar%lat1 = lafirst
3159
3160! Keys for equatorial projections
3161CASE ('mercator')
3162
3163 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
3164 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
3165 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3166 this%grid%proj%lov = 0.0d0 ! not in grib
3167
3168 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3169 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3170
3172 IF (iscansnegatively == 0) THEN
3173 this%grid%grid%xmin = x1
3174 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
3175 ELSE
3176 this%grid%grid%xmax = x1
3177 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
3178 ENDIF
3179 IF (jscanspositively == 0) THEN
3180 this%grid%grid%ymax = y1
3181 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
3182 ELSE
3183 this%grid%grid%ymin = y1
3184 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
3185 ENDIF
3186
3187 IF (editionnumber == 2) THEN
3188 CALL grib_get(gaid,'orientationOfTheGridInDegrees',orient)
3189 IF (orient /= 0.0d0) THEN
3191 "griddim_import_gribapi, Mercator grid orientation != 0 not supported")
3192 CALL raise_error()
3193 ENDIF
3194 ENDIF
3195
3196#ifdef DEBUG
3199 "griddim_import_gribapi, unprojected first point "//t2c(lofirst)//" "// &
3200 t2c(lafirst))
3201
3202 CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3203 CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3206 "griddim_import_gribapi, extremes from grib "//t2c(x1)//" "//t2c(y1))
3208 "griddim_import_gribapi, extremes from proj "//t2c(this%grid%grid%xmin)// &
3209 " "//t2c(this%grid%grid%ymin)//" "//t2c(this%grid%grid%xmax)//" "// &
3210 t2c(this%grid%grid%ymax))
3211#endif
3212
3213CASE ('UTM')
3214
3215 CALL grib_get(gaid,'zone',zone)
3216
3217 CALL grib_get(gaid,'datum',datum)
3218 IF (datum == 0) THEN
3219 CALL grib_get(gaid,'referenceLongitude',reflon)
3220 CALL grib_get(gaid,'falseEasting',this%grid%proj%xoff)
3221 CALL grib_get(gaid,'falseNorthing',this%grid%proj%yoff)
3223 ELSE
3225 CALL raise_fatal_error()
3226 ENDIF
3227
3228 CALL grib_get(gaid,'eastingOfFirstGridPoint',lofirst)
3229 CALL grib_get(gaid,'eastingOfLastGridPoint',lolast)
3230 CALL grib_get(gaid,'northingOfFirstGridPoint',lafirst)
3231 CALL grib_get(gaid,'northingOfLastGridPoint',lalast)
3232
3233 IF (iscansnegatively == 0) THEN
3234 this%grid%grid%xmin = lofirst
3235 this%grid%grid%xmax = lolast
3236 ELSE
3237 this%grid%grid%xmax = lofirst
3238 this%grid%grid%xmin = lolast
3239 ENDIF
3240 IF (jscanspositively == 0) THEN
3241 this%grid%grid%ymax = lafirst
3242 this%grid%grid%ymin = lalast
3243 ELSE
3244 this%grid%grid%ymin = lafirst
3245 this%grid%grid%ymax = lalast
3246 ENDIF
3247
3248! compute dx and dy (should we get them from grib?)
3249 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3250
3251CASE default
3253 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
3254 CALL raise_error()
3255
3256END SELECT
3257
3258CONTAINS
3259! utilities routines for grib_api, is there a better place?
3260SUBROUTINE grib_get_dmiss(gaid, key, value)
3261INTEGER,INTENT(in) :: gaid
3262CHARACTER(len=*),INTENT(in) :: key
3263DOUBLE PRECISION,INTENT(out) :: value
3264
3265INTEGER :: ierr
3266
3267CALL grib_get(gaid, key, value, ierr)
3268IF (ierr /= grib_success) value = dmiss
3269
3270END SUBROUTINE grib_get_dmiss
3271
3272SUBROUTINE grib_get_imiss(gaid, key, value)
3273INTEGER,INTENT(in) :: gaid
3274CHARACTER(len=*),INTENT(in) :: key
3275INTEGER,INTENT(out) :: value
3276
3277INTEGER :: ierr
3278
3279CALL grib_get(gaid, key, value, ierr)
3280IF (ierr /= grib_success) value = imiss
3281
3282END SUBROUTINE grib_get_imiss
3283
3284
3285SUBROUTINE griddim_import_ellipsoid(this, gaid)
3286TYPE(griddim_def),INTENT(inout) :: this
3287INTEGER,INTENT(in) :: gaid
3288
3289INTEGER :: shapeofearth, iv, is
3290DOUBLE PRECISION :: r1, r2
3291
3292IF (editionnumber == 2) THEN
3293 CALL grib_get(gaid, 'shapeOfTheEarth', shapeofearth)
3294 SELECT CASE(shapeofearth)
3295 CASE(0) ! spherical
3297 CASE(1) ! spherical generic
3298 CALL grib_get(gaid, 'scaleFactorOfRadiusOfSphericalEarth', is)
3299 CALL grib_get(gaid, 'scaledValueOfRadiusOfSphericalEarth', iv)
3300 r1 = dble(iv) / 10**is
3302 CASE(2) ! iau65
3304 CASE(3,7) ! ellipsoidal generic
3305 CALL grib_get(gaid, 'scaleFactorOfEarthMajorAxis', is)
3306 CALL grib_get(gaid, 'scaledValueOfEarthMajorAxis', iv)
3307 r1 = dble(iv) / 10**is
3308 CALL grib_get(gaid, 'scaleFactorOfEarthMinorAxis', is)
3309 CALL grib_get(gaid, 'scaledValueOfEarthMinorAxis', iv)
3310 r2 = dble(iv) / 10**is
3311 IF (shapeofearth == 3) THEN ! km->m
3312 r1 = r1*1000.0d0
3313 r2 = r2*1000.0d0
3314 ENDIF
3315 IF (abs(r1) < 1.0d-6) THEN ! suspicious data read from grib
3317 'read from grib, going on with spherical Earth but the results may be wrong')
3319 ELSE
3321 ENDIF
3322 CASE(4) ! iag-grs80
3324 CASE(5) ! wgs84
3326 CASE(6) ! spherical
3328! CASE(7) ! google earth-like?
3329 CASE default
3331 t2c(shapeofearth)//' not supported in grib2')
3332 CALL raise_fatal_error()
3333
3334 END SELECT
3335
3336ELSE
3337
3338 CALL grib_get(gaid, 'earthIsOblate', shapeofearth)
3339 IF (shapeofearth == 0) THEN ! spherical
3341 ELSE ! iau65
3343 ENDIF
3344
3345ENDIF
3346
3347END SUBROUTINE griddim_import_ellipsoid
3348
3349
3350END SUBROUTINE griddim_import_gribapi
3351
3352
3353! grib_api driver
3354SUBROUTINE griddim_export_gribapi(this, gaid)
3355USE grib_api
3356TYPE(griddim_def),INTENT(in) :: this ! griddim object
3357INTEGER, INTENT(inout) :: gaid ! grib_api id of the grib loaded in memory to export
3358
3359INTEGER :: EditionNumber, iScansNegatively, jScansPositively, nv, pvl, zone, ierr
3360DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, reflon
3361DOUBLE PRECISION :: sdx, sdy, ratio, tol, ex, ey
3362
3363
3364! Generic keys
3365CALL grib_get(gaid,'GRIBEditionNumber', editionnumber)
3366! the following required since eccodes
3367IF (editionnumber == 2) CALL grib_set(gaid,'shapeOfTheEarth', 6)
3368CALL grib_set(gaid,'typeOfGrid', this%grid%proj%proj_type)
3369#ifdef DEBUG
3371 "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
3372#endif
3373
3374IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
3375 CALL grib_set(gaid,'numberOfDataPoints', this%dim%nx)
3376! reenable when possible
3377! CALL griddim_export_ellipsoid(this, gaid)
3378 RETURN
3379ENDIF
3380
3381
3382! Edition dependent setup
3383IF (editionnumber == 1) THEN
3384 ratio = 1.d3
3385ELSE IF (editionnumber == 2) THEN
3386 ratio = 1.d6
3387ELSE
3388 ratio = 0.0d0 ! signal error?!
3389ENDIF
3390
3391! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
3392CALL griddim_export_ellipsoid(this, gaid)
3393CALL grib_set(gaid, 'Ni', this%dim%nx)
3394CALL grib_set(gaid, 'Nj', this%dim%ny)
3395CALL grib_set(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
3396
3397CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
3398CALL grib_get(gaid,'jScansPositively',jscanspositively)
3399
3400! Keys for rotated grids (checked through missing values and/or error code)
3401!SELECT CASE (this%grid%proj%proj_type)
3402!CASE ('rotated_ll', 'stretched_rotated_ll', 'polar_stereographic', 'lambert', 'albers')
3403CALL grib_set_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
3404 this%grid%proj%rotated%longitude_south_pole, 0.0d0)
3405CALL grib_set_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
3406 this%grid%proj%rotated%latitude_south_pole, -90.0d0)
3407IF (editionnumber == 1) THEN
3408 CALL grib_set_dmiss(gaid,'angleOfRotationInDegrees', &
3409 this%grid%proj%rotated%angle_rotation, 0.0d0)
3410ELSE IF (editionnumber == 2)THEN
3411 CALL grib_set_dmiss(gaid,'angleOfRotationOfProjectionInDegrees', &
3412 this%grid%proj%rotated%angle_rotation, 0.0d0)
3413ENDIF
3414
3415! Keys for stretched grids (checked through missing values and/or error code)
3416! units must be verified, still experimental in grib_api
3417! # TODO: Is it a float? Is it signed?
3418IF (editionnumber == 1) THEN
3419 CALL grib_set_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
3420 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
3421 CALL grib_set_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
3422 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
3423 CALL grib_set_dmiss(gaid,'stretchingFactor', &
3424 this%grid%proj%stretched%stretch_factor, 1.0d0)
3425ELSE IF (editionnumber == 2) THEN
3426 CALL grib_set_dmiss(gaid,'longitudeOfThePoleOfStretching', &
3427 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
3428 CALL grib_set_dmiss(gaid,'latitudeOfThePoleOfStretching', &
3429 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
3430 CALL grib_set_dmiss(gaid,'stretchingFactor', &
3431 this%grid%proj%stretched%stretch_factor, 1.0d6, 1.0d6)
3432ENDIF
3433
3434! Projection-dependent keys
3435SELECT CASE (this%grid%proj%proj_type)
3436
3437! Keys for sphaerical coordinate systems
3438CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
3439
3440 IF (iscansnegatively == 0) THEN
3441 lofirst = this%grid%grid%xmin
3442 lolast = this%grid%grid%xmax
3443 ELSE
3444 lofirst = this%grid%grid%xmax
3445 lolast = this%grid%grid%xmin
3446 ENDIF
3447 IF (jscanspositively == 0) THEN
3448 lafirst = this%grid%grid%ymax
3449 lalast = this%grid%grid%ymin
3450 ELSE
3451 lafirst = this%grid%grid%ymin
3452 lalast = this%grid%grid%ymax
3453 ENDIF
3454
3455! reset lon in standard grib 2 definition [0,360]
3456 IF (editionnumber == 1) THEN
3457 CALL long_reset_m180_360(lofirst)
3458 CALL long_reset_m180_360(lolast)
3459 ELSE IF (editionnumber == 2) THEN
3460 CALL long_reset_0_360(lofirst)
3461 CALL long_reset_0_360(lolast)
3462 ENDIF
3463
3464 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3465 CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3466 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3467 CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3468
3469! test relative coordinate truncation error with respect to tol
3470! tol should be tuned
3471 sdx = this%grid%grid%dx*ratio
3472 sdy = this%grid%grid%dy*ratio
3473! error is computed relatively to the whole grid extension
3474 ex = abs(((sdx - nint(sdx))/ratio*this%dim%nx)/ &
3475 (this%grid%grid%xmax-this%grid%grid%xmin))
3476 ey = abs(((sdy - nint(sdy))/ratio*this%dim%ny)/ &
3477 (this%grid%grid%ymax-this%grid%grid%ymin))
3478 tol = 1.0d-3
3479 IF (ex > tol .OR. ey > tol) THEN
3480#ifdef DEBUG
3482 "griddim_export_gribapi, error on x "//&
3483 trim(to_char(ex))//"/"//t2c(tol))
3485 "griddim_export_gribapi, error on y "//&
3486 trim(to_char(ey))//"/"//t2c(tol))
3487#endif
3488! previous frmula relative to a single grid step
3489! tol = 1.0d0/ratio
3490! IF (ABS(NINT(sdx)/sdx - 1.0d0) > tol .OR. ABS(NINT(sdy)/sdy - 1.0d0) > tol) THEN
3491!#ifdef DEBUG
3492! CALL l4f_category_log(this%category,L4F_DEBUG, &
3493! "griddim_export_gribapi, dlon relative error: "//&
3494! TRIM(to_char(ABS(NINT(sdx)/sdx - 1.0d0)))//">"//TRIM(to_char(tol)))
3495! CALL l4f_category_log(this%category,L4F_DEBUG, &
3496! "griddim_export_gribapi, dlat relative error: "//&
3497! TRIM(to_char(ABS(NINT(sdy)/sdy - 1.0d0)))//">"//TRIM(to_char(tol)))
3498!#endif
3500 "griddim_export_gribapi, increments not given: inaccurate!")
3501 CALL grib_set_missing(gaid,'Di')
3502 CALL grib_set_missing(gaid,'Dj')
3503 CALL grib_set(gaid,'ijDirectionIncrementGiven',0)
3504 ELSE
3505#ifdef DEBUG
3507 "griddim_export_gribapi, setting increments: "// &
3508 trim(to_char(this%grid%grid%dx))//' '//trim(to_char(this%grid%grid%dy)))
3509#endif
3510 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
3511 CALL grib_set(gaid,'iDirectionIncrement',nint(this%grid%grid%dx*ratio))
3512 CALL grib_set(gaid,'jDirectionIncrement',nint(this%grid%grid%dy*ratio))
3513! this does not work in grib_set
3514! CALL grib_set(gaid,'iDirectionIncrementInDegrees',this%grid%grid%dx)
3515! CALL grib_set(gaid,'jDirectionIncrementInDegrees',this%grid%grid%dy)
3516 ENDIF
3517
3518! Keys for polar projections
3519CASE ('polar_stereographic', 'lambert', 'albers')
3520! increments are required
3521 CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
3522 CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
3523 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
3524! latin1/latin2 may be missing (e.g. stereographic)
3525 CALL grib_set_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
3526 CALL grib_set_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
3527! projection center flag, aka hemisphere
3528 CALL grib_set(gaid,'projectionCenterFlag',&
3529 this%grid%proj%polar%projection_center_flag, ierr)
3530 IF (ierr /= grib_success) THEN ! try center/centre
3531 CALL grib_set(gaid,'projectionCentreFlag',&
3532 this%grid%proj%polar%projection_center_flag)
3533 ENDIF
3534
3535
3536! line of view, aka central meridian
3537 CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
3538! latitude at which dx and dy are valid
3539 IF (editionnumber == 2) THEN
3540 CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3541 ENDIF
3542
3543! compute lon and lat of first point from projected extremes
3544 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
3545 lofirst, lafirst)
3546! reset lon in standard grib 2 definition [0,360]
3547 IF (editionnumber == 1) THEN
3548 CALL long_reset_m180_360(lofirst)
3549 ELSE IF (editionnumber == 2) THEN
3550 CALL long_reset_0_360(lofirst)
3551 ENDIF
3552 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3553 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3554
3555! Keys for equatorial projections
3556CASE ('mercator')
3557
3558! increments are required
3559 CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
3560 CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
3561 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
3562
3563! line of view, aka central meridian, not in grib
3564! CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
3565! latitude at which dx and dy are valid
3566 CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3567
3568! compute lon and lat of first and last points from projected extremes
3569 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
3570 lofirst, lafirst)
3571 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3572 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3573 CALL get_unproj_last(this, iscansnegatively, jscanspositively, &
3574 lolast, lalast)
3575 CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3576 CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3577
3578 IF (editionnumber == 2) THEN
3579 CALL grib_set(gaid,'orientationOfTheGridInDegrees',0.)
3580 ENDIF
3581
3582CASE ('UTM')
3583
3584 CALL grib_set(gaid,'datum',0)
3586 CALL grib_set(gaid,'referenceLongitude',nint(reflon/1.0d6))
3587 CALL grib_set(gaid,'falseEasting',this%grid%proj%xoff)
3588 CALL grib_set(gaid,'falseNorthing',this%grid%proj%yoff)
3589
3590 CALL grib_set(gaid,'iDirectionIncrement',this%grid%grid%dx)
3591 CALL grib_set(gaid,'jDirectionIncrement',this%grid%grid%dy)
3592
3593!res/scann ??
3594
3595 CALL grib_set(gaid,'zone',zone)
3596
3597 IF (iscansnegatively == 0) THEN
3598 lofirst = this%grid%grid%xmin
3599 lolast = this%grid%grid%xmax
3600 ELSE
3601 lofirst = this%grid%grid%xmax
3602 lolast = this%grid%grid%xmin
3603 ENDIF
3604 IF (jscanspositively == 0) THEN
3605 lafirst = this%grid%grid%ymax
3606 lalast = this%grid%grid%ymin
3607 ELSE
3608 lafirst = this%grid%grid%ymin
3609 lalast = this%grid%grid%ymax
3610 ENDIF
3611
3612 CALL grib_set(gaid,'eastingOfFirstGridPoint',lofirst)
3613 CALL grib_set(gaid,'eastingOfLastGridPoint',lolast)
3614 CALL grib_set(gaid,'northingOfFirstGridPoint',lafirst)
3615 CALL grib_set(gaid,'northingOfLastGridPoint',lalast)
3616
3617CASE default
3619 "griddim_export_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
3620 CALL raise_error()
3621
3622END SELECT
3623
3624! hack for position of vertical coordinate parameters
3625! buggy in grib_api
3626IF (editionnumber == 1) THEN
3627! CALL grib_get(gaid,"PVPresent",pvp) ! alias, probably useless
3628 CALL grib_get(gaid,"NV",nv)
3629#ifdef DEBUG
3631 trim(to_char(nv))//" vertical coordinate parameters")
3632#endif
3633
3634 IF (nv == 0) THEN
3635 pvl = 255
3636 ELSE
3637 SELECT CASE (this%grid%proj%proj_type)
3638 CASE ('regular_ll') ! check whether "29-32 Set to zero (reserved)" are required
3639 pvl = 33
3640 CASE ('polar_stereographic')
3641 pvl = 33
3642 CASE ('rotated_ll', 'stretched_ll', 'lambert', 'albers')
3643 pvl = 43
3644 CASE ('stretched_rotated_ll')
3645 pvl = 43
3646 CASE DEFAULT
3647 pvl = 43 !?
3648 END SELECT
3649 ENDIF
3650
3651 CALL grib_set(gaid,"pvlLocation",pvl)
3652#ifdef DEBUG
3654 trim(to_char(pvl))//" as vertical coordinate parameter location")
3655#endif
3656
3657ENDIF
3658
3659
3660CONTAINS
3661! utilities routines for grib_api, is there a better place?
3662SUBROUTINE grib_set_dmiss(gaid, key, val, default, factor)
3663INTEGER,INTENT(in) :: gaid
3664CHARACTER(len=*),INTENT(in) :: key
3665DOUBLE PRECISION,INTENT(in) :: val
3666DOUBLE PRECISION,INTENT(in),OPTIONAL :: default
3667DOUBLE PRECISION,INTENT(in),OPTIONAL :: factor
3668
3669INTEGER :: ierr
3670
3672 IF (PRESENT(factor)) THEN
3673 CALL grib_set(gaid, key, val*factor, ierr)
3674 ELSE
3675 CALL grib_set(gaid, key, val, ierr)
3676 ENDIF
3677ELSE IF (PRESENT(default)) THEN
3678 CALL grib_set(gaid, key, default, ierr)
3679ENDIF
3680
3681END SUBROUTINE grib_set_dmiss
3682
3683SUBROUTINE grib_set_imiss(gaid, key, value, default)
3684INTEGER,INTENT(in) :: gaid
3685CHARACTER(len=*),INTENT(in) :: key
3686INTEGER,INTENT(in) :: value
3687INTEGER,INTENT(in),OPTIONAL :: default
3688
3689INTEGER :: ierr
3690
3692 CALL grib_set(gaid, key, value, ierr)
3693ELSE IF (PRESENT(default)) THEN
3694 CALL grib_set(gaid, key, default, ierr)
3695ENDIF
3696
3697END SUBROUTINE grib_set_imiss
3698
3699SUBROUTINE griddim_export_ellipsoid(this, gaid)
3700TYPE(griddim_def),INTENT(in) :: this
3701INTEGER,INTENT(in) :: gaid
3702
3703INTEGER :: ellips_type, ierr
3704DOUBLE PRECISION :: r1, r2, f
3705
3707
3708IF (editionnumber == 2) THEN
3709
3710! why does it produce a message even with ierr?
3711 CALL grib_set_missing(gaid, 'scaleFactorOfRadiusOfSphericalEarth', ierr)
3712 CALL grib_set_missing(gaid, 'scaledValueOfRadiusOfSphericalEarth', ierr)
3713 CALL grib_set_missing(gaid, 'scaleFactorOfEarthMajorAxis', ierr)
3714 CALL grib_set_missing(gaid, 'scaledValueOfEarthMajorAxis', ierr)
3715 CALL grib_set_missing(gaid, 'scaleFactorOfEarthMinorAxis', ierr)
3716 CALL grib_set_missing(gaid, 'scaledValueOfEarthMinorAxis', ierr)
3717
3718 SELECT CASE(ellips_type)
3719 CASE(ellips_grs80) ! iag-grs80
3720 CALL grib_set(gaid, 'shapeOfTheEarth', 4)
3721 CASE(ellips_wgs84) ! wgs84
3722 CALL grib_set(gaid, 'shapeOfTheEarth', 5)
3723 CASE default
3724 IF (f == 0.0d0) THEN ! spherical Earth
3725 IF (r1 == 6367470.0d0) THEN ! spherical
3726 CALL grib_set(gaid, 'shapeOfTheEarth', 0)
3727 ELSE IF (r1 == 6371229.0d0) THEN ! spherical
3728 CALL grib_set(gaid, 'shapeOfTheEarth', 6)
3729 ELSE ! spherical generic
3730 CALL grib_set(gaid, 'shapeOfTheEarth', 1)
3731 CALL grib_set(gaid, 'scaleFactorOfRadiusOfSphericalEarth', 2)
3732 CALL grib_set(gaid, 'scaledValueOfRadiusOfSphericalEarth', int(r1*100.0d0))
3733 ENDIF
3734 ELSE ! ellipsoidal
3735 IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
3736 CALL grib_set(gaid, 'shapeOfTheEarth', 2)
3737 ELSE ! ellipsoidal generic
3738 CALL grib_set(gaid, 'shapeOfTheEarth', 3)
3739 r2 = r1*(1.0d0 - f)
3740 CALL grib_set(gaid, 'scaleFactorOfEarthMajorAxis', 5)
3741 CALL grib_set(gaid, 'scaledValueOfEarthMajorAxis', &
3742 int(r1*100.0d0))
3743 CALL grib_set(gaid, 'scaleFactorOfEarthMinorAxis', 5)
3744 CALL grib_set(gaid, 'scaledValueOfEarthMinorAxis', &
3745 int(r2*100.0d0))
3746 ENDIF
3747 ENDIF
3748 END SELECT
3749
3750ELSE
3751
3752 IF (r1 == 6367470.0d0 .AND. f == 0.0d0) THEN ! spherical
3753 CALL grib_set(gaid, 'earthIsOblate', 0)
3754 ELSE IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
3755 CALL grib_set(gaid, 'earthIsOblate', 1)
3756 ELSE IF (f == 0.0d0) THEN ! generic spherical
3757 CALL grib_set(gaid, 'earthIsOblate', 0)
3759 ¬ supported in grib 1, coding with default radius of 6367470 m')
3760 ELSE ! generic ellipsoidal
3761 CALL grib_set(gaid, 'earthIsOblate', 1)
3763 ¬ supported in grib 1, coding with default iau65 ellipsoid')
3764 ENDIF
3765
3766ENDIF
3767
3768END SUBROUTINE griddim_export_ellipsoid
3769
3770SUBROUTINE get_unproj_first(this, iScansNegatively, jScansPositively, &
3771 loFirst, laFirst)
3772TYPE(griddim_def),INTENT(in) :: this ! griddim object
3773INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
3774DOUBLE PRECISION,INTENT(out) :: loFirst, laFirst
3775
3776! compute lon and lat of first point from projected extremes
3777IF (iscansnegatively == 0) THEN
3778 IF (jscanspositively == 0) THEN
3780 ELSE
3782 ENDIF
3783ELSE
3784 IF (jscanspositively == 0) THEN
3786 ELSE
3788 ENDIF
3789ENDIF
3790! use the values kept for personal pleasure ?
3791! loFirst = this%grid%proj%polar%lon1
3792! laFirst = this%grid%proj%polar%lat1
3793END SUBROUTINE get_unproj_first
3794
3795SUBROUTINE get_unproj_last(this, iScansNegatively, jScansPositively, &
3796 loLast, laLast)
3797TYPE(griddim_def),INTENT(in) :: this ! griddim object
3798INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
3799DOUBLE PRECISION,INTENT(out) :: loLast, laLast
3800
3801! compute lon and lat of last point from projected extremes
3802IF (iscansnegatively == 0) THEN
3803 IF (jscanspositively == 0) THEN
3805 ELSE
3807 ENDIF
3808ELSE
3809 IF (jscanspositively == 0) THEN
3811 ELSE
3813 ENDIF
3814ENDIF
3815! use the values kept for personal pleasure ?
3816! loLast = this%grid%proj%polar%lon?
3817! laLast = this%grid%proj%polar%lat?
3818END SUBROUTINE get_unproj_last
3819
3820END SUBROUTINE griddim_export_gribapi
3821#endif
3822
3823
3824#ifdef HAVE_LIBGDAL
3825! gdal driver
3826SUBROUTINE griddim_import_gdal(this, gdalid, gdal_options)
3827USE gdal
3828TYPE(griddim_def),INTENT(inout) :: this ! griddim object
3829TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal rasterband pointer
3830TYPE(gdal_file_id_options),INTENT(in) :: gdal_options
3831
3832TYPE(gdaldataseth) :: hds
3833REAL(kind=c_double) :: geotrans(6) !, invgeotrans(6), x1, y1, x2, y2
3834INTEGER(kind=c_int) :: offsetx, offsety
3835INTEGER :: ier
3836
3837hds = gdalgetbanddataset(gdalid) ! go back to dataset
3838ier = gdalgetgeotransform(hds, geotrans)
3839
3840IF (ier /= 0) THEN
3842 'griddim_import_gdal: error in accessing gdal dataset')
3843 CALL raise_error()
3844 RETURN
3845ENDIF
3846IF (geotrans(3) /= 0.0_c_double .OR. geotrans(5) /= 0.0_c_double) THEN ! dataset has not a diagonal transformation
3848 'griddim_import_gdal: dataset has a non-diagonal transformation matrix, unsupported')
3849 CALL raise_error()
3850 RETURN
3851ENDIF
3852
3853CALL gdaldatasetbbsize_f(hds, gdal_options%xmin, gdal_options%ymin, &
3854 gdal_options%xmax, gdal_options%ymax, &
3855 this%dim%nx, this%dim%ny, offsetx, offsety, &
3856 this%grid%grid%xmin, this%grid%grid%ymin, this%grid%grid%xmax, this%grid%grid%ymax)
3857
3858IF (this%dim%nx == 0 .OR. this%dim%ny == 0) THEN
3860 'griddim_import_gdal: requested bounding box '//t2c(gdal_options%xmin)//','// &
3861 t2c(gdal_options%ymin)//','//t2c(gdal_options%xmax)//','//&
3862 t2c(gdal_options%ymax))
3864 'determines an empty dataset '//t2c(this%dim%nx)//'x'//t2c(this%dim%ny))
3865ENDIF
3866
3867! get grid corners
3868!CALL gdalapplygeotransform(geotrans, 0.5_c_double, 0.5_c_double, x1, y1)
3869!CALL gdalapplygeotransform(geotrans, gdalgetrasterbandxsize(gdalid)-0.5_c_double, &
3870! gdalgetrasterbandysize(gdalid)-0.5_c_double, x2, y2)
3871
3872!IF (geotrans(3) == 0.0_c_double .AND. geotrans(5) == 0.0_c_double) THEN ! transformation is diagonal, no transposing
3873! this%dim%nx = gdalgetrasterbandxsize(gdalid)
3874! this%dim%ny = gdalgetrasterbandysize(gdalid)
3875! this%grid%grid%xmin = MIN(x1, x2)
3876! this%grid%grid%xmax = MAX(x1, x2)
3877! this%grid%grid%ymin = MIN(y1, y2)
3878! this%grid%grid%ymax = MAX(y1, y2)
3879!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
3880!
3881! this%dim%nx = gdalgetrasterbandysize(gdalid)
3882! this%dim%ny = gdalgetrasterbandxsize(gdalid)
3883! this%grid%grid%xmin = MIN(y1, y2)
3884! this%grid%grid%xmax = MAX(y1, y2)
3885! this%grid%grid%ymin = MIN(x1, x2)
3886! this%grid%grid%ymax = MAX(x1, x2)
3887!
3888!ELSE ! transformation is a rotation, not supported
3889!ENDIF
3890
3891this%grid%proj%proj_type = 'regular_ll' ! forced, only one supported (check?)
3892
3893CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3894this%grid%grid%component_flag = 0
3895
3896END SUBROUTINE griddim_import_gdal
3897#endif
3898
3899
3901SUBROUTINE griddim_display(this)
3902TYPE(griddim_def),INTENT(in) :: this
3903
3904print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
3905
3909
3910print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
3911
3912END SUBROUTINE griddim_display
3913
3914
3915#define VOL7D_POLY_TYPE TYPE(grid_def)
3916#define VOL7D_POLY_TYPES _grid
3917#include "array_utilities_inc.F90"
3918#undef VOL7D_POLY_TYPE
3919#undef VOL7D_POLY_TYPES
3920
3921#define VOL7D_POLY_TYPE TYPE(griddim_def)
3922#define VOL7D_POLY_TYPES _griddim
3923#include "array_utilities_inc.F90"
3924#undef VOL7D_POLY_TYPE
3925#undef VOL7D_POLY_TYPES
3926
3927
3939SUBROUTINE griddim_wind_unrot(this, rot_mat)
3941TYPE(griddim_def), INTENT(in) :: this
3942DOUBLE PRECISION, POINTER :: rot_mat(:,:,:)
3943
3944DOUBLE PRECISION :: dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j
3945DOUBLE PRECISION :: lat_factor
3946INTEGER :: i, j, ip1, im1, jp1, jm1;
3947
3948IF (this%dim%nx <= 0 .OR. this%dim%ny <= 0 .OR. &
3949 .NOT. ASSOCIATED(this%dim%lon) .OR. .NOT. ASSOCIATED(this%dim%lat)) THEN
3950 CALL l4f_category_log(this%category, l4f_error, 'rotate_uv must be called after correct init of griddim object')
3951 NULLIFY(rot_mat)
3952 RETURN
3953ENDIF
3954
3955ALLOCATE(rot_mat(this%dim%nx, this%dim%ny, 4))
3956
3957DO j = 1, this%dim%ny
3958 jp1 = min(j+1, this%dim%ny)
3959 jm1 = max(j-1, 1)
3960 DO i = 1, this%dim%nx
3961 ip1 = min(i+1, this%dim%nx)
3962 im1 = max(i-1, 1)
3963
3964 dlat_i = this%dim%lat(ip1,j) - this%dim%lat(im1,j)
3965 dlat_j = this%dim%lat(i,jp1) - this%dim%lat(i,jm1)
3966
3967 dlon_i = this%dim%lon(ip1,j) - this%dim%lon(im1,j)
3968! if ( dlon_i > pi ) dlon_i -= 2*pi;
3969! if ( dlon_i < (-pi) ) dlon_i += 2*pi;
3970 dlon_j = this%dim%lon(i,jp1) - this%dim%lon(i,jm1)
3971! if ( dlon_j > pi ) dlon_j -= 2*pi;
3972! if ( dlon_j < (-pi) ) dlon_j += 2*pi;
3973
3974! check whether this is really necessary !!!!
3975 lat_factor = cos(degrad*this%dim%lat(i,j))
3976 dlon_i = dlon_i * lat_factor
3977 dlon_j = dlon_j * lat_factor
3978
3979 dist_i = sqrt(dlon_i*dlon_i + dlat_i*dlat_i)
3980 dist_j = sqrt(dlon_j*dlon_j + dlat_j*dlat_j)
3981
3982 IF (dist_i > 0.d0) THEN
3983 rot_mat(i,j,1) = dlon_i/dist_i
3984 rot_mat(i,j,2) = dlat_i/dist_i
3985 ELSE
3986 rot_mat(i,j,1) = 0.0d0
3987 rot_mat(i,j,2) = 0.0d0
3988 ENDIF
3989 IF (dist_j > 0.d0) THEN
3990 rot_mat(i,j,3) = dlon_j/dist_j
3991 rot_mat(i,j,4) = dlat_j/dist_j
3992 ELSE
3993 rot_mat(i,j,3) = 0.0d0
3994 rot_mat(i,j,4) = 0.0d0
3995 ENDIF
3996
3997 ENDDO
3998ENDDO
3999
4000END SUBROUTINE griddim_wind_unrot
4001
4002
4003! compute zoom indices from geographical zoom coordinates
4004SUBROUTINE griddim_zoom_coord(this, ilon, ilat, flon, flat, ix, iy, fx, fy)
4005TYPE(griddim_def),INTENT(in) :: this
4006DOUBLE PRECISION,INTENT(in) :: ilon, ilat, flon, flat
4007INTEGER,INTENT(out) :: ix, iy, fx, fy
4008
4009DOUBLE PRECISION :: ix1, iy1, fx1, fy1
4010
4011! compute projected coordinates of vertices of desired lonlat rectangle
4014
4015CALL griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
4016
4017END SUBROUTINE griddim_zoom_coord
4018
4019
4020! compute zoom indices from projected zoom coordinates
4021SUBROUTINE griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
4022TYPE(griddim_def),INTENT(in) :: this
4023DOUBLE PRECISION,INTENT(in) :: ix1, iy1, fx1, fy1
4024INTEGER,INTENT(out) :: ix, iy, fx, fy
4025
4026INTEGER :: lix, liy, lfx, lfy
4027
4028! compute projected indices
4029lix = nint((ix1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
4030liy = nint((iy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
4031lfx = nint((fx1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
4032lfy = nint((fy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
4033! swap projected indices, in case grid is "strongly rotated"
4034ix = min(lix, lfx)
4035fx = max(lix, lfx)
4036iy = min(liy, lfy)
4037fy = max(liy, lfy)
4038
4039END SUBROUTINE griddim_zoom_projcoord
4040
4041
4045SUBROUTINE long_reset_0_360(lon)
4046DOUBLE PRECISION,INTENT(inout) :: lon
4047
4049DO WHILE(lon < 0.0d0)
4050 lon = lon + 360.0d0
4051END DO
4052DO WHILE(lon >= 360.0d0)
4053 lon = lon - 360.0d0
4054END DO
4055
4056END SUBROUTINE long_reset_0_360
4057
4058
4062SUBROUTINE long_reset_m180_360(lon)
4063DOUBLE PRECISION,INTENT(inout) :: lon
4064
4066DO WHILE(lon < -180.0d0)
4067 lon = lon + 360.0d0
4068END DO
4069DO WHILE(lon >= 360.0d0)
4070 lon = lon - 360.0d0
4071END DO
4072
4073END SUBROUTINE long_reset_m180_360
4074
4075
4079!SUBROUTINE long_reset_m90_270(lon)
4080!DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
4081!
4082!IF (.NOT.c_e(lon)) RETURN
4083!DO WHILE(lon < -90.0D0)
4084! lon = lon + 360.0D0
4085!END DO
4086!DO WHILE(lon >= 270.0D0)
4087! lon = lon - 360.0D0
4088!END DO
4089!
4090!END SUBROUTINE long_reset_m90_270
4091
4092
4096SUBROUTINE long_reset_m180_180(lon)
4097DOUBLE PRECISION,INTENT(inout) :: lon
4098
4100DO WHILE(lon < -180.0d0)
4101 lon = lon + 360.0d0
4102END DO
4103DO WHILE(lon >= 180.0d0)
4104 lon = lon - 360.0d0
4105END DO
4106
4107END SUBROUTINE long_reset_m180_180
4108
4109
4110SUBROUTINE long_reset_to_cart_closest(lon, lonref)
4111DOUBLE PRECISION,INTENT(inout) :: lon
4112DOUBLE PRECISION,INTENT(in) :: lonref
4113
4115IF (abs(lon-lonref) < 180.0d0) RETURN ! nothing to do
4116lon = lon - nint((lon-lonref)/360.0d0)*360.0d0 ! se non e` vera e` ben trovata
4117
4118END SUBROUTINE long_reset_to_cart_closest
4119
4120
4122
Method for testing the existence of the object. Definition: geo_proj_class.F90:274 Compute forward coordinate transformation from geographical system to projected system. Definition: grid_class.F90:314 Read the object from a formatted or unformatted file. Definition: grid_class.F90:340 Compute backward coordinate transformation from projected system to geographical system. Definition: grid_class.F90:320 Write the object on a formatted or unformatted file. Definition: grid_class.F90:335 Emit log message for a category with specific priority. Definition: log4fortran.F90:463 Module for describing geographically referenced regular grids. Definition: grid_class.F90:243 Module for defining the extension and coordinates of a rectangular georeferenced grid. Definition: grid_dim_class.F90:217 This module defines an abstract interface to different drivers for access to files containing gridded... Definition: grid_id_class.F90:255 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:263 This object completely describes a grid on a geographic projection. Definition: grid_class.F90:276 Derived type describing the extension of a grid and the geographical coordinates of each point. Definition: grid_dim_class.F90:227 |