libsim Versione 7.2.1
|
◆ pack_distinct_grid()
compatta gli elementi distinti di vect in un array Definizione alla linea 2058 del file grid_class.F90. 2060! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2061! authors:
2062! Davide Cesari <dcesari@arpa.emr.it>
2063! Paolo Patruno <ppatruno@arpa.emr.it>
2064
2065! This program is free software; you can redistribute it and/or
2066! modify it under the terms of the GNU General Public License as
2067! published by the Free Software Foundation; either version 2 of
2068! the License, or (at your option) any later version.
2069
2070! This program is distributed in the hope that it will be useful,
2071! but WITHOUT ANY WARRANTY; without even the implied warranty of
2072! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2073! GNU General Public License for more details.
2074
2075! You should have received a copy of the GNU General Public License
2076! along with this program. If not, see <http://www.gnu.org/licenses/>.
2077#include "config.h"
2078
2109use geo_proj_class
2111use grid_rect_class
2117implicit none
2118
2119
2120character (len=255),parameter:: subcategory="grid_class"
2121
2122
2129 private
2130 type(geo_proj) :: proj
2131 type(grid_rect) :: grid
2132 integer :: category = 0
2134
2135
2142 type(grid_def) :: grid
2143 type(grid_dim) :: dim
2144 integer :: category = 0
2146
2147
2151INTERFACE OPERATOR (==)
2152 MODULE PROCEDURE grid_eq, griddim_eq
2153END INTERFACE
2154
2158INTERFACE OPERATOR (/=)
2159 MODULE PROCEDURE grid_ne, griddim_ne
2160END INTERFACE
2161
2164 MODULE PROCEDURE griddim_init
2165END INTERFACE
2166
2169 MODULE PROCEDURE griddim_delete
2170END INTERFACE
2171
2174 MODULE PROCEDURE griddim_copy
2175END INTERFACE
2176
2180 MODULE PROCEDURE griddim_coord_proj!, griddim_proj
2181END INTERFACE
2182
2186 MODULE PROCEDURE griddim_coord_unproj, griddim_unproj
2187END INTERFACE
2188
2191 MODULE PROCEDURE griddim_get_val
2192END INTERFACE
2193
2196 MODULE PROCEDURE griddim_set_val
2197END INTERFACE
2198
2201 MODULE PROCEDURE griddim_write_unit
2202END INTERFACE
2203
2206 MODULE PROCEDURE griddim_read_unit
2207END INTERFACE
2208
2210INTERFACE import
2211 MODULE PROCEDURE griddim_import_grid_id
2212END INTERFACE
2213
2216 MODULE PROCEDURE griddim_export_grid_id
2217END INTERFACE
2218
2221 MODULE PROCEDURE griddim_display
2222END INTERFACE
2223
2224#define VOL7D_POLY_TYPE TYPE(grid_def)
2225#define VOL7D_POLY_TYPES _grid
2226#include "array_utilities_pre.F90"
2227#undef VOL7D_POLY_TYPE
2228#undef VOL7D_POLY_TYPES
2229
2230#define VOL7D_POLY_TYPE TYPE(griddim_def)
2231#define VOL7D_POLY_TYPES _griddim
2232#include "array_utilities_pre.F90"
2233#undef VOL7D_POLY_TYPE
2234#undef VOL7D_POLY_TYPES
2235
2236INTERFACE wind_unrot
2237 MODULE PROCEDURE griddim_wind_unrot
2238END INTERFACE
2239
2240
2241PRIVATE
2242
2244 griddim_zoom_coord, griddim_zoom_projcoord, &
2248PUBLIC OPERATOR(==),OPERATOR(/=)
2249PUBLIC count_distinct, pack_distinct, count_and_pack_distinct, &
2250 map_distinct, map_inv_distinct,index
2252PUBLIC griddim_central_lon, griddim_set_central_lon
2253CONTAINS
2254
2256SUBROUTINE griddim_init(this, nx, ny, &
2257 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
2258 proj_type, lov, zone, xoff, yoff, &
2259 longitude_south_pole, latitude_south_pole, angle_rotation, &
2260 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
2261 latin1, latin2, lad, projection_center_flag, &
2262 ellips_smaj_axis, ellips_flatt, ellips_type, &
2263 categoryappend)
2264TYPE(griddim_def),INTENT(inout) :: this
2265INTEGER,INTENT(in),OPTIONAL :: nx
2266INTEGER,INTENT(in),OPTIONAL :: ny
2267DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin
2268DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmax
2269DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymin
2270DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymax
2271DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx
2272DOUBLE PRECISION,INTENT(in),OPTIONAL :: dy
2275INTEGER,INTENT(in),OPTIONAL :: component_flag
2276CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
2277DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
2278INTEGER,INTENT(in),OPTIONAL :: zone
2279DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
2280DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
2281DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
2282DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
2283DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
2284DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
2285DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
2286DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
2287DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
2288DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
2289DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
2290INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
2291DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
2292DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
2293INTEGER,INTENT(in),OPTIONAL :: ellips_type
2294CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2295
2296CHARACTER(len=512) :: a_name
2297
2298IF (PRESENT(categoryappend)) THEN
2299 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
2300 trim(categoryappend))
2301ELSE
2302 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2303ENDIF
2304this%category=l4f_category_get(a_name)
2305
2306! geographical projection
2307this%grid%proj = geo_proj_new( &
2308 proj_type=proj_type, lov=lov, zone=zone, xoff=xoff, yoff=yoff, &
2309 longitude_south_pole=longitude_south_pole, &
2310 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
2311 longitude_stretch_pole=longitude_stretch_pole, &
2312 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
2313 latin1=latin1, latin2=latin2, lad=lad, projection_center_flag=projection_center_flag, &
2314 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, ellips_type=ellips_type)
2315! grid extension
2316this%grid%grid = grid_rect_new( &
2317 xmin, xmax, ymin, ymax, dx, dy, component_flag)
2318! grid size
2319this%dim = grid_dim_new(nx, ny)
2320
2321#ifdef DEBUG
2323#endif
2324
2325END SUBROUTINE griddim_init
2326
2327
2329SUBROUTINE griddim_delete(this)
2330TYPE(griddim_def),INTENT(inout) :: this
2331
2335
2336call l4f_category_delete(this%category)
2337
2338END SUBROUTINE griddim_delete
2339
2340
2342SUBROUTINE griddim_copy(this, that, categoryappend)
2343TYPE(griddim_def),INTENT(in) :: this
2344TYPE(griddim_def),INTENT(out) :: that
2345CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2346
2347CHARACTER(len=512) :: a_name
2348
2350
2354
2355! new category
2356IF (PRESENT(categoryappend)) THEN
2357 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
2358 trim(categoryappend))
2359ELSE
2360 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2361ENDIF
2362that%category=l4f_category_get(a_name)
2363
2364END SUBROUTINE griddim_copy
2365
2366
2369ELEMENTAL SUBROUTINE griddim_coord_proj(this, lon, lat, x, y)
2370TYPE(griddim_def),INTENT(in) :: this
2372DOUBLE PRECISION,INTENT(in) :: lon, lat
2374DOUBLE PRECISION,INTENT(out) :: x, y
2375
2377
2378END SUBROUTINE griddim_coord_proj
2379
2380
2383ELEMENTAL SUBROUTINE griddim_coord_unproj(this, x, y, lon, lat)
2384TYPE(griddim_def),INTENT(in) :: this
2386DOUBLE PRECISION,INTENT(in) :: x, y
2388DOUBLE PRECISION,INTENT(out) :: lon, lat
2389
2391
2392END SUBROUTINE griddim_coord_unproj
2393
2394
2395! Computes and sets the grid parameters required to compute
2396! coordinates of grid points in the projected system.
2397! probably meaningless
2398!SUBROUTINE griddim_proj(this)
2399!TYPE(griddim_def),INTENT(inout) :: this !< definition of projection and grid
2400!
2401!CALL proj(this, this%dim%lon(1,1), this%dim%lat(1,1), &
2402! this%grid%grid%xmin, this%grid%grid%ymin)
2403!
2404!CALL proj(this, this%dim%lon(this%dim%nx,this%dim%ny), &
2405! this%dim%lat(this%dim%nx,this%dim%ny), &
2406! this%grid%grid%xmax, this%grid%grid%ymax)
2407!
2408!END SUBROUTINE griddim_proj
2409
2417SUBROUTINE griddim_unproj(this)
2418TYPE(griddim_def),INTENT(inout) :: this
2419
2421CALL alloc(this%dim)
2422CALL griddim_unproj_internal(this)
2423
2424END SUBROUTINE griddim_unproj
2425
2426! internal subroutine needed for allocating automatic arrays
2427SUBROUTINE griddim_unproj_internal(this)
2428TYPE(griddim_def),INTENT(inout) ::this ! definition of projection and grid
2429
2430DOUBLE PRECISION :: x(this%dim%nx,this%dim%ny), y(this%dim%nx,this%dim%ny)
2431
2432CALL grid_rect_coordinates(this%grid%grid, x, y)
2434
2435END SUBROUTINE griddim_unproj_internal
2436
2437
2439SUBROUTINE griddim_get_val(this, nx, ny, &
2440 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
2441 proj, proj_type, lov, zone, xoff, yoff, &
2442 longitude_south_pole, latitude_south_pole, angle_rotation, &
2443 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
2444 latin1, latin2, lad, projection_center_flag, &
2445 ellips_smaj_axis, ellips_flatt, ellips_type)
2446TYPE(griddim_def),INTENT(in) :: this
2447INTEGER,INTENT(out),OPTIONAL :: nx
2448INTEGER,INTENT(out),OPTIONAL :: ny
2450DOUBLE PRECISION,INTENT(out),OPTIONAL :: xmin, xmax, ymin, ymax
2452DOUBLE PRECISION,INTENT(out),OPTIONAL :: dx, dy
2455INTEGER,INTENT(out),OPTIONAL :: component_flag
2456TYPE(geo_proj),INTENT(out),OPTIONAL :: proj
2457CHARACTER(len=*),INTENT(out),OPTIONAL :: proj_type
2458DOUBLE PRECISION,INTENT(out),OPTIONAL :: lov
2459INTEGER,INTENT(out),OPTIONAL :: zone
2460DOUBLE PRECISION,INTENT(out),OPTIONAL :: xoff
2461DOUBLE PRECISION,INTENT(out),OPTIONAL :: yoff
2462DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_south_pole
2463DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_south_pole
2464DOUBLE PRECISION,INTENT(out),OPTIONAL :: angle_rotation
2465DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_stretch_pole
2466DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_stretch_pole
2467DOUBLE PRECISION,INTENT(out),OPTIONAL :: stretch_factor
2468DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin1
2469DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin2
2470DOUBLE PRECISION,INTENT(out),OPTIONAL :: lad
2471INTEGER,INTENT(out),OPTIONAL :: projection_center_flag
2472DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_smaj_axis
2473DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_flatt
2474INTEGER,INTENT(out),OPTIONAL :: ellips_type
2475
2476IF (PRESENT(nx)) nx = this%dim%nx
2477IF (PRESENT(ny)) ny = this%dim%ny
2478
2480
2482 xoff=xoff, yoff=yoff, &
2483 longitude_south_pole=longitude_south_pole, &
2484 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
2485 longitude_stretch_pole=longitude_stretch_pole, &
2486 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
2487 latin1=latin1, latin2=latin2, lad=lad, &
2488 projection_center_flag=projection_center_flag, &
2489 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
2490 ellips_type=ellips_type)
2491
2493 xmin, xmax, ymin, ymax, dx, dy, component_flag)
2494
2495END SUBROUTINE griddim_get_val
2496
2497
2499SUBROUTINE griddim_set_val(this, nx, ny, &
2500 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
2501 proj_type, lov, zone, xoff, yoff, &
2502 longitude_south_pole, latitude_south_pole, angle_rotation, &
2503 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
2504 latin1, latin2, lad, projection_center_flag, &
2505 ellips_smaj_axis, ellips_flatt, ellips_type)
2506TYPE(griddim_def),INTENT(inout) :: this
2507INTEGER,INTENT(in),OPTIONAL :: nx
2508INTEGER,INTENT(in),OPTIONAL :: ny
2510DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin, xmax, ymin, ymax
2512DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx, dy
2515INTEGER,INTENT(in),OPTIONAL :: component_flag
2516CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
2517DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
2518INTEGER,INTENT(in),OPTIONAL :: zone
2519DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
2520DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
2521DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
2522DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
2523DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
2524DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
2525DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
2526DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
2527DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
2528DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
2529DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
2530INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
2531DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
2532DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
2533INTEGER,INTENT(in),OPTIONAL :: ellips_type
2534
2535IF (PRESENT(nx)) this%dim%nx = nx
2536IF (PRESENT(ny)) this%dim%ny = ny
2537
2539 xoff=xoff, yoff=yoff, longitude_south_pole=longitude_south_pole, &
2540 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
2541 longitude_stretch_pole=longitude_stretch_pole, &
2542 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
2543 latin1=latin1, latin2=latin2, lad=lad, &
2544 projection_center_flag=projection_center_flag, &
2545 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
2546 ellips_type=ellips_type)
2547
2549 xmin, xmax, ymin, ymax, dx, dy, component_flag)
2550
2551END SUBROUTINE griddim_set_val
2552
2553
2558SUBROUTINE griddim_read_unit(this, unit)
2559TYPE(griddim_def),INTENT(out) :: this
2560INTEGER, INTENT(in) :: unit
2561
2562
2566
2567END SUBROUTINE griddim_read_unit
2568
2569
2574SUBROUTINE griddim_write_unit(this, unit)
2575TYPE(griddim_def),INTENT(in) :: this
2576INTEGER, INTENT(in) :: unit
2577
2578
2582
2583END SUBROUTINE griddim_write_unit
2584
2585
2589FUNCTION griddim_central_lon(this) RESULT(lon)
2590TYPE(griddim_def),INTENT(inout) :: this
2591
2592DOUBLE PRECISION :: lon
2593
2594CALL griddim_pistola_central_lon(this, lon)
2595
2596END FUNCTION griddim_central_lon
2597
2598
2602SUBROUTINE griddim_set_central_lon(this, lonref)
2603TYPE(griddim_def),INTENT(inout) :: this
2604DOUBLE PRECISION,INTENT(in) :: lonref
2605
2606DOUBLE PRECISION :: lon
2607
2608CALL griddim_pistola_central_lon(this, lon, lonref)
2609
2610END SUBROUTINE griddim_set_central_lon
2611
2612
2613! internal subroutine for performing tasks common to the prevous two
2614SUBROUTINE griddim_pistola_central_lon(this, lon, lonref)
2615TYPE(griddim_def),INTENT(inout) :: this ! grid descriptor
2616DOUBLE PRECISION,INTENT(inout) :: lon ! central longitude
2617DOUBLE PRECISION,INTENT(in),OPTIONAL :: lonref ! reference longitude
2618
2619INTEGER :: unit
2620DOUBLE PRECISION :: lonsp, latsp, londelta, lov, lonrot
2621CHARACTER(len=80) :: ptype
2622
2623lon = dmiss
2625IF (unit == geo_proj_unit_meter) THEN ! it is a plane projection
2627 IF (PRESENT(lonref)) THEN
2628 CALL long_reset_to_cart_closest(lov, lonref)
2630 ENDIF
2631
2632ELSE IF (unit == geo_proj_unit_degree) THEN ! it is a spheric projection
2634 longitude_south_pole=lonsp, latitude_south_pole=latsp)
2635 SELECT CASE(ptype)
2636 CASE('rotated_ll','stretched_rotated_ll') ! use origin of rotated system
2637 IF (latsp < 0.0d0) THEN
2638 lon = lonsp
2639 IF (PRESENT(lonref)) THEN
2640 CALL long_reset_to_cart_closest(lon, lonref)
2642! now reset rotated coordinates around zero
2644 lonrot = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
2645 ENDIF
2646 londelta = lonrot
2647 CALL long_reset_to_cart_closest(londelta, 0.0d0)
2648 londelta = londelta - lonrot
2649 this%grid%grid%xmin = this%grid%grid%xmin + londelta
2650 this%grid%grid%xmax = this%grid%grid%xmax + londelta
2651 ENDIF
2652 ELSE ! this part to be checked
2653 lon = modulo(lonsp + 180.0d0, 360.0d0)
2654! IF (PRESENT(lonref)) THEN
2655! CALL long_reset_to_cart_closest(lov, lonref)
2656! CALL set_val(this%grid%proj, longitude_south_pole=lonref)
2657! ENDIF
2658 ENDIF
2659 CASE default ! use real grid limits
2661 lon = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
2662 ENDIF
2663 IF (PRESENT(lonref)) THEN
2664 londelta = lon
2665 CALL long_reset_to_cart_closest(londelta, lonref)
2666 londelta = londelta - lon
2667 this%grid%grid%xmin = this%grid%grid%xmin + londelta
2668 this%grid%grid%xmax = this%grid%grid%xmax + londelta
2669 ENDIF
2670 END SELECT
2671ENDIF
2672
2673END SUBROUTINE griddim_pistola_central_lon
2674
2675
2679SUBROUTINE griddim_gen_coord(this, x, y)
2680TYPE(griddim_def),INTENT(in) :: this
2681DOUBLE PRECISION,INTENT(out) :: x(:,:)
2682DOUBLE PRECISION,INTENT(out) :: y(:,:)
2683
2684
2685CALL grid_rect_coordinates(this%grid%grid, x, y)
2686
2687END SUBROUTINE griddim_gen_coord
2688
2689
2691SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
2692TYPE(griddim_def), INTENT(in) :: this
2693INTEGER,INTENT(in) :: nx
2694INTEGER,INTENT(in) :: ny
2695DOUBLE PRECISION,INTENT(out) :: dx
2696DOUBLE PRECISION,INTENT(out) :: dy
2697
2698CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
2699
2700END SUBROUTINE griddim_steps
2701
2702
2704SUBROUTINE griddim_setsteps(this)
2705TYPE(griddim_def), INTENT(inout) :: this
2706
2707CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
2708
2709END SUBROUTINE griddim_setsteps
2710
2711
2712! TODO
2713! bisogna sviluppare gli altri operatori
2714ELEMENTAL FUNCTION grid_eq(this, that) RESULT(res)
2715TYPE(grid_def),INTENT(IN) :: this, that
2716LOGICAL :: res
2717
2718res = this%proj == that%proj .AND. &
2719 this%grid == that%grid
2720
2721END FUNCTION grid_eq
2722
2723
2724ELEMENTAL FUNCTION griddim_eq(this, that) RESULT(res)
2725TYPE(griddim_def),INTENT(IN) :: this, that
2726LOGICAL :: res
2727
2728res = this%grid == that%grid .AND. &
2729 this%dim == that%dim
2730
2731END FUNCTION griddim_eq
2732
2733
2734ELEMENTAL FUNCTION grid_ne(this, that) RESULT(res)
2735TYPE(grid_def),INTENT(IN) :: this, that
2736LOGICAL :: res
2737
2738res = .NOT.(this == that)
2739
2740END FUNCTION grid_ne
2741
2742
2743ELEMENTAL FUNCTION griddim_ne(this, that) RESULT(res)
2744TYPE(griddim_def),INTENT(IN) :: this, that
2745LOGICAL :: res
2746
2747res = .NOT.(this == that)
2748
2749END FUNCTION griddim_ne
2750
2751
2757SUBROUTINE griddim_import_grid_id(this, ingrid_id)
2758#ifdef HAVE_LIBGDAL
2759USE gdal
2760#endif
2761TYPE(griddim_def),INTENT(inout) :: this
2762TYPE(grid_id),INTENT(in) :: ingrid_id
2763
2764#ifdef HAVE_LIBGRIBAPI
2765INTEGER :: gaid
2766#endif
2767#ifdef HAVE_LIBGDAL
2768TYPE(gdalrasterbandh) :: gdalid
2769#endif
2771
2772#ifdef HAVE_LIBGRIBAPI
2773gaid = grid_id_get_gaid(ingrid_id)
2775#endif
2776#ifdef HAVE_LIBGDAL
2777gdalid = grid_id_get_gdalid(ingrid_id)
2778IF (gdalassociated(gdalid)) CALL griddim_import_gdal(this, gdalid, &
2779 grid_id_get_gdal_options(ingrid_id))
2780#endif
2781
2782END SUBROUTINE griddim_import_grid_id
2783
2784
2789SUBROUTINE griddim_export_grid_id(this, outgrid_id)
2790#ifdef HAVE_LIBGDAL
2791USE gdal
2792#endif
2793TYPE(griddim_def),INTENT(in) :: this
2794TYPE(grid_id),INTENT(inout) :: outgrid_id
2795
2796#ifdef HAVE_LIBGRIBAPI
2797INTEGER :: gaid
2798#endif
2799#ifdef HAVE_LIBGDAL
2800TYPE(gdalrasterbandh) :: gdalid
2801#endif
2802
2803#ifdef HAVE_LIBGRIBAPI
2804gaid = grid_id_get_gaid(outgrid_id)
2806#endif
2807#ifdef HAVE_LIBGDAL
2808gdalid = grid_id_get_gdalid(outgrid_id)
2809!IF (gdalassociated(gdalid)
2810! export for gdal not implemented, log?
2811#endif
2812
2813END SUBROUTINE griddim_export_grid_id
2814
2815
2816#ifdef HAVE_LIBGRIBAPI
2817! grib_api driver
2818SUBROUTINE griddim_import_gribapi(this, gaid)
2819USE grib_api
2820TYPE(griddim_def),INTENT(inout) :: this ! griddim object
2821INTEGER, INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
2822
2823DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, x1, y1, orient
2824INTEGER :: EditionNumber, iScansNegatively, jScansPositively, zone, datum, &
2825 reflon, ierr
2826
2827! Generic keys
2828CALL grib_get(gaid, 'typeOfGrid', this%grid%proj%proj_type)
2829#ifdef DEBUG
2831 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type))
2832#endif
2833CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
2834
2835IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
2836 CALL grib_get(gaid,'numberOfDataPoints', this%dim%nx)
2837 this%dim%ny = 1
2838 this%grid%grid%component_flag = 0
2839 CALL griddim_import_ellipsoid(this, gaid)
2840 RETURN
2841ENDIF
2842
2843! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
2844CALL grib_get(gaid, 'Ni', this%dim%nx)
2845CALL grib_get(gaid, 'Nj', this%dim%ny)
2846CALL griddim_import_ellipsoid(this, gaid) ! placed here, not valid for utm datum /= 1
2847
2848CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
2849CALL grib_get(gaid,'jScansPositively',jscanspositively)
2850CALL grib_get(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
2851
2852! Keys for rotated grids (checked through missing values)
2853CALL grib_get_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
2854 this%grid%proj%rotated%longitude_south_pole)
2855CALL grib_get_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
2856 this%grid%proj%rotated%latitude_south_pole)
2857CALL grib_get_dmiss(gaid,'angleOfRotationInDegrees', &
2858 this%grid%proj%rotated%angle_rotation)
2859
2860! Keys for stretched grids (checked through missing values)
2861! units must be verified, still experimental in grib_api
2862! # TODO: Is it a float? Is it signed?
2863IF (editionnumber == 1) THEN
2864 CALL grib_get_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
2865 this%grid%proj%stretched%longitude_stretch_pole)
2866 CALL grib_get_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
2867 this%grid%proj%stretched%latitude_stretch_pole)
2868 CALL grib_get_dmiss(gaid,'stretchingFactor', &
2869 this%grid%proj%stretched%stretch_factor)
2870ELSE IF (editionnumber == 2) THEN
2871 CALL grib_get_dmiss(gaid,'longitudeOfThePoleOfStretching', &
2872 this%grid%proj%stretched%longitude_stretch_pole)
2873 CALL grib_get_dmiss(gaid,'latitudeOfThePoleOfStretching', &
2874 this%grid%proj%stretched%latitude_stretch_pole)
2875 CALL grib_get_dmiss(gaid,'stretchingFactor', &
2876 this%grid%proj%stretched%stretch_factor)
2878 this%grid%proj%stretched%stretch_factor = this%grid%proj%stretched%stretch_factor*1.0d-6
2879ENDIF
2880
2881! Projection-dependent keys
2882SELECT CASE (this%grid%proj%proj_type)
2883
2884! Keys for spherical coordinate systems
2885CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
2886
2887 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
2888 CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
2889 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
2890 CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
2891
2892! longitudes are sometimes wrongly coded even in grib2 and even by the
2893! Metoffice!
2894! longitudeOfFirstGridPointInDegrees = 354.911;
2895! longitudeOfLastGridPointInDegrees = 363.311;
2896 CALL long_reset_0_360(lofirst)
2897 CALL long_reset_0_360(lolast)
2898
2899 IF (iscansnegatively == 0) THEN
2900 this%grid%grid%xmin = lofirst
2901 this%grid%grid%xmax = lolast
2902 ELSE
2903 this%grid%grid%xmax = lofirst
2904 this%grid%grid%xmin = lolast
2905 ENDIF
2906 IF (jscanspositively == 0) THEN
2907 this%grid%grid%ymax = lafirst
2908 this%grid%grid%ymin = lalast
2909 ELSE
2910 this%grid%grid%ymin = lafirst
2911 this%grid%grid%ymax = lalast
2912 ENDIF
2913
2914! reset longitudes in order to have a Cartesian plane
2915 IF (this%grid%grid%xmax < this%grid%grid%xmin) &
2916 this%grid%grid%xmin = this%grid%grid%xmin - 360.d0
2917
2918! compute dx and dy (should we get them from grib?)
2919 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
2920
2921! Keys for polar projections
2922CASE ('polar_stereographic', 'lambert', 'albers')
2923
2924 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
2925 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
2926! latin1/latin2 may be missing (e.g. stereographic)
2927 CALL grib_get_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
2928 CALL grib_get_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
2929#ifdef DEBUG
2931 "griddim_import_gribapi, latin1/2 "// &
2932 trim(to_char(this%grid%proj%polar%latin1))//" "// &
2933 trim(to_char(this%grid%proj%polar%latin2)))
2934#endif
2935! projection center flag, aka hemisphere
2936 CALL grib_get(gaid,'projectionCenterFlag',&
2937 this%grid%proj%polar%projection_center_flag, ierr)
2938 IF (ierr /= grib_success) THEN ! try center/centre
2939 CALL grib_get(gaid,'projectionCentreFlag',&
2940 this%grid%proj%polar%projection_center_flag)
2941 ENDIF
2942
2943 IF (iand(this%grid%proj%polar%projection_center_flag,64) == 1) THEN
2945 "griddim_import_gribapi, bi-polar projections not supported")
2946 CALL raise_error()
2947 ENDIF
2948! line of view, aka central meridian
2949 CALL grib_get(gaid,'LoVInDegrees',this%grid%proj%lov)
2950#ifdef DEBUG
2952 "griddim_import_gribapi, central meridian "//trim(to_char(this%grid%proj%lov)))
2953#endif
2954
2955! latitude at which dx and dy are valid
2956 IF (editionnumber == 1) THEN
2957! ECMWF (gribex/grib_api) says: Grid lengths are in metres, at the
2958! 60-degree parallel nearest to the pole on the projection plane.
2959! IF (IAND(this%projection_center_flag, 128) == 0) THEN
2960! this%grid%proj%polar%lad = 60.D0
2961! ELSE
2962! this%grid%proj%polar%lad = -60.D0
2963! ENDIF
2964! WMO says: Grid lengths are in units of metres, at the secant cone
2965! intersection parallel nearest to the pole on the projection plane.
2966 this%grid%proj%polar%lad = this%grid%proj%polar%latin1
2967 ELSE IF (editionnumber == 2) THEN
2968 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
2969 ENDIF
2970#ifdef DEBUG
2972 "griddim_import_gribapi, lad "//trim(to_char(this%grid%proj%polar%lad)))
2973#endif
2974
2975! compute projected extremes from lon and lat of first point
2976 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
2977 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
2978 CALL long_reset_0_360(lofirst)
2979 CALL long_reset_to_cart_closest(this%grid%proj%lov, lofirst)
2980#ifdef DEBUG
2982 "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
2984 "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
2985#endif
2986
2988 IF (iscansnegatively == 0) THEN
2989 this%grid%grid%xmin = x1
2990 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
2991 ELSE
2992 this%grid%grid%xmax = x1
2993 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
2994 ENDIF
2995 IF (jscanspositively == 0) THEN
2996 this%grid%grid%ymax = y1
2997 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
2998 ELSE
2999 this%grid%grid%ymin = y1
3000 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
3001 ENDIF
3002! keep these values for personal pleasure
3003 this%grid%proj%polar%lon1 = lofirst
3004 this%grid%proj%polar%lat1 = lafirst
3005
3006! Keys for equatorial projections
3007CASE ('mercator')
3008
3009 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
3010 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
3011 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3012 this%grid%proj%lov = 0.0d0 ! not in grib
3013
3014 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3015 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3016
3018 IF (iscansnegatively == 0) THEN
3019 this%grid%grid%xmin = x1
3020 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
3021 ELSE
3022 this%grid%grid%xmax = x1
3023 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
3024 ENDIF
3025 IF (jscanspositively == 0) THEN
3026 this%grid%grid%ymax = y1
3027 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
3028 ELSE
3029 this%grid%grid%ymin = y1
3030 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
3031 ENDIF
3032
3033 IF (editionnumber == 2) THEN
3034 CALL grib_get(gaid,'orientationOfTheGridInDegrees',orient)
3035 IF (orient /= 0.0d0) THEN
3037 "griddim_import_gribapi, Mercator grid orientation != 0 not supported")
3038 CALL raise_error()
3039 ENDIF
3040 ENDIF
3041
3042#ifdef DEBUG
3045 "griddim_import_gribapi, unprojected first point "//t2c(lofirst)//" "// &
3046 t2c(lafirst))
3047
3048 CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3049 CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3052 "griddim_import_gribapi, extremes from grib "//t2c(x1)//" "//t2c(y1))
3054 "griddim_import_gribapi, extremes from proj "//t2c(this%grid%grid%xmin)// &
3055 " "//t2c(this%grid%grid%ymin)//" "//t2c(this%grid%grid%xmax)//" "// &
3056 t2c(this%grid%grid%ymax))
3057#endif
3058
3059CASE ('UTM')
3060
3061 CALL grib_get(gaid,'zone',zone)
3062
3063 CALL grib_get(gaid,'datum',datum)
3064 IF (datum == 0) THEN
3065 CALL grib_get(gaid,'referenceLongitude',reflon)
3066 CALL grib_get(gaid,'falseEasting',this%grid%proj%xoff)
3067 CALL grib_get(gaid,'falseNorthing',this%grid%proj%yoff)
3069 ELSE
3071 CALL raise_fatal_error()
3072 ENDIF
3073
3074 CALL grib_get(gaid,'eastingOfFirstGridPoint',lofirst)
3075 CALL grib_get(gaid,'eastingOfLastGridPoint',lolast)
3076 CALL grib_get(gaid,'northingOfFirstGridPoint',lafirst)
3077 CALL grib_get(gaid,'northingOfLastGridPoint',lalast)
3078
3079 IF (iscansnegatively == 0) THEN
3080 this%grid%grid%xmin = lofirst
3081 this%grid%grid%xmax = lolast
3082 ELSE
3083 this%grid%grid%xmax = lofirst
3084 this%grid%grid%xmin = lolast
3085 ENDIF
3086 IF (jscanspositively == 0) THEN
3087 this%grid%grid%ymax = lafirst
3088 this%grid%grid%ymin = lalast
3089 ELSE
3090 this%grid%grid%ymin = lafirst
3091 this%grid%grid%ymax = lalast
3092 ENDIF
3093
3094! compute dx and dy (should we get them from grib?)
3095 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3096
3097CASE default
3099 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
3100 CALL raise_error()
3101
3102END SELECT
3103
3104CONTAINS
3105! utilities routines for grib_api, is there a better place?
3106SUBROUTINE grib_get_dmiss(gaid, key, value)
3107INTEGER,INTENT(in) :: gaid
3108CHARACTER(len=*),INTENT(in) :: key
3109DOUBLE PRECISION,INTENT(out) :: value
3110
3111INTEGER :: ierr
3112
3113CALL grib_get(gaid, key, value, ierr)
3114IF (ierr /= grib_success) value = dmiss
3115
3116END SUBROUTINE grib_get_dmiss
3117
3118SUBROUTINE grib_get_imiss(gaid, key, value)
3119INTEGER,INTENT(in) :: gaid
3120CHARACTER(len=*),INTENT(in) :: key
3121INTEGER,INTENT(out) :: value
3122
3123INTEGER :: ierr
3124
3125CALL grib_get(gaid, key, value, ierr)
3126IF (ierr /= grib_success) value = imiss
3127
3128END SUBROUTINE grib_get_imiss
3129
3130
3131SUBROUTINE griddim_import_ellipsoid(this, gaid)
3132TYPE(griddim_def),INTENT(inout) :: this
3133INTEGER,INTENT(in) :: gaid
3134
3135INTEGER :: shapeofearth, iv, is
3136DOUBLE PRECISION :: r1, r2
3137
3138IF (editionnumber == 2) THEN
3139 CALL grib_get(gaid, 'shapeOfTheEarth', shapeofearth)
3140 SELECT CASE(shapeofearth)
3141 CASE(0) ! spherical
3143 CASE(1) ! spherical generic
3144 CALL grib_get(gaid, 'scaleFactorOfRadiusOfSphericalEarth', is)
3145 CALL grib_get(gaid, 'scaledValueOfRadiusOfSphericalEarth', iv)
3146 r1 = dble(iv) / 10**is
3148 CASE(2) ! iau65
3150 CASE(3,7) ! ellipsoidal generic
3151 CALL grib_get(gaid, 'scaleFactorOfEarthMajorAxis', is)
3152 CALL grib_get(gaid, 'scaledValueOfEarthMajorAxis', iv)
3153 r1 = dble(iv) / 10**is
3154 CALL grib_get(gaid, 'scaleFactorOfEarthMinorAxis', is)
3155 CALL grib_get(gaid, 'scaledValueOfEarthMinorAxis', iv)
3156 r2 = dble(iv) / 10**is
3157 IF (shapeofearth == 3) THEN ! km->m
3158 r1 = r1*1000.0d0
3159 r2 = r2*1000.0d0
3160 ENDIF
3161 IF (abs(r1) < 1.0d-6) THEN ! suspicious data read from grib
3163 'read from grib, going on with spherical Earth but the results may be wrong')
3165 ELSE
3167 ENDIF
3168 CASE(4) ! iag-grs80
3170 CASE(5) ! wgs84
3172 CASE(6) ! spherical
3174! CASE(7) ! google earth-like?
3175 CASE default
3177 t2c(shapeofearth)//' not supported in grib2')
3178 CALL raise_fatal_error()
3179
3180 END SELECT
3181
3182ELSE
3183
3184 CALL grib_get(gaid, 'earthIsOblate', shapeofearth)
3185 IF (shapeofearth == 0) THEN ! spherical
3187 ELSE ! iau65
3189 ENDIF
3190
3191ENDIF
3192
3193END SUBROUTINE griddim_import_ellipsoid
3194
3195
3196END SUBROUTINE griddim_import_gribapi
3197
3198
3199! grib_api driver
3200SUBROUTINE griddim_export_gribapi(this, gaid)
3201USE grib_api
3202TYPE(griddim_def),INTENT(in) :: this ! griddim object
3203INTEGER, INTENT(inout) :: gaid ! grib_api id of the grib loaded in memory to export
3204
3205INTEGER :: EditionNumber, iScansNegatively, jScansPositively, nv, pvl, zone, ierr
3206DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, reflon
3207DOUBLE PRECISION :: sdx, sdy, ratio, tol, ex, ey
3208
3209
3210! Generic keys
3211CALL grib_get(gaid,'GRIBEditionNumber', editionnumber)
3212! the following required since eccodes
3213IF (editionnumber == 2) CALL grib_set(gaid,'shapeOfTheEarth', 6)
3214CALL grib_set(gaid,'typeOfGrid', this%grid%proj%proj_type)
3215#ifdef DEBUG
3217 "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
3218#endif
3219
3220IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
3221 CALL grib_set(gaid,'numberOfDataPoints', this%dim%nx)
3222! reenable when possible
3223! CALL griddim_export_ellipsoid(this, gaid)
3224 RETURN
3225ENDIF
3226
3227
3228! Edition dependent setup
3229IF (editionnumber == 1) THEN
3230 ratio = 1.d3
3231ELSE IF (editionnumber == 2) THEN
3232 ratio = 1.d6
3233ELSE
3234 ratio = 0.0d0 ! signal error?!
3235ENDIF
3236
3237! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
3238CALL griddim_export_ellipsoid(this, gaid)
3239CALL grib_set(gaid, 'Ni', this%dim%nx)
3240CALL grib_set(gaid, 'Nj', this%dim%ny)
3241CALL grib_set(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
3242
3243CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
3244CALL grib_get(gaid,'jScansPositively',jscanspositively)
3245
3246! Keys for rotated grids (checked through missing values and/or error code)
3247!SELECT CASE (this%grid%proj%proj_type)
3248!CASE ('rotated_ll', 'stretched_rotated_ll', 'polar_stereographic', 'lambert', 'albers')
3249CALL grib_set_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
3250 this%grid%proj%rotated%longitude_south_pole, 0.0d0)
3251CALL grib_set_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
3252 this%grid%proj%rotated%latitude_south_pole, -90.0d0)
3253IF (editionnumber == 1) THEN
3254 CALL grib_set_dmiss(gaid,'angleOfRotationInDegrees', &
3255 this%grid%proj%rotated%angle_rotation, 0.0d0)
3256ELSE IF (editionnumber == 2)THEN
3257 CALL grib_set_dmiss(gaid,'angleOfRotationOfProjectionInDegrees', &
3258 this%grid%proj%rotated%angle_rotation, 0.0d0)
3259ENDIF
3260
3261! Keys for stretched grids (checked through missing values and/or error code)
3262! units must be verified, still experimental in grib_api
3263! # TODO: Is it a float? Is it signed?
3264IF (editionnumber == 1) THEN
3265 CALL grib_set_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
3266 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
3267 CALL grib_set_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
3268 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
3269 CALL grib_set_dmiss(gaid,'stretchingFactor', &
3270 this%grid%proj%stretched%stretch_factor, 1.0d0)
3271ELSE IF (editionnumber == 2) THEN
3272 CALL grib_set_dmiss(gaid,'longitudeOfThePoleOfStretching', &
3273 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
3274 CALL grib_set_dmiss(gaid,'latitudeOfThePoleOfStretching', &
3275 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
3276 CALL grib_set_dmiss(gaid,'stretchingFactor', &
3277 this%grid%proj%stretched%stretch_factor, 1.0d6, 1.0d6)
3278ENDIF
3279
3280! Projection-dependent keys
3281SELECT CASE (this%grid%proj%proj_type)
3282
3283! Keys for sphaerical coordinate systems
3284CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
3285
3286 IF (iscansnegatively == 0) THEN
3287 lofirst = this%grid%grid%xmin
3288 lolast = this%grid%grid%xmax
3289 ELSE
3290 lofirst = this%grid%grid%xmax
3291 lolast = this%grid%grid%xmin
3292 ENDIF
3293 IF (jscanspositively == 0) THEN
3294 lafirst = this%grid%grid%ymax
3295 lalast = this%grid%grid%ymin
3296 ELSE
3297 lafirst = this%grid%grid%ymin
3298 lalast = this%grid%grid%ymax
3299 ENDIF
3300
3301! reset lon in standard grib 2 definition [0,360]
3302 IF (editionnumber == 1) THEN
3303 CALL long_reset_m180_360(lofirst)
3304 CALL long_reset_m180_360(lolast)
3305 ELSE IF (editionnumber == 2) THEN
3306 CALL long_reset_0_360(lofirst)
3307 CALL long_reset_0_360(lolast)
3308 ENDIF
3309
3310 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3311 CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3312 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3313 CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3314
3315! test relative coordinate truncation error with respect to tol
3316! tol should be tuned
3317 sdx = this%grid%grid%dx*ratio
3318 sdy = this%grid%grid%dy*ratio
3319! error is computed relatively to the whole grid extension
3320 ex = abs(((sdx - nint(sdx))/ratio*this%dim%nx)/ &
3321 (this%grid%grid%xmax-this%grid%grid%xmin))
3322 ey = abs(((sdy - nint(sdy))/ratio*this%dim%ny)/ &
3323 (this%grid%grid%ymax-this%grid%grid%ymin))
3324 tol = 1.0d-3
3325 IF (ex > tol .OR. ey > tol) THEN
3326#ifdef DEBUG
3328 "griddim_export_gribapi, error on x "//&
3329 trim(to_char(ex))//"/"//t2c(tol))
3331 "griddim_export_gribapi, error on y "//&
3332 trim(to_char(ey))//"/"//t2c(tol))
3333#endif
3334! previous frmula relative to a single grid step
3335! tol = 1.0d0/ratio
3336! IF (ABS(NINT(sdx)/sdx - 1.0d0) > tol .OR. ABS(NINT(sdy)/sdy - 1.0d0) > tol) THEN
3337!#ifdef DEBUG
3338! CALL l4f_category_log(this%category,L4F_DEBUG, &
3339! "griddim_export_gribapi, dlon relative error: "//&
3340! TRIM(to_char(ABS(NINT(sdx)/sdx - 1.0d0)))//">"//TRIM(to_char(tol)))
3341! CALL l4f_category_log(this%category,L4F_DEBUG, &
3342! "griddim_export_gribapi, dlat relative error: "//&
3343! TRIM(to_char(ABS(NINT(sdy)/sdy - 1.0d0)))//">"//TRIM(to_char(tol)))
3344!#endif
3346 "griddim_export_gribapi, increments not given: inaccurate!")
3347 CALL grib_set_missing(gaid,'Di')
3348 CALL grib_set_missing(gaid,'Dj')
3349 CALL grib_set(gaid,'ijDirectionIncrementGiven',0)
3350 ELSE
3351#ifdef DEBUG
3353 "griddim_export_gribapi, setting increments: "// &
3354 trim(to_char(this%grid%grid%dx))//' '//trim(to_char(this%grid%grid%dy)))
3355#endif
3356 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
3357 CALL grib_set(gaid,'iDirectionIncrement',nint(this%grid%grid%dx*ratio))
3358 CALL grib_set(gaid,'jDirectionIncrement',nint(this%grid%grid%dy*ratio))
3359! this does not work in grib_set
3360! CALL grib_set(gaid,'iDirectionIncrementInDegrees',this%grid%grid%dx)
3361! CALL grib_set(gaid,'jDirectionIncrementInDegrees',this%grid%grid%dy)
3362 ENDIF
3363
3364! Keys for polar projections
3365CASE ('polar_stereographic', 'lambert', 'albers')
3366! increments are required
3367 CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
3368 CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
3369 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
3370! latin1/latin2 may be missing (e.g. stereographic)
3371 CALL grib_set_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
3372 CALL grib_set_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
3373! projection center flag, aka hemisphere
3374 CALL grib_set(gaid,'projectionCenterFlag',&
3375 this%grid%proj%polar%projection_center_flag, ierr)
3376 IF (ierr /= grib_success) THEN ! try center/centre
3377 CALL grib_set(gaid,'projectionCentreFlag',&
3378 this%grid%proj%polar%projection_center_flag)
3379 ENDIF
3380
3381
3382! line of view, aka central meridian
3383 CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
3384! latitude at which dx and dy are valid
3385 IF (editionnumber == 2) THEN
3386 CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3387 ENDIF
3388
3389! compute lon and lat of first point from projected extremes
3390 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
3391 lofirst, lafirst)
3392! reset lon in standard grib 2 definition [0,360]
3393 IF (editionnumber == 1) THEN
3394 CALL long_reset_m180_360(lofirst)
3395 ELSE IF (editionnumber == 2) THEN
3396 CALL long_reset_0_360(lofirst)
3397 ENDIF
3398 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3399 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3400
3401! Keys for equatorial projections
3402CASE ('mercator')
3403
3404! increments are required
3405 CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
3406 CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
3407 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
3408
3409! line of view, aka central meridian, not in grib
3410! CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
3411! latitude at which dx and dy are valid
3412 CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3413
3414! compute lon and lat of first and last points from projected extremes
3415 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
3416 lofirst, lafirst)
3417 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3418 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3419 CALL get_unproj_last(this, iscansnegatively, jscanspositively, &
3420 lolast, lalast)
3421 CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3422 CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3423
3424 IF (editionnumber == 2) THEN
3425 CALL grib_set(gaid,'orientationOfTheGridInDegrees',0.)
3426 ENDIF
3427
3428CASE ('UTM')
3429
3430 CALL grib_set(gaid,'datum',0)
3432 CALL grib_set(gaid,'referenceLongitude',nint(reflon/1.0d6))
3433 CALL grib_set(gaid,'falseEasting',this%grid%proj%xoff)
3434 CALL grib_set(gaid,'falseNorthing',this%grid%proj%yoff)
3435
3436 CALL grib_set(gaid,'iDirectionIncrement',this%grid%grid%dx)
3437 CALL grib_set(gaid,'jDirectionIncrement',this%grid%grid%dy)
3438
3439!res/scann ??
3440
3441 CALL grib_set(gaid,'zone',zone)
3442
3443 IF (iscansnegatively == 0) THEN
3444 lofirst = this%grid%grid%xmin
3445 lolast = this%grid%grid%xmax
3446 ELSE
3447 lofirst = this%grid%grid%xmax
3448 lolast = this%grid%grid%xmin
3449 ENDIF
3450 IF (jscanspositively == 0) THEN
3451 lafirst = this%grid%grid%ymax
3452 lalast = this%grid%grid%ymin
3453 ELSE
3454 lafirst = this%grid%grid%ymin
3455 lalast = this%grid%grid%ymax
3456 ENDIF
3457
3458 CALL grib_set(gaid,'eastingOfFirstGridPoint',lofirst)
3459 CALL grib_set(gaid,'eastingOfLastGridPoint',lolast)
3460 CALL grib_set(gaid,'northingOfFirstGridPoint',lafirst)
3461 CALL grib_set(gaid,'northingOfLastGridPoint',lalast)
3462
3463CASE default
3465 "griddim_export_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
3466 CALL raise_error()
3467
3468END SELECT
3469
3470! hack for position of vertical coordinate parameters
3471! buggy in grib_api
3472IF (editionnumber == 1) THEN
3473! CALL grib_get(gaid,"PVPresent",pvp) ! alias, probably useless
3474 CALL grib_get(gaid,"NV",nv)
3475#ifdef DEBUG
3477 trim(to_char(nv))//" vertical coordinate parameters")
3478#endif
3479
3480 IF (nv == 0) THEN
3481 pvl = 255
3482 ELSE
3483 SELECT CASE (this%grid%proj%proj_type)
3484 CASE ('regular_ll') ! check whether "29-32 Set to zero (reserved)" are required
3485 pvl = 33
3486 CASE ('polar_stereographic')
3487 pvl = 33
3488 CASE ('rotated_ll', 'stretched_ll', 'lambert', 'albers')
3489 pvl = 43
3490 CASE ('stretched_rotated_ll')
3491 pvl = 43
3492 CASE DEFAULT
3493 pvl = 43 !?
3494 END SELECT
3495 ENDIF
3496
3497 CALL grib_set(gaid,"pvlLocation",pvl)
3498#ifdef DEBUG
3500 trim(to_char(pvl))//" as vertical coordinate parameter location")
3501#endif
3502
3503ENDIF
3504
3505
3506CONTAINS
3507! utilities routines for grib_api, is there a better place?
3508SUBROUTINE grib_set_dmiss(gaid, key, val, default, factor)
3509INTEGER,INTENT(in) :: gaid
3510CHARACTER(len=*),INTENT(in) :: key
3511DOUBLE PRECISION,INTENT(in) :: val
3512DOUBLE PRECISION,INTENT(in),OPTIONAL :: default
3513DOUBLE PRECISION,INTENT(in),OPTIONAL :: factor
3514
3515INTEGER :: ierr
3516
3518 IF (PRESENT(factor)) THEN
3519 CALL grib_set(gaid, key, val*factor, ierr)
3520 ELSE
3521 CALL grib_set(gaid, key, val, ierr)
3522 ENDIF
3523ELSE IF (PRESENT(default)) THEN
3524 CALL grib_set(gaid, key, default, ierr)
3525ENDIF
3526
3527END SUBROUTINE grib_set_dmiss
3528
3529SUBROUTINE grib_set_imiss(gaid, key, value, default)
3530INTEGER,INTENT(in) :: gaid
3531CHARACTER(len=*),INTENT(in) :: key
3532INTEGER,INTENT(in) :: value
3533INTEGER,INTENT(in),OPTIONAL :: default
3534
3535INTEGER :: ierr
3536
3538 CALL grib_set(gaid, key, value, ierr)
3539ELSE IF (PRESENT(default)) THEN
3540 CALL grib_set(gaid, key, default, ierr)
3541ENDIF
3542
3543END SUBROUTINE grib_set_imiss
3544
3545SUBROUTINE griddim_export_ellipsoid(this, gaid)
3546TYPE(griddim_def),INTENT(in) :: this
3547INTEGER,INTENT(in) :: gaid
3548
3549INTEGER :: ellips_type, ierr
3550DOUBLE PRECISION :: r1, r2, f
3551
3553
3554IF (editionnumber == 2) THEN
3555
3556! why does it produce a message even with ierr?
3557 CALL grib_set_missing(gaid, 'scaleFactorOfRadiusOfSphericalEarth', ierr)
3558 CALL grib_set_missing(gaid, 'scaledValueOfRadiusOfSphericalEarth', ierr)
3559 CALL grib_set_missing(gaid, 'scaleFactorOfEarthMajorAxis', ierr)
3560 CALL grib_set_missing(gaid, 'scaledValueOfEarthMajorAxis', ierr)
3561 CALL grib_set_missing(gaid, 'scaleFactorOfEarthMinorAxis', ierr)
3562 CALL grib_set_missing(gaid, 'scaledValueOfEarthMinorAxis', ierr)
3563
3564 SELECT CASE(ellips_type)
3565 CASE(ellips_grs80) ! iag-grs80
3566 CALL grib_set(gaid, 'shapeOfTheEarth', 4)
3567 CASE(ellips_wgs84) ! wgs84
3568 CALL grib_set(gaid, 'shapeOfTheEarth', 5)
3569 CASE default
3570 IF (f == 0.0d0) THEN ! spherical Earth
3571 IF (r1 == 6367470.0d0) THEN ! spherical
3572 CALL grib_set(gaid, 'shapeOfTheEarth', 0)
3573 ELSE IF (r1 == 6371229.0d0) THEN ! spherical
3574 CALL grib_set(gaid, 'shapeOfTheEarth', 6)
3575 ELSE ! spherical generic
3576 CALL grib_set(gaid, 'shapeOfTheEarth', 1)
3577 CALL grib_set(gaid, 'scaleFactorOfRadiusOfSphericalEarth', 2)
3578 CALL grib_set(gaid, 'scaledValueOfRadiusOfSphericalEarth', int(r1*100.0d0))
3579 ENDIF
3580 ELSE ! ellipsoidal
3581 IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
3582 CALL grib_set(gaid, 'shapeOfTheEarth', 2)
3583 ELSE ! ellipsoidal generic
3584 CALL grib_set(gaid, 'shapeOfTheEarth', 3)
3585 r2 = r1*(1.0d0 - f)
3586 CALL grib_set(gaid, 'scaleFactorOfEarthMajorAxis', 5)
3587 CALL grib_set(gaid, 'scaledValueOfEarthMajorAxis', &
3588 int(r1*100.0d0))
3589 CALL grib_set(gaid, 'scaleFactorOfEarthMinorAxis', 5)
3590 CALL grib_set(gaid, 'scaledValueOfEarthMinorAxis', &
3591 int(r2*100.0d0))
3592 ENDIF
3593 ENDIF
3594 END SELECT
3595
3596ELSE
3597
3598 IF (r1 == 6367470.0d0 .AND. f == 0.0d0) THEN ! spherical
3599 CALL grib_set(gaid, 'earthIsOblate', 0)
3600 ELSE IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
3601 CALL grib_set(gaid, 'earthIsOblate', 1)
3602 ELSE IF (f == 0.0d0) THEN ! generic spherical
3603 CALL grib_set(gaid, 'earthIsOblate', 0)
3605 ¬ supported in grib 1, coding with default radius of 6367470 m')
3606 ELSE ! generic ellipsoidal
3607 CALL grib_set(gaid, 'earthIsOblate', 1)
3609 ¬ supported in grib 1, coding with default iau65 ellipsoid')
3610 ENDIF
3611
3612ENDIF
3613
3614END SUBROUTINE griddim_export_ellipsoid
3615
3616SUBROUTINE get_unproj_first(this, iScansNegatively, jScansPositively, &
3617 loFirst, laFirst)
3618TYPE(griddim_def),INTENT(in) :: this ! griddim object
3619INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
3620DOUBLE PRECISION,INTENT(out) :: loFirst, laFirst
3621
3622! compute lon and lat of first point from projected extremes
3623IF (iscansnegatively == 0) THEN
3624 IF (jscanspositively == 0) THEN
3626 ELSE
3628 ENDIF
3629ELSE
3630 IF (jscanspositively == 0) THEN
3632 ELSE
3634 ENDIF
3635ENDIF
3636! use the values kept for personal pleasure ?
3637! loFirst = this%grid%proj%polar%lon1
3638! laFirst = this%grid%proj%polar%lat1
3639END SUBROUTINE get_unproj_first
3640
3641SUBROUTINE get_unproj_last(this, iScansNegatively, jScansPositively, &
3642 loLast, laLast)
3643TYPE(griddim_def),INTENT(in) :: this ! griddim object
3644INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
3645DOUBLE PRECISION,INTENT(out) :: loLast, laLast
3646
3647! compute lon and lat of last point from projected extremes
3648IF (iscansnegatively == 0) THEN
3649 IF (jscanspositively == 0) THEN
3651 ELSE
3653 ENDIF
3654ELSE
3655 IF (jscanspositively == 0) THEN
3657 ELSE
3659 ENDIF
3660ENDIF
3661! use the values kept for personal pleasure ?
3662! loLast = this%grid%proj%polar%lon?
3663! laLast = this%grid%proj%polar%lat?
3664END SUBROUTINE get_unproj_last
3665
3666END SUBROUTINE griddim_export_gribapi
3667#endif
3668
3669
3670#ifdef HAVE_LIBGDAL
3671! gdal driver
3672SUBROUTINE griddim_import_gdal(this, gdalid, gdal_options)
3673USE gdal
3674TYPE(griddim_def),INTENT(inout) :: this ! griddim object
3675TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal rasterband pointer
3676TYPE(gdal_file_id_options),INTENT(in) :: gdal_options
3677
3678TYPE(gdaldataseth) :: hds
3679REAL(kind=c_double) :: geotrans(6) !, invgeotrans(6), x1, y1, x2, y2
3680INTEGER(kind=c_int) :: offsetx, offsety
3681INTEGER :: ier
3682
3683hds = gdalgetbanddataset(gdalid) ! go back to dataset
3684ier = gdalgetgeotransform(hds, geotrans)
3685
3686IF (ier /= 0) THEN
3688 'griddim_import_gdal: error in accessing gdal dataset')
3689 CALL raise_error()
3690 RETURN
3691ENDIF
3692IF (geotrans(3) /= 0.0_c_double .OR. geotrans(5) /= 0.0_c_double) THEN ! dataset has not a diagonal transformation
3694 'griddim_import_gdal: dataset has a non-diagonal transformation matrix, unsupported')
3695 CALL raise_error()
3696 RETURN
3697ENDIF
3698
3699CALL gdaldatasetbbsize_f(hds, gdal_options%xmin, gdal_options%ymin, &
3700 gdal_options%xmax, gdal_options%ymax, &
3701 this%dim%nx, this%dim%ny, offsetx, offsety, &
3702 this%grid%grid%xmin, this%grid%grid%ymin, this%grid%grid%xmax, this%grid%grid%ymax)
3703
3704IF (this%dim%nx == 0 .OR. this%dim%ny == 0) THEN
3706 'griddim_import_gdal: requested bounding box '//t2c(gdal_options%xmin)//','// &
3707 t2c(gdal_options%ymin)//','//t2c(gdal_options%xmax)//','//&
3708 t2c(gdal_options%ymax))
3710 'determines an empty dataset '//t2c(this%dim%nx)//'x'//t2c(this%dim%ny))
3711ENDIF
3712
3713! get grid corners
3714!CALL gdalapplygeotransform(geotrans, 0.5_c_double, 0.5_c_double, x1, y1)
3715!CALL gdalapplygeotransform(geotrans, gdalgetrasterbandxsize(gdalid)-0.5_c_double, &
3716! gdalgetrasterbandysize(gdalid)-0.5_c_double, x2, y2)
3717
3718!IF (geotrans(3) == 0.0_c_double .AND. geotrans(5) == 0.0_c_double) THEN ! transformation is diagonal, no transposing
3719! this%dim%nx = gdalgetrasterbandxsize(gdalid)
3720! this%dim%ny = gdalgetrasterbandysize(gdalid)
3721! this%grid%grid%xmin = MIN(x1, x2)
3722! this%grid%grid%xmax = MAX(x1, x2)
3723! this%grid%grid%ymin = MIN(y1, y2)
3724! this%grid%grid%ymax = MAX(y1, y2)
3725!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
3726!
3727! this%dim%nx = gdalgetrasterbandysize(gdalid)
3728! this%dim%ny = gdalgetrasterbandxsize(gdalid)
3729! this%grid%grid%xmin = MIN(y1, y2)
3730! this%grid%grid%xmax = MAX(y1, y2)
3731! this%grid%grid%ymin = MIN(x1, x2)
3732! this%grid%grid%ymax = MAX(x1, x2)
3733!
3734!ELSE ! transformation is a rotation, not supported
3735!ENDIF
3736
3737this%grid%proj%proj_type = 'regular_ll' ! forced, only one supported (check?)
3738
3739CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3740this%grid%grid%component_flag = 0
3741
3742END SUBROUTINE griddim_import_gdal
3743#endif
3744
3745
3747SUBROUTINE griddim_display(this)
3748TYPE(griddim_def),INTENT(in) :: this
3749
3750print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
3751
3755
3756print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
3757
3758END SUBROUTINE griddim_display
3759
3760
3761#define VOL7D_POLY_TYPE TYPE(grid_def)
3762#define VOL7D_POLY_TYPES _grid
3763#include "array_utilities_inc.F90"
3764#undef VOL7D_POLY_TYPE
3765#undef VOL7D_POLY_TYPES
3766
3767#define VOL7D_POLY_TYPE TYPE(griddim_def)
3768#define VOL7D_POLY_TYPES _griddim
3769#include "array_utilities_inc.F90"
3770#undef VOL7D_POLY_TYPE
3771#undef VOL7D_POLY_TYPES
3772
3773
3785SUBROUTINE griddim_wind_unrot(this, rot_mat)
3787TYPE(griddim_def), INTENT(in) :: this
3788DOUBLE PRECISION, POINTER :: rot_mat(:,:,:)
3789
3790DOUBLE PRECISION :: dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j
3791DOUBLE PRECISION :: lat_factor
3792INTEGER :: i, j, ip1, im1, jp1, jm1;
3793
3794IF (this%dim%nx <= 0 .OR. this%dim%ny <= 0 .OR. &
3795 .NOT. ASSOCIATED(this%dim%lon) .OR. .NOT. ASSOCIATED(this%dim%lat)) THEN
3796 CALL l4f_category_log(this%category, l4f_error, 'rotate_uv must be called after correct init of griddim object')
3797 NULLIFY(rot_mat)
3798 RETURN
3799ENDIF
3800
3801ALLOCATE(rot_mat(this%dim%nx, this%dim%ny, 4))
3802
3803DO j = 1, this%dim%ny
3804 jp1 = min(j+1, this%dim%ny)
3805 jm1 = max(j-1, 1)
3806 DO i = 1, this%dim%nx
3807 ip1 = min(i+1, this%dim%nx)
3808 im1 = max(i-1, 1)
3809
3810 dlat_i = this%dim%lat(ip1,j) - this%dim%lat(im1,j)
3811 dlat_j = this%dim%lat(i,jp1) - this%dim%lat(i,jm1)
3812
3813 dlon_i = this%dim%lon(ip1,j) - this%dim%lon(im1,j)
3814! if ( dlon_i > pi ) dlon_i -= 2*pi;
3815! if ( dlon_i < (-pi) ) dlon_i += 2*pi;
3816 dlon_j = this%dim%lon(i,jp1) - this%dim%lon(i,jm1)
3817! if ( dlon_j > pi ) dlon_j -= 2*pi;
3818! if ( dlon_j < (-pi) ) dlon_j += 2*pi;
3819
3820! check whether this is really necessary !!!!
3821 lat_factor = cos(degrad*this%dim%lat(i,j))
3822 dlon_i = dlon_i * lat_factor
3823 dlon_j = dlon_j * lat_factor
3824
3825 dist_i = sqrt(dlon_i*dlon_i + dlat_i*dlat_i)
3826 dist_j = sqrt(dlon_j*dlon_j + dlat_j*dlat_j)
3827
3828 IF (dist_i > 0.d0) THEN
3829 rot_mat(i,j,1) = dlon_i/dist_i
3830 rot_mat(i,j,2) = dlat_i/dist_i
3831 ELSE
3832 rot_mat(i,j,1) = 0.0d0
3833 rot_mat(i,j,2) = 0.0d0
3834 ENDIF
3835 IF (dist_j > 0.d0) THEN
3836 rot_mat(i,j,3) = dlon_j/dist_j
3837 rot_mat(i,j,4) = dlat_j/dist_j
3838 ELSE
3839 rot_mat(i,j,3) = 0.0d0
3840 rot_mat(i,j,4) = 0.0d0
3841 ENDIF
3842
3843 ENDDO
3844ENDDO
3845
3846END SUBROUTINE griddim_wind_unrot
3847
3848
3849! compute zoom indices from geographical zoom coordinates
3850SUBROUTINE griddim_zoom_coord(this, ilon, ilat, flon, flat, ix, iy, fx, fy)
3851TYPE(griddim_def),INTENT(in) :: this
3852DOUBLE PRECISION,INTENT(in) :: ilon, ilat, flon, flat
3853INTEGER,INTENT(out) :: ix, iy, fx, fy
3854
3855DOUBLE PRECISION :: ix1, iy1, fx1, fy1
3856
3857! compute projected coordinates of vertices of desired lonlat rectangle
3860
3861CALL griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
3862
3863END SUBROUTINE griddim_zoom_coord
3864
3865
3866! compute zoom indices from projected zoom coordinates
3867SUBROUTINE griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
3868TYPE(griddim_def),INTENT(in) :: this
3869DOUBLE PRECISION,INTENT(in) :: ix1, iy1, fx1, fy1
3870INTEGER,INTENT(out) :: ix, iy, fx, fy
3871
3872INTEGER :: lix, liy, lfx, lfy
3873
3874! compute projected indices
3875lix = nint((ix1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
3876liy = nint((iy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
3877lfx = nint((fx1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
3878lfy = nint((fy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
3879! swap projected indices, in case grid is "strongly rotated"
3880ix = min(lix, lfx)
3881fx = max(lix, lfx)
3882iy = min(liy, lfy)
3883fy = max(liy, lfy)
3884
3885END SUBROUTINE griddim_zoom_projcoord
3886
3887
3891SUBROUTINE long_reset_0_360(lon)
3892DOUBLE PRECISION,INTENT(inout) :: lon
3893
3895DO WHILE(lon < 0.0d0)
3896 lon = lon + 360.0d0
3897END DO
3898DO WHILE(lon >= 360.0d0)
3899 lon = lon - 360.0d0
3900END DO
3901
3902END SUBROUTINE long_reset_0_360
3903
3904
3908SUBROUTINE long_reset_m180_360(lon)
3909DOUBLE PRECISION,INTENT(inout) :: lon
3910
3912DO WHILE(lon < -180.0d0)
3913 lon = lon + 360.0d0
3914END DO
3915DO WHILE(lon >= 360.0d0)
3916 lon = lon - 360.0d0
3917END DO
3918
3919END SUBROUTINE long_reset_m180_360
3920
3921
3925!SUBROUTINE long_reset_m90_270(lon)
3926!DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
3927!
3928!IF (.NOT.c_e(lon)) RETURN
3929!DO WHILE(lon < -90.0D0)
3930! lon = lon + 360.0D0
3931!END DO
3932!DO WHILE(lon >= 270.0D0)
3933! lon = lon - 360.0D0
3934!END DO
3935!
3936!END SUBROUTINE long_reset_m90_270
3937
3938
3942SUBROUTINE long_reset_m180_180(lon)
3943DOUBLE PRECISION,INTENT(inout) :: lon
3944
3946DO WHILE(lon < -180.0d0)
3947 lon = lon + 360.0d0
3948END DO
3949DO WHILE(lon >= 180.0d0)
3950 lon = lon - 360.0d0
3951END DO
3952
3953END SUBROUTINE long_reset_m180_180
3954
3955
3956SUBROUTINE long_reset_to_cart_closest(lon, lonref)
3957DOUBLE PRECISION,INTENT(inout) :: lon
3958DOUBLE PRECISION,INTENT(in) :: lonref
3959
3961IF (abs(lon-lonref) < 180.0d0) RETURN ! nothing to do
3962lon = lon - nint((lon-lonref)/360.0d0)*360.0d0 ! se non e` vera e` ben trovata
3963
3964END SUBROUTINE long_reset_to_cart_closest
3965
3966
3968
Compute forward coordinate transformation from geographical system to projected system. Definition grid_class.F90:308 Read the object from a formatted or unformatted file. Definition grid_class.F90:334 Compute backward coordinate transformation from projected system to geographical system. Definition grid_class.F90:314 Write the object on a formatted or unformatted file. Definition grid_class.F90:329 Emit log message for a category with specific priority. Definition log4fortran.F90:457 Module for describing geographically referenced regular grids. Definition grid_class.F90:237 Module for defining the extension and coordinates of a rectangular georeferenced grid. Definition grid_dim_class.F90:211 This module defines an abstract interface to different drivers for access to files containing gridded... Definition grid_id_class.F90:249 Definitions of constants and functions for working with missing values. Definition missing_values.f90:50 Module for quickly interpreting the OPTIONAL parameters passed to a subprogram. Definition optional_values.f90:28 This object, mainly for internal use, describes a grid on a geographical projection,... Definition grid_class.F90:257 This object completely describes a grid on a geographic projection. Definition grid_class.F90:270 Derived type describing the extension of a grid and the geographical coordinates of each point. Definition grid_dim_class.F90:221 |