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