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