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