libsim Versione 7.1.11
|
◆ pack_distinct_grid()
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
2115use geo_proj_class
2117use grid_rect_class
2123implicit none
2124
2125
2126character (len=255),parameter:: subcategory="grid_class"
2127
2128
2135 private
2136 type(geo_proj) :: proj
2137 type(grid_rect) :: grid
2138 integer :: category = 0
2140
2141
2148 type(grid_def) :: grid
2149 type(grid_dim) :: dim
2150 integer :: category = 0
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
2170 MODULE PROCEDURE griddim_init
2171END INTERFACE
2172
2175 MODULE PROCEDURE griddim_delete
2176END INTERFACE
2177
2180 MODULE PROCEDURE griddim_copy
2181END INTERFACE
2182
2186 MODULE PROCEDURE griddim_coord_proj!, griddim_proj
2187END INTERFACE
2188
2192 MODULE PROCEDURE griddim_coord_unproj, griddim_unproj
2193END INTERFACE
2194
2197 MODULE PROCEDURE griddim_get_val
2198END INTERFACE
2199
2202 MODULE PROCEDURE griddim_set_val
2203END INTERFACE
2204
2207 MODULE PROCEDURE griddim_write_unit
2208END INTERFACE
2209
2212 MODULE PROCEDURE griddim_read_unit
2213END INTERFACE
2214
2216INTERFACE import
2217 MODULE PROCEDURE griddim_import_grid_id
2218END INTERFACE
2219
2222 MODULE PROCEDURE griddim_export_grid_id
2223END INTERFACE
2224
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
2250 griddim_zoom_coord, griddim_zoom_projcoord, &
2254PUBLIC OPERATOR(==),OPERATOR(/=)
2255PUBLIC count_distinct, pack_distinct, count_and_pack_distinct, &
2256 map_distinct, map_inv_distinct,index
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
2329#endif
2330
2331END SUBROUTINE griddim_init
2332
2333
2335SUBROUTINE griddim_delete(this)
2336TYPE(griddim_def),INTENT(inout) :: this
2337
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
2356
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
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
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
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)
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
2486
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
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
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
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
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
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
2631IF (unit == geo_proj_unit_meter) THEN ! it is a plane projection
2633 IF (PRESENT(lonref)) THEN
2634 CALL long_reset_to_cart_closest(lov, lonref)
2636 ENDIF
2637
2638ELSE IF (unit == geo_proj_unit_degree) THEN ! it is a spheric projection
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)
2648! now reset rotated coordinates around zero
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
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
2777
2778#ifdef HAVE_LIBGRIBAPI
2779gaid = grid_id_get_gaid(ingrid_id)
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)
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
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)
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
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
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
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
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
2988 "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
2990 "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
2991#endif
2992
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
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
3043 "griddim_import_gribapi, Mercator grid orientation != 0 not supported")
3044 CALL raise_error()
3045 ENDIF
3046 ENDIF
3047
3048#ifdef 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)
3058 "griddim_import_gribapi, extremes from grib "//t2c(x1)//" "//t2c(y1))
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)
3075 ELSE
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
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
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
3154 CASE(2) ! iau65
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
3169 'read from grib, going on with spherical Earth but the results may be wrong')
3171 ELSE
3173 ENDIF
3174 CASE(4) ! iag-grs80
3176 CASE(5) ! wgs84
3178 CASE(6) ! spherical
3180! CASE(7) ! google earth-like?
3181 CASE default
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
3193 ELSE ! iau65
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
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
3334 "griddim_export_gribapi, error on x "//&
3335 trim(to_char(ex))//"/"//t2c(tol))
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
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
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)
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
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
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
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
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
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
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)
3611 ¬ supported in grib 1, coding with default radius of 6367470 m')
3612 ELSE ! generic ellipsoidal
3613 CALL grib_set(gaid, 'earthIsOblate', 1)
3615 ¬ 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
3632 ELSE
3634 ENDIF
3635ELSE
3636 IF (jscanspositively == 0) THEN
3638 ELSE
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
3657 ELSE
3659 ENDIF
3660ELSE
3661 IF (jscanspositively == 0) THEN
3663 ELSE
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
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
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
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))
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
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
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
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
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
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
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
3974
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 |