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