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