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