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