libsim Versione 7.1.11

◆ pack_distinct_grid()

type(grid_def) function, dimension(dim) pack_distinct_grid ( type(grid_def), dimension(:), intent(in)  vect,
integer, intent(in)  dim,
logical, dimension(:), intent(in), optional  mask,
logical, intent(in), optional  back 
)
private

compatta gli elementi distinti di vect in un array

Definizione alla linea 2064 del file grid_class.F90.

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

Generated with Doxygen.