libsim Versione 7.2.0
|
◆ pack_distinct_griddim()
compatta gli elementi distinti di vect in un array Definizione alla linea 2572 del file grid_class.F90. 2574! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2575! authors:
2576! Davide Cesari <dcesari@arpa.emr.it>
2577! Paolo Patruno <ppatruno@arpa.emr.it>
2578
2579! This program is free software; you can redistribute it and/or
2580! modify it under the terms of the GNU General Public License as
2581! published by the Free Software Foundation; either version 2 of
2582! the License, or (at your option) any later version.
2583
2584! This program is distributed in the hope that it will be useful,
2585! but WITHOUT ANY WARRANTY; without even the implied warranty of
2586! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2587! GNU General Public License for more details.
2588
2589! You should have received a copy of the GNU General Public License
2590! along with this program. If not, see <http://www.gnu.org/licenses/>.
2591#include "config.h"
2592
2623use geo_proj_class
2625use grid_rect_class
2631implicit none
2632
2633
2634character (len=255),parameter:: subcategory="grid_class"
2635
2636
2643 private
2644 type(geo_proj) :: proj
2645 type(grid_rect) :: grid
2646 integer :: category = 0
2648
2649
2656 type(grid_def) :: grid
2657 type(grid_dim) :: dim
2658 integer :: category = 0
2660
2661
2665INTERFACE OPERATOR (==)
2666 MODULE PROCEDURE grid_eq, griddim_eq
2667END INTERFACE
2668
2672INTERFACE OPERATOR (/=)
2673 MODULE PROCEDURE grid_ne, griddim_ne
2674END INTERFACE
2675
2678 MODULE PROCEDURE griddim_init
2679END INTERFACE
2680
2683 MODULE PROCEDURE griddim_delete
2684END INTERFACE
2685
2688 MODULE PROCEDURE griddim_copy
2689END INTERFACE
2690
2694 MODULE PROCEDURE griddim_coord_proj!, griddim_proj
2695END INTERFACE
2696
2700 MODULE PROCEDURE griddim_coord_unproj, griddim_unproj
2701END INTERFACE
2702
2705 MODULE PROCEDURE griddim_get_val
2706END INTERFACE
2707
2710 MODULE PROCEDURE griddim_set_val
2711END INTERFACE
2712
2715 MODULE PROCEDURE griddim_write_unit
2716END INTERFACE
2717
2720 MODULE PROCEDURE griddim_read_unit
2721END INTERFACE
2722
2724INTERFACE import
2725 MODULE PROCEDURE griddim_import_grid_id
2726END INTERFACE
2727
2730 MODULE PROCEDURE griddim_export_grid_id
2731END INTERFACE
2732
2735 MODULE PROCEDURE griddim_display
2736END INTERFACE
2737
2738#define VOL7D_POLY_TYPE TYPE(grid_def)
2739#define VOL7D_POLY_TYPES _grid
2740#include "array_utilities_pre.F90"
2741#undef VOL7D_POLY_TYPE
2742#undef VOL7D_POLY_TYPES
2743
2744#define VOL7D_POLY_TYPE TYPE(griddim_def)
2745#define VOL7D_POLY_TYPES _griddim
2746#include "array_utilities_pre.F90"
2747#undef VOL7D_POLY_TYPE
2748#undef VOL7D_POLY_TYPES
2749
2750INTERFACE wind_unrot
2751 MODULE PROCEDURE griddim_wind_unrot
2752END INTERFACE
2753
2754
2755PRIVATE
2756
2758 griddim_zoom_coord, griddim_zoom_projcoord, &
2762PUBLIC OPERATOR(==),OPERATOR(/=)
2763PUBLIC count_distinct, pack_distinct, count_and_pack_distinct, &
2764 map_distinct, map_inv_distinct,index
2766PUBLIC griddim_central_lon, griddim_set_central_lon
2767CONTAINS
2768
2770SUBROUTINE griddim_init(this, nx, ny, &
2771 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
2772 proj_type, lov, zone, xoff, yoff, &
2773 longitude_south_pole, latitude_south_pole, angle_rotation, &
2774 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
2775 latin1, latin2, lad, projection_center_flag, &
2776 ellips_smaj_axis, ellips_flatt, ellips_type, &
2777 categoryappend)
2778TYPE(griddim_def),INTENT(inout) :: this
2779INTEGER,INTENT(in),OPTIONAL :: nx
2780INTEGER,INTENT(in),OPTIONAL :: ny
2781DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin
2782DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmax
2783DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymin
2784DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymax
2785DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx
2786DOUBLE PRECISION,INTENT(in),OPTIONAL :: dy
2789INTEGER,INTENT(in),OPTIONAL :: component_flag
2790CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
2791DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
2792INTEGER,INTENT(in),OPTIONAL :: zone
2793DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
2794DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
2795DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
2796DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
2797DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
2798DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
2799DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
2800DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
2801DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
2802DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
2803DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
2804INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
2805DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
2806DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
2807INTEGER,INTENT(in),OPTIONAL :: ellips_type
2808CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2809
2810CHARACTER(len=512) :: a_name
2811
2812IF (PRESENT(categoryappend)) THEN
2813 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
2814 trim(categoryappend))
2815ELSE
2816 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2817ENDIF
2818this%category=l4f_category_get(a_name)
2819
2820! geographical projection
2821this%grid%proj = geo_proj_new( &
2822 proj_type=proj_type, lov=lov, zone=zone, xoff=xoff, yoff=yoff, &
2823 longitude_south_pole=longitude_south_pole, &
2824 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
2825 longitude_stretch_pole=longitude_stretch_pole, &
2826 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
2827 latin1=latin1, latin2=latin2, lad=lad, projection_center_flag=projection_center_flag, &
2828 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, ellips_type=ellips_type)
2829! grid extension
2830this%grid%grid = grid_rect_new( &
2831 xmin, xmax, ymin, ymax, dx, dy, component_flag)
2832! grid size
2833this%dim = grid_dim_new(nx, ny)
2834
2835#ifdef DEBUG
2837#endif
2838
2839END SUBROUTINE griddim_init
2840
2841
2843SUBROUTINE griddim_delete(this)
2844TYPE(griddim_def),INTENT(inout) :: this
2845
2849
2850call l4f_category_delete(this%category)
2851
2852END SUBROUTINE griddim_delete
2853
2854
2856SUBROUTINE griddim_copy(this, that, categoryappend)
2857TYPE(griddim_def),INTENT(in) :: this
2858TYPE(griddim_def),INTENT(out) :: that
2859CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2860
2861CHARACTER(len=512) :: a_name
2862
2864
2868
2869! new category
2870IF (PRESENT(categoryappend)) THEN
2871 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
2872 trim(categoryappend))
2873ELSE
2874 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2875ENDIF
2876that%category=l4f_category_get(a_name)
2877
2878END SUBROUTINE griddim_copy
2879
2880
2883ELEMENTAL SUBROUTINE griddim_coord_proj(this, lon, lat, x, y)
2884TYPE(griddim_def),INTENT(in) :: this
2886DOUBLE PRECISION,INTENT(in) :: lon, lat
2888DOUBLE PRECISION,INTENT(out) :: x, y
2889
2891
2892END SUBROUTINE griddim_coord_proj
2893
2894
2897ELEMENTAL SUBROUTINE griddim_coord_unproj(this, x, y, lon, lat)
2898TYPE(griddim_def),INTENT(in) :: this
2900DOUBLE PRECISION,INTENT(in) :: x, y
2902DOUBLE PRECISION,INTENT(out) :: lon, lat
2903
2905
2906END SUBROUTINE griddim_coord_unproj
2907
2908
2909! Computes and sets the grid parameters required to compute
2910! coordinates of grid points in the projected system.
2911! probably meaningless
2912!SUBROUTINE griddim_proj(this)
2913!TYPE(griddim_def),INTENT(inout) :: this !< definition of projection and grid
2914!
2915!CALL proj(this, this%dim%lon(1,1), this%dim%lat(1,1), &
2916! this%grid%grid%xmin, this%grid%grid%ymin)
2917!
2918!CALL proj(this, this%dim%lon(this%dim%nx,this%dim%ny), &
2919! this%dim%lat(this%dim%nx,this%dim%ny), &
2920! this%grid%grid%xmax, this%grid%grid%ymax)
2921!
2922!END SUBROUTINE griddim_proj
2923
2931SUBROUTINE griddim_unproj(this)
2932TYPE(griddim_def),INTENT(inout) :: this
2933
2935CALL alloc(this%dim)
2936CALL griddim_unproj_internal(this)
2937
2938END SUBROUTINE griddim_unproj
2939
2940! internal subroutine needed for allocating automatic arrays
2941SUBROUTINE griddim_unproj_internal(this)
2942TYPE(griddim_def),INTENT(inout) ::this ! definition of projection and grid
2943
2944DOUBLE PRECISION :: x(this%dim%nx,this%dim%ny), y(this%dim%nx,this%dim%ny)
2945
2946CALL grid_rect_coordinates(this%grid%grid, x, y)
2948
2949END SUBROUTINE griddim_unproj_internal
2950
2951
2953SUBROUTINE griddim_get_val(this, nx, ny, &
2954 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
2955 proj, proj_type, lov, zone, xoff, yoff, &
2956 longitude_south_pole, latitude_south_pole, angle_rotation, &
2957 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
2958 latin1, latin2, lad, projection_center_flag, &
2959 ellips_smaj_axis, ellips_flatt, ellips_type)
2960TYPE(griddim_def),INTENT(in) :: this
2961INTEGER,INTENT(out),OPTIONAL :: nx
2962INTEGER,INTENT(out),OPTIONAL :: ny
2964DOUBLE PRECISION,INTENT(out),OPTIONAL :: xmin, xmax, ymin, ymax
2966DOUBLE PRECISION,INTENT(out),OPTIONAL :: dx, dy
2969INTEGER,INTENT(out),OPTIONAL :: component_flag
2970TYPE(geo_proj),INTENT(out),OPTIONAL :: proj
2971CHARACTER(len=*),INTENT(out),OPTIONAL :: proj_type
2972DOUBLE PRECISION,INTENT(out),OPTIONAL :: lov
2973INTEGER,INTENT(out),OPTIONAL :: zone
2974DOUBLE PRECISION,INTENT(out),OPTIONAL :: xoff
2975DOUBLE PRECISION,INTENT(out),OPTIONAL :: yoff
2976DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_south_pole
2977DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_south_pole
2978DOUBLE PRECISION,INTENT(out),OPTIONAL :: angle_rotation
2979DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_stretch_pole
2980DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_stretch_pole
2981DOUBLE PRECISION,INTENT(out),OPTIONAL :: stretch_factor
2982DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin1
2983DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin2
2984DOUBLE PRECISION,INTENT(out),OPTIONAL :: lad
2985INTEGER,INTENT(out),OPTIONAL :: projection_center_flag
2986DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_smaj_axis
2987DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_flatt
2988INTEGER,INTENT(out),OPTIONAL :: ellips_type
2989
2990IF (PRESENT(nx)) nx = this%dim%nx
2991IF (PRESENT(ny)) ny = this%dim%ny
2992
2994
2996 xoff=xoff, yoff=yoff, &
2997 longitude_south_pole=longitude_south_pole, &
2998 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
2999 longitude_stretch_pole=longitude_stretch_pole, &
3000 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
3001 latin1=latin1, latin2=latin2, lad=lad, &
3002 projection_center_flag=projection_center_flag, &
3003 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
3004 ellips_type=ellips_type)
3005
3007 xmin, xmax, ymin, ymax, dx, dy, component_flag)
3008
3009END SUBROUTINE griddim_get_val
3010
3011
3013SUBROUTINE griddim_set_val(this, nx, ny, &
3014 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
3015 proj_type, lov, zone, xoff, yoff, &
3016 longitude_south_pole, latitude_south_pole, angle_rotation, &
3017 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
3018 latin1, latin2, lad, projection_center_flag, &
3019 ellips_smaj_axis, ellips_flatt, ellips_type)
3020TYPE(griddim_def),INTENT(inout) :: this
3021INTEGER,INTENT(in),OPTIONAL :: nx
3022INTEGER,INTENT(in),OPTIONAL :: ny
3024DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin, xmax, ymin, ymax
3026DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx, dy
3029INTEGER,INTENT(in),OPTIONAL :: component_flag
3030CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
3031DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
3032INTEGER,INTENT(in),OPTIONAL :: zone
3033DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
3034DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
3035DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
3036DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
3037DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
3038DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
3039DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
3040DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
3041DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
3042DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
3043DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
3044INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
3045DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
3046DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
3047INTEGER,INTENT(in),OPTIONAL :: ellips_type
3048
3049IF (PRESENT(nx)) this%dim%nx = nx
3050IF (PRESENT(ny)) this%dim%ny = ny
3051
3053 xoff=xoff, yoff=yoff, longitude_south_pole=longitude_south_pole, &
3054 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
3055 longitude_stretch_pole=longitude_stretch_pole, &
3056 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
3057 latin1=latin1, latin2=latin2, lad=lad, &
3058 projection_center_flag=projection_center_flag, &
3059 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
3060 ellips_type=ellips_type)
3061
3063 xmin, xmax, ymin, ymax, dx, dy, component_flag)
3064
3065END SUBROUTINE griddim_set_val
3066
3067
3072SUBROUTINE griddim_read_unit(this, unit)
3073TYPE(griddim_def),INTENT(out) :: this
3074INTEGER, INTENT(in) :: unit
3075
3076
3080
3081END SUBROUTINE griddim_read_unit
3082
3083
3088SUBROUTINE griddim_write_unit(this, unit)
3089TYPE(griddim_def),INTENT(in) :: this
3090INTEGER, INTENT(in) :: unit
3091
3092
3096
3097END SUBROUTINE griddim_write_unit
3098
3099
3103FUNCTION griddim_central_lon(this) RESULT(lon)
3104TYPE(griddim_def),INTENT(inout) :: this
3105
3106DOUBLE PRECISION :: lon
3107
3108CALL griddim_pistola_central_lon(this, lon)
3109
3110END FUNCTION griddim_central_lon
3111
3112
3116SUBROUTINE griddim_set_central_lon(this, lonref)
3117TYPE(griddim_def),INTENT(inout) :: this
3118DOUBLE PRECISION,INTENT(in) :: lonref
3119
3120DOUBLE PRECISION :: lon
3121
3122CALL griddim_pistola_central_lon(this, lon, lonref)
3123
3124END SUBROUTINE griddim_set_central_lon
3125
3126
3127! internal subroutine for performing tasks common to the prevous two
3128SUBROUTINE griddim_pistola_central_lon(this, lon, lonref)
3129TYPE(griddim_def),INTENT(inout) :: this ! grid descriptor
3130DOUBLE PRECISION,INTENT(inout) :: lon ! central longitude
3131DOUBLE PRECISION,INTENT(in),OPTIONAL :: lonref ! reference longitude
3132
3133INTEGER :: unit
3134DOUBLE PRECISION :: lonsp, latsp, londelta, lov, lonrot
3135CHARACTER(len=80) :: ptype
3136
3137lon = dmiss
3139IF (unit == geo_proj_unit_meter) THEN ! it is a plane projection
3141 IF (PRESENT(lonref)) THEN
3142 CALL long_reset_to_cart_closest(lov, lonref)
3144 ENDIF
3145
3146ELSE IF (unit == geo_proj_unit_degree) THEN ! it is a spheric projection
3148 longitude_south_pole=lonsp, latitude_south_pole=latsp)
3149 SELECT CASE(ptype)
3150 CASE('rotated_ll','stretched_rotated_ll') ! use origin of rotated system
3151 IF (latsp < 0.0d0) THEN
3152 lon = lonsp
3153 IF (PRESENT(lonref)) THEN
3154 CALL long_reset_to_cart_closest(lon, lonref)
3156! now reset rotated coordinates around zero
3158 lonrot = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
3159 ENDIF
3160 londelta = lonrot
3161 CALL long_reset_to_cart_closest(londelta, 0.0d0)
3162 londelta = londelta - lonrot
3163 this%grid%grid%xmin = this%grid%grid%xmin + londelta
3164 this%grid%grid%xmax = this%grid%grid%xmax + londelta
3165 ENDIF
3166 ELSE ! this part to be checked
3167 lon = modulo(lonsp + 180.0d0, 360.0d0)
3168! IF (PRESENT(lonref)) THEN
3169! CALL long_reset_to_cart_closest(lov, lonref)
3170! CALL set_val(this%grid%proj, longitude_south_pole=lonref)
3171! ENDIF
3172 ENDIF
3173 CASE default ! use real grid limits
3175 lon = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
3176 ENDIF
3177 IF (PRESENT(lonref)) THEN
3178 londelta = lon
3179 CALL long_reset_to_cart_closest(londelta, lonref)
3180 londelta = londelta - lon
3181 this%grid%grid%xmin = this%grid%grid%xmin + londelta
3182 this%grid%grid%xmax = this%grid%grid%xmax + londelta
3183 ENDIF
3184 END SELECT
3185ENDIF
3186
3187END SUBROUTINE griddim_pistola_central_lon
3188
3189
3193SUBROUTINE griddim_gen_coord(this, x, y)
3194TYPE(griddim_def),INTENT(in) :: this
3195DOUBLE PRECISION,INTENT(out) :: x(:,:)
3196DOUBLE PRECISION,INTENT(out) :: y(:,:)
3197
3198
3199CALL grid_rect_coordinates(this%grid%grid, x, y)
3200
3201END SUBROUTINE griddim_gen_coord
3202
3203
3205SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
3206TYPE(griddim_def), INTENT(in) :: this
3207INTEGER,INTENT(in) :: nx
3208INTEGER,INTENT(in) :: ny
3209DOUBLE PRECISION,INTENT(out) :: dx
3210DOUBLE PRECISION,INTENT(out) :: dy
3211
3212CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
3213
3214END SUBROUTINE griddim_steps
3215
3216
3218SUBROUTINE griddim_setsteps(this)
3219TYPE(griddim_def), INTENT(inout) :: this
3220
3221CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3222
3223END SUBROUTINE griddim_setsteps
3224
3225
3226! TODO
3227! bisogna sviluppare gli altri operatori
3228ELEMENTAL FUNCTION grid_eq(this, that) RESULT(res)
3229TYPE(grid_def),INTENT(IN) :: this, that
3230LOGICAL :: res
3231
3232res = this%proj == that%proj .AND. &
3233 this%grid == that%grid
3234
3235END FUNCTION grid_eq
3236
3237
3238ELEMENTAL FUNCTION griddim_eq(this, that) RESULT(res)
3239TYPE(griddim_def),INTENT(IN) :: this, that
3240LOGICAL :: res
3241
3242res = this%grid == that%grid .AND. &
3243 this%dim == that%dim
3244
3245END FUNCTION griddim_eq
3246
3247
3248ELEMENTAL FUNCTION grid_ne(this, that) RESULT(res)
3249TYPE(grid_def),INTENT(IN) :: this, that
3250LOGICAL :: res
3251
3252res = .NOT.(this == that)
3253
3254END FUNCTION grid_ne
3255
3256
3257ELEMENTAL FUNCTION griddim_ne(this, that) RESULT(res)
3258TYPE(griddim_def),INTENT(IN) :: this, that
3259LOGICAL :: res
3260
3261res = .NOT.(this == that)
3262
3263END FUNCTION griddim_ne
3264
3265
3271SUBROUTINE griddim_import_grid_id(this, ingrid_id)
3272#ifdef HAVE_LIBGDAL
3273USE gdal
3274#endif
3275TYPE(griddim_def),INTENT(inout) :: this
3276TYPE(grid_id),INTENT(in) :: ingrid_id
3277
3278#ifdef HAVE_LIBGRIBAPI
3279INTEGER :: gaid
3280#endif
3281#ifdef HAVE_LIBGDAL
3282TYPE(gdalrasterbandh) :: gdalid
3283#endif
3285
3286#ifdef HAVE_LIBGRIBAPI
3287gaid = grid_id_get_gaid(ingrid_id)
3289#endif
3290#ifdef HAVE_LIBGDAL
3291gdalid = grid_id_get_gdalid(ingrid_id)
3292IF (gdalassociated(gdalid)) CALL griddim_import_gdal(this, gdalid, &
3293 grid_id_get_gdal_options(ingrid_id))
3294#endif
3295
3296END SUBROUTINE griddim_import_grid_id
3297
3298
3303SUBROUTINE griddim_export_grid_id(this, outgrid_id)
3304#ifdef HAVE_LIBGDAL
3305USE gdal
3306#endif
3307TYPE(griddim_def),INTENT(in) :: this
3308TYPE(grid_id),INTENT(inout) :: outgrid_id
3309
3310#ifdef HAVE_LIBGRIBAPI
3311INTEGER :: gaid
3312#endif
3313#ifdef HAVE_LIBGDAL
3314TYPE(gdalrasterbandh) :: gdalid
3315#endif
3316
3317#ifdef HAVE_LIBGRIBAPI
3318gaid = grid_id_get_gaid(outgrid_id)
3320#endif
3321#ifdef HAVE_LIBGDAL
3322gdalid = grid_id_get_gdalid(outgrid_id)
3323!IF (gdalassociated(gdalid)
3324! export for gdal not implemented, log?
3325#endif
3326
3327END SUBROUTINE griddim_export_grid_id
3328
3329
3330#ifdef HAVE_LIBGRIBAPI
3331! grib_api driver
3332SUBROUTINE griddim_import_gribapi(this, gaid)
3333USE grib_api
3334TYPE(griddim_def),INTENT(inout) :: this ! griddim object
3335INTEGER, INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
3336
3337DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, x1, y1, orient
3338INTEGER :: EditionNumber, iScansNegatively, jScansPositively, zone, datum, &
3339 reflon, ierr
3340
3341! Generic keys
3342CALL grib_get(gaid, 'typeOfGrid', this%grid%proj%proj_type)
3343#ifdef DEBUG
3345 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type))
3346#endif
3347CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
3348
3349IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
3350 CALL grib_get(gaid,'numberOfDataPoints', this%dim%nx)
3351 this%dim%ny = 1
3352 this%grid%grid%component_flag = 0
3353 CALL griddim_import_ellipsoid(this, gaid)
3354 RETURN
3355ENDIF
3356
3357! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
3358CALL grib_get(gaid, 'Ni', this%dim%nx)
3359CALL grib_get(gaid, 'Nj', this%dim%ny)
3360CALL griddim_import_ellipsoid(this, gaid) ! placed here, not valid for utm datum /= 1
3361
3362CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
3363CALL grib_get(gaid,'jScansPositively',jscanspositively)
3364CALL grib_get(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
3365
3366! Keys for rotated grids (checked through missing values)
3367CALL grib_get_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
3368 this%grid%proj%rotated%longitude_south_pole)
3369CALL grib_get_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
3370 this%grid%proj%rotated%latitude_south_pole)
3371CALL grib_get_dmiss(gaid,'angleOfRotationInDegrees', &
3372 this%grid%proj%rotated%angle_rotation)
3373
3374! Keys for stretched grids (checked through missing values)
3375! units must be verified, still experimental in grib_api
3376! # TODO: Is it a float? Is it signed?
3377IF (editionnumber == 1) THEN
3378 CALL grib_get_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
3379 this%grid%proj%stretched%longitude_stretch_pole)
3380 CALL grib_get_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
3381 this%grid%proj%stretched%latitude_stretch_pole)
3382 CALL grib_get_dmiss(gaid,'stretchingFactor', &
3383 this%grid%proj%stretched%stretch_factor)
3384ELSE IF (editionnumber == 2) THEN
3385 CALL grib_get_dmiss(gaid,'longitudeOfThePoleOfStretching', &
3386 this%grid%proj%stretched%longitude_stretch_pole)
3387 CALL grib_get_dmiss(gaid,'latitudeOfThePoleOfStretching', &
3388 this%grid%proj%stretched%latitude_stretch_pole)
3389 CALL grib_get_dmiss(gaid,'stretchingFactor', &
3390 this%grid%proj%stretched%stretch_factor)
3392 this%grid%proj%stretched%stretch_factor = this%grid%proj%stretched%stretch_factor*1.0d-6
3393ENDIF
3394
3395! Projection-dependent keys
3396SELECT CASE (this%grid%proj%proj_type)
3397
3398! Keys for spherical coordinate systems
3399CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
3400
3401 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3402 CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3403 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3404 CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3405
3406! longitudes are sometimes wrongly coded even in grib2 and even by the
3407! Metoffice!
3408! longitudeOfFirstGridPointInDegrees = 354.911;
3409! longitudeOfLastGridPointInDegrees = 363.311;
3410 CALL long_reset_0_360(lofirst)
3411 CALL long_reset_0_360(lolast)
3412
3413 IF (iscansnegatively == 0) THEN
3414 this%grid%grid%xmin = lofirst
3415 this%grid%grid%xmax = lolast
3416 ELSE
3417 this%grid%grid%xmax = lofirst
3418 this%grid%grid%xmin = lolast
3419 ENDIF
3420 IF (jscanspositively == 0) THEN
3421 this%grid%grid%ymax = lafirst
3422 this%grid%grid%ymin = lalast
3423 ELSE
3424 this%grid%grid%ymin = lafirst
3425 this%grid%grid%ymax = lalast
3426 ENDIF
3427
3428! reset longitudes in order to have a Cartesian plane
3429 IF (this%grid%grid%xmax < this%grid%grid%xmin) &
3430 this%grid%grid%xmin = this%grid%grid%xmin - 360.d0
3431
3432! compute dx and dy (should we get them from grib?)
3433 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3434
3435! Keys for polar projections
3436CASE ('polar_stereographic', 'lambert', 'albers')
3437
3438 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
3439 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
3440! latin1/latin2 may be missing (e.g. stereographic)
3441 CALL grib_get_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
3442 CALL grib_get_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
3443#ifdef DEBUG
3445 "griddim_import_gribapi, latin1/2 "// &
3446 trim(to_char(this%grid%proj%polar%latin1))//" "// &
3447 trim(to_char(this%grid%proj%polar%latin2)))
3448#endif
3449! projection center flag, aka hemisphere
3450 CALL grib_get(gaid,'projectionCenterFlag',&
3451 this%grid%proj%polar%projection_center_flag, ierr)
3452 IF (ierr /= grib_success) THEN ! try center/centre
3453 CALL grib_get(gaid,'projectionCentreFlag',&
3454 this%grid%proj%polar%projection_center_flag)
3455 ENDIF
3456
3457 IF (iand(this%grid%proj%polar%projection_center_flag,64) == 1) THEN
3459 "griddim_import_gribapi, bi-polar projections not supported")
3460 CALL raise_error()
3461 ENDIF
3462! line of view, aka central meridian
3463 CALL grib_get(gaid,'LoVInDegrees',this%grid%proj%lov)
3464#ifdef DEBUG
3466 "griddim_import_gribapi, central meridian "//trim(to_char(this%grid%proj%lov)))
3467#endif
3468
3469! latitude at which dx and dy are valid
3470 IF (editionnumber == 1) THEN
3471! ECMWF (gribex/grib_api) says: Grid lengths are in metres, at the
3472! 60-degree parallel nearest to the pole on the projection plane.
3473! IF (IAND(this%projection_center_flag, 128) == 0) THEN
3474! this%grid%proj%polar%lad = 60.D0
3475! ELSE
3476! this%grid%proj%polar%lad = -60.D0
3477! ENDIF
3478! WMO says: Grid lengths are in units of metres, at the secant cone
3479! intersection parallel nearest to the pole on the projection plane.
3480 this%grid%proj%polar%lad = this%grid%proj%polar%latin1
3481 ELSE IF (editionnumber == 2) THEN
3482 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3483 ENDIF
3484#ifdef DEBUG
3486 "griddim_import_gribapi, lad "//trim(to_char(this%grid%proj%polar%lad)))
3487#endif
3488
3489! compute projected extremes from lon and lat of first point
3490 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3491 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3492 CALL long_reset_0_360(lofirst)
3493 CALL long_reset_to_cart_closest(this%grid%proj%lov, lofirst)
3494#ifdef DEBUG
3496 "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
3498 "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
3499#endif
3500
3502 IF (iscansnegatively == 0) THEN
3503 this%grid%grid%xmin = x1
3504 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
3505 ELSE
3506 this%grid%grid%xmax = x1
3507 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
3508 ENDIF
3509 IF (jscanspositively == 0) THEN
3510 this%grid%grid%ymax = y1
3511 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
3512 ELSE
3513 this%grid%grid%ymin = y1
3514 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
3515 ENDIF
3516! keep these values for personal pleasure
3517 this%grid%proj%polar%lon1 = lofirst
3518 this%grid%proj%polar%lat1 = lafirst
3519
3520! Keys for equatorial projections
3521CASE ('mercator')
3522
3523 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
3524 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
3525 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3526 this%grid%proj%lov = 0.0d0 ! not in grib
3527
3528 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3529 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3530
3532 IF (iscansnegatively == 0) THEN
3533 this%grid%grid%xmin = x1
3534 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
3535 ELSE
3536 this%grid%grid%xmax = x1
3537 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
3538 ENDIF
3539 IF (jscanspositively == 0) THEN
3540 this%grid%grid%ymax = y1
3541 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
3542 ELSE
3543 this%grid%grid%ymin = y1
3544 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
3545 ENDIF
3546
3547 IF (editionnumber == 2) THEN
3548 CALL grib_get(gaid,'orientationOfTheGridInDegrees',orient)
3549 IF (orient /= 0.0d0) THEN
3551 "griddim_import_gribapi, Mercator grid orientation != 0 not supported")
3552 CALL raise_error()
3553 ENDIF
3554 ENDIF
3555
3556#ifdef DEBUG
3559 "griddim_import_gribapi, unprojected first point "//t2c(lofirst)//" "// &
3560 t2c(lafirst))
3561
3562 CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3563 CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3566 "griddim_import_gribapi, extremes from grib "//t2c(x1)//" "//t2c(y1))
3568 "griddim_import_gribapi, extremes from proj "//t2c(this%grid%grid%xmin)// &
3569 " "//t2c(this%grid%grid%ymin)//" "//t2c(this%grid%grid%xmax)//" "// &
3570 t2c(this%grid%grid%ymax))
3571#endif
3572
3573CASE ('UTM')
3574
3575 CALL grib_get(gaid,'zone',zone)
3576
3577 CALL grib_get(gaid,'datum',datum)
3578 IF (datum == 0) THEN
3579 CALL grib_get(gaid,'referenceLongitude',reflon)
3580 CALL grib_get(gaid,'falseEasting',this%grid%proj%xoff)
3581 CALL grib_get(gaid,'falseNorthing',this%grid%proj%yoff)
3583 ELSE
3585 CALL raise_fatal_error()
3586 ENDIF
3587
3588 CALL grib_get(gaid,'eastingOfFirstGridPoint',lofirst)
3589 CALL grib_get(gaid,'eastingOfLastGridPoint',lolast)
3590 CALL grib_get(gaid,'northingOfFirstGridPoint',lafirst)
3591 CALL grib_get(gaid,'northingOfLastGridPoint',lalast)
3592
3593 IF (iscansnegatively == 0) THEN
3594 this%grid%grid%xmin = lofirst
3595 this%grid%grid%xmax = lolast
3596 ELSE
3597 this%grid%grid%xmax = lofirst
3598 this%grid%grid%xmin = lolast
3599 ENDIF
3600 IF (jscanspositively == 0) THEN
3601 this%grid%grid%ymax = lafirst
3602 this%grid%grid%ymin = lalast
3603 ELSE
3604 this%grid%grid%ymin = lafirst
3605 this%grid%grid%ymax = lalast
3606 ENDIF
3607
3608! compute dx and dy (should we get them from grib?)
3609 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3610
3611CASE default
3613 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
3614 CALL raise_error()
3615
3616END SELECT
3617
3618CONTAINS
3619! utilities routines for grib_api, is there a better place?
3620SUBROUTINE grib_get_dmiss(gaid, key, value)
3621INTEGER,INTENT(in) :: gaid
3622CHARACTER(len=*),INTENT(in) :: key
3623DOUBLE PRECISION,INTENT(out) :: value
3624
3625INTEGER :: ierr
3626
3627CALL grib_get(gaid, key, value, ierr)
3628IF (ierr /= grib_success) value = dmiss
3629
3630END SUBROUTINE grib_get_dmiss
3631
3632SUBROUTINE grib_get_imiss(gaid, key, value)
3633INTEGER,INTENT(in) :: gaid
3634CHARACTER(len=*),INTENT(in) :: key
3635INTEGER,INTENT(out) :: value
3636
3637INTEGER :: ierr
3638
3639CALL grib_get(gaid, key, value, ierr)
3640IF (ierr /= grib_success) value = imiss
3641
3642END SUBROUTINE grib_get_imiss
3643
3644
3645SUBROUTINE griddim_import_ellipsoid(this, gaid)
3646TYPE(griddim_def),INTENT(inout) :: this
3647INTEGER,INTENT(in) :: gaid
3648
3649INTEGER :: shapeofearth, iv, is
3650DOUBLE PRECISION :: r1, r2
3651
3652IF (editionnumber == 2) THEN
3653 CALL grib_get(gaid, 'shapeOfTheEarth', shapeofearth)
3654 SELECT CASE(shapeofearth)
3655 CASE(0) ! spherical
3657 CASE(1) ! spherical generic
3658 CALL grib_get(gaid, 'scaleFactorOfRadiusOfSphericalEarth', is)
3659 CALL grib_get(gaid, 'scaledValueOfRadiusOfSphericalEarth', iv)
3660 r1 = dble(iv) / 10**is
3662 CASE(2) ! iau65
3664 CASE(3,7) ! ellipsoidal generic
3665 CALL grib_get(gaid, 'scaleFactorOfEarthMajorAxis', is)
3666 CALL grib_get(gaid, 'scaledValueOfEarthMajorAxis', iv)
3667 r1 = dble(iv) / 10**is
3668 CALL grib_get(gaid, 'scaleFactorOfEarthMinorAxis', is)
3669 CALL grib_get(gaid, 'scaledValueOfEarthMinorAxis', iv)
3670 r2 = dble(iv) / 10**is
3671 IF (shapeofearth == 3) THEN ! km->m
3672 r1 = r1*1000.0d0
3673 r2 = r2*1000.0d0
3674 ENDIF
3675 IF (abs(r1) < 1.0d-6) THEN ! suspicious data read from grib
3677 'read from grib, going on with spherical Earth but the results may be wrong')
3679 ELSE
3681 ENDIF
3682 CASE(4) ! iag-grs80
3684 CASE(5) ! wgs84
3686 CASE(6) ! spherical
3688! CASE(7) ! google earth-like?
3689 CASE default
3691 t2c(shapeofearth)//' not supported in grib2')
3692 CALL raise_fatal_error()
3693
3694 END SELECT
3695
3696ELSE
3697
3698 CALL grib_get(gaid, 'earthIsOblate', shapeofearth)
3699 IF (shapeofearth == 0) THEN ! spherical
3701 ELSE ! iau65
3703 ENDIF
3704
3705ENDIF
3706
3707END SUBROUTINE griddim_import_ellipsoid
3708
3709
3710END SUBROUTINE griddim_import_gribapi
3711
3712
3713! grib_api driver
3714SUBROUTINE griddim_export_gribapi(this, gaid)
3715USE grib_api
3716TYPE(griddim_def),INTENT(in) :: this ! griddim object
3717INTEGER, INTENT(inout) :: gaid ! grib_api id of the grib loaded in memory to export
3718
3719INTEGER :: EditionNumber, iScansNegatively, jScansPositively, nv, pvl, zone, ierr
3720DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, reflon
3721DOUBLE PRECISION :: sdx, sdy, ratio, tol, ex, ey
3722
3723
3724! Generic keys
3725CALL grib_get(gaid,'GRIBEditionNumber', editionnumber)
3726! the following required since eccodes
3727IF (editionnumber == 2) CALL grib_set(gaid,'shapeOfTheEarth', 6)
3728CALL grib_set(gaid,'typeOfGrid', this%grid%proj%proj_type)
3729#ifdef DEBUG
3731 "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
3732#endif
3733
3734IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
3735 CALL grib_set(gaid,'numberOfDataPoints', this%dim%nx)
3736! reenable when possible
3737! CALL griddim_export_ellipsoid(this, gaid)
3738 RETURN
3739ENDIF
3740
3741
3742! Edition dependent setup
3743IF (editionnumber == 1) THEN
3744 ratio = 1.d3
3745ELSE IF (editionnumber == 2) THEN
3746 ratio = 1.d6
3747ELSE
3748 ratio = 0.0d0 ! signal error?!
3749ENDIF
3750
3751! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
3752CALL griddim_export_ellipsoid(this, gaid)
3753CALL grib_set(gaid, 'Ni', this%dim%nx)
3754CALL grib_set(gaid, 'Nj', this%dim%ny)
3755CALL grib_set(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
3756
3757CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
3758CALL grib_get(gaid,'jScansPositively',jscanspositively)
3759
3760! Keys for rotated grids (checked through missing values and/or error code)
3761!SELECT CASE (this%grid%proj%proj_type)
3762!CASE ('rotated_ll', 'stretched_rotated_ll', 'polar_stereographic', 'lambert', 'albers')
3763CALL grib_set_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
3764 this%grid%proj%rotated%longitude_south_pole, 0.0d0)
3765CALL grib_set_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
3766 this%grid%proj%rotated%latitude_south_pole, -90.0d0)
3767IF (editionnumber == 1) THEN
3768 CALL grib_set_dmiss(gaid,'angleOfRotationInDegrees', &
3769 this%grid%proj%rotated%angle_rotation, 0.0d0)
3770ELSE IF (editionnumber == 2)THEN
3771 CALL grib_set_dmiss(gaid,'angleOfRotationOfProjectionInDegrees', &
3772 this%grid%proj%rotated%angle_rotation, 0.0d0)
3773ENDIF
3774
3775! Keys for stretched grids (checked through missing values and/or error code)
3776! units must be verified, still experimental in grib_api
3777! # TODO: Is it a float? Is it signed?
3778IF (editionnumber == 1) THEN
3779 CALL grib_set_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
3780 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
3781 CALL grib_set_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
3782 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
3783 CALL grib_set_dmiss(gaid,'stretchingFactor', &
3784 this%grid%proj%stretched%stretch_factor, 1.0d0)
3785ELSE IF (editionnumber == 2) THEN
3786 CALL grib_set_dmiss(gaid,'longitudeOfThePoleOfStretching', &
3787 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
3788 CALL grib_set_dmiss(gaid,'latitudeOfThePoleOfStretching', &
3789 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
3790 CALL grib_set_dmiss(gaid,'stretchingFactor', &
3791 this%grid%proj%stretched%stretch_factor, 1.0d6, 1.0d6)
3792ENDIF
3793
3794! Projection-dependent keys
3795SELECT CASE (this%grid%proj%proj_type)
3796
3797! Keys for sphaerical coordinate systems
3798CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
3799
3800 IF (iscansnegatively == 0) THEN
3801 lofirst = this%grid%grid%xmin
3802 lolast = this%grid%grid%xmax
3803 ELSE
3804 lofirst = this%grid%grid%xmax
3805 lolast = this%grid%grid%xmin
3806 ENDIF
3807 IF (jscanspositively == 0) THEN
3808 lafirst = this%grid%grid%ymax
3809 lalast = this%grid%grid%ymin
3810 ELSE
3811 lafirst = this%grid%grid%ymin
3812 lalast = this%grid%grid%ymax
3813 ENDIF
3814
3815! reset lon in standard grib 2 definition [0,360]
3816 IF (editionnumber == 1) THEN
3817 CALL long_reset_m180_360(lofirst)
3818 CALL long_reset_m180_360(lolast)
3819 ELSE IF (editionnumber == 2) THEN
3820 CALL long_reset_0_360(lofirst)
3821 CALL long_reset_0_360(lolast)
3822 ENDIF
3823
3824 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3825 CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3826 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3827 CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3828
3829! test relative coordinate truncation error with respect to tol
3830! tol should be tuned
3831 sdx = this%grid%grid%dx*ratio
3832 sdy = this%grid%grid%dy*ratio
3833! error is computed relatively to the whole grid extension
3834 ex = abs(((sdx - nint(sdx))/ratio*this%dim%nx)/ &
3835 (this%grid%grid%xmax-this%grid%grid%xmin))
3836 ey = abs(((sdy - nint(sdy))/ratio*this%dim%ny)/ &
3837 (this%grid%grid%ymax-this%grid%grid%ymin))
3838 tol = 1.0d-3
3839 IF (ex > tol .OR. ey > tol) THEN
3840#ifdef DEBUG
3842 "griddim_export_gribapi, error on x "//&
3843 trim(to_char(ex))//"/"//t2c(tol))
3845 "griddim_export_gribapi, error on y "//&
3846 trim(to_char(ey))//"/"//t2c(tol))
3847#endif
3848! previous frmula relative to a single grid step
3849! tol = 1.0d0/ratio
3850! IF (ABS(NINT(sdx)/sdx - 1.0d0) > tol .OR. ABS(NINT(sdy)/sdy - 1.0d0) > tol) THEN
3851!#ifdef DEBUG
3852! CALL l4f_category_log(this%category,L4F_DEBUG, &
3853! "griddim_export_gribapi, dlon relative error: "//&
3854! TRIM(to_char(ABS(NINT(sdx)/sdx - 1.0d0)))//">"//TRIM(to_char(tol)))
3855! CALL l4f_category_log(this%category,L4F_DEBUG, &
3856! "griddim_export_gribapi, dlat relative error: "//&
3857! TRIM(to_char(ABS(NINT(sdy)/sdy - 1.0d0)))//">"//TRIM(to_char(tol)))
3858!#endif
3860 "griddim_export_gribapi, increments not given: inaccurate!")
3861 CALL grib_set_missing(gaid,'Di')
3862 CALL grib_set_missing(gaid,'Dj')
3863 CALL grib_set(gaid,'ijDirectionIncrementGiven',0)
3864 ELSE
3865#ifdef DEBUG
3867 "griddim_export_gribapi, setting increments: "// &
3868 trim(to_char(this%grid%grid%dx))//' '//trim(to_char(this%grid%grid%dy)))
3869#endif
3870 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
3871 CALL grib_set(gaid,'iDirectionIncrement',nint(this%grid%grid%dx*ratio))
3872 CALL grib_set(gaid,'jDirectionIncrement',nint(this%grid%grid%dy*ratio))
3873! this does not work in grib_set
3874! CALL grib_set(gaid,'iDirectionIncrementInDegrees',this%grid%grid%dx)
3875! CALL grib_set(gaid,'jDirectionIncrementInDegrees',this%grid%grid%dy)
3876 ENDIF
3877
3878! Keys for polar projections
3879CASE ('polar_stereographic', 'lambert', 'albers')
3880! increments are required
3881 CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
3882 CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
3883 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
3884! latin1/latin2 may be missing (e.g. stereographic)
3885 CALL grib_set_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
3886 CALL grib_set_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
3887! projection center flag, aka hemisphere
3888 CALL grib_set(gaid,'projectionCenterFlag',&
3889 this%grid%proj%polar%projection_center_flag, ierr)
3890 IF (ierr /= grib_success) THEN ! try center/centre
3891 CALL grib_set(gaid,'projectionCentreFlag',&
3892 this%grid%proj%polar%projection_center_flag)
3893 ENDIF
3894
3895
3896! line of view, aka central meridian
3897 CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
3898! latitude at which dx and dy are valid
3899 IF (editionnumber == 2) THEN
3900 CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3901 ENDIF
3902
3903! compute lon and lat of first point from projected extremes
3904 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
3905 lofirst, lafirst)
3906! reset lon in standard grib 2 definition [0,360]
3907 IF (editionnumber == 1) THEN
3908 CALL long_reset_m180_360(lofirst)
3909 ELSE IF (editionnumber == 2) THEN
3910 CALL long_reset_0_360(lofirst)
3911 ENDIF
3912 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3913 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3914
3915! Keys for equatorial projections
3916CASE ('mercator')
3917
3918! increments are required
3919 CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
3920 CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
3921 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
3922
3923! line of view, aka central meridian, not in grib
3924! CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
3925! latitude at which dx and dy are valid
3926 CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3927
3928! compute lon and lat of first and last points from projected extremes
3929 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
3930 lofirst, lafirst)
3931 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3932 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3933 CALL get_unproj_last(this, iscansnegatively, jscanspositively, &
3934 lolast, lalast)
3935 CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3936 CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3937
3938 IF (editionnumber == 2) THEN
3939 CALL grib_set(gaid,'orientationOfTheGridInDegrees',0.)
3940 ENDIF
3941
3942CASE ('UTM')
3943
3944 CALL grib_set(gaid,'datum',0)
3946 CALL grib_set(gaid,'referenceLongitude',nint(reflon/1.0d6))
3947 CALL grib_set(gaid,'falseEasting',this%grid%proj%xoff)
3948 CALL grib_set(gaid,'falseNorthing',this%grid%proj%yoff)
3949
3950 CALL grib_set(gaid,'iDirectionIncrement',this%grid%grid%dx)
3951 CALL grib_set(gaid,'jDirectionIncrement',this%grid%grid%dy)
3952
3953!res/scann ??
3954
3955 CALL grib_set(gaid,'zone',zone)
3956
3957 IF (iscansnegatively == 0) THEN
3958 lofirst = this%grid%grid%xmin
3959 lolast = this%grid%grid%xmax
3960 ELSE
3961 lofirst = this%grid%grid%xmax
3962 lolast = this%grid%grid%xmin
3963 ENDIF
3964 IF (jscanspositively == 0) THEN
3965 lafirst = this%grid%grid%ymax
3966 lalast = this%grid%grid%ymin
3967 ELSE
3968 lafirst = this%grid%grid%ymin
3969 lalast = this%grid%grid%ymax
3970 ENDIF
3971
3972 CALL grib_set(gaid,'eastingOfFirstGridPoint',lofirst)
3973 CALL grib_set(gaid,'eastingOfLastGridPoint',lolast)
3974 CALL grib_set(gaid,'northingOfFirstGridPoint',lafirst)
3975 CALL grib_set(gaid,'northingOfLastGridPoint',lalast)
3976
3977CASE default
3979 "griddim_export_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
3980 CALL raise_error()
3981
3982END SELECT
3983
3984! hack for position of vertical coordinate parameters
3985! buggy in grib_api
3986IF (editionnumber == 1) THEN
3987! CALL grib_get(gaid,"PVPresent",pvp) ! alias, probably useless
3988 CALL grib_get(gaid,"NV",nv)
3989#ifdef DEBUG
3991 trim(to_char(nv))//" vertical coordinate parameters")
3992#endif
3993
3994 IF (nv == 0) THEN
3995 pvl = 255
3996 ELSE
3997 SELECT CASE (this%grid%proj%proj_type)
3998 CASE ('regular_ll') ! check whether "29-32 Set to zero (reserved)" are required
3999 pvl = 33
4000 CASE ('polar_stereographic')
4001 pvl = 33
4002 CASE ('rotated_ll', 'stretched_ll', 'lambert', 'albers')
4003 pvl = 43
4004 CASE ('stretched_rotated_ll')
4005 pvl = 43
4006 CASE DEFAULT
4007 pvl = 43 !?
4008 END SELECT
4009 ENDIF
4010
4011 CALL grib_set(gaid,"pvlLocation",pvl)
4012#ifdef DEBUG
4014 trim(to_char(pvl))//" as vertical coordinate parameter location")
4015#endif
4016
4017ENDIF
4018
4019
4020CONTAINS
4021! utilities routines for grib_api, is there a better place?
4022SUBROUTINE grib_set_dmiss(gaid, key, val, default, factor)
4023INTEGER,INTENT(in) :: gaid
4024CHARACTER(len=*),INTENT(in) :: key
4025DOUBLE PRECISION,INTENT(in) :: val
4026DOUBLE PRECISION,INTENT(in),OPTIONAL :: default
4027DOUBLE PRECISION,INTENT(in),OPTIONAL :: factor
4028
4029INTEGER :: ierr
4030
4032 IF (PRESENT(factor)) THEN
4033 CALL grib_set(gaid, key, val*factor, ierr)
4034 ELSE
4035 CALL grib_set(gaid, key, val, ierr)
4036 ENDIF
4037ELSE IF (PRESENT(default)) THEN
4038 CALL grib_set(gaid, key, default, ierr)
4039ENDIF
4040
4041END SUBROUTINE grib_set_dmiss
4042
4043SUBROUTINE grib_set_imiss(gaid, key, value, default)
4044INTEGER,INTENT(in) :: gaid
4045CHARACTER(len=*),INTENT(in) :: key
4046INTEGER,INTENT(in) :: value
4047INTEGER,INTENT(in),OPTIONAL :: default
4048
4049INTEGER :: ierr
4050
4052 CALL grib_set(gaid, key, value, ierr)
4053ELSE IF (PRESENT(default)) THEN
4054 CALL grib_set(gaid, key, default, ierr)
4055ENDIF
4056
4057END SUBROUTINE grib_set_imiss
4058
4059SUBROUTINE griddim_export_ellipsoid(this, gaid)
4060TYPE(griddim_def),INTENT(in) :: this
4061INTEGER,INTENT(in) :: gaid
4062
4063INTEGER :: ellips_type, ierr
4064DOUBLE PRECISION :: r1, r2, f
4065
4067
4068IF (editionnumber == 2) THEN
4069
4070! why does it produce a message even with ierr?
4071 CALL grib_set_missing(gaid, 'scaleFactorOfRadiusOfSphericalEarth', ierr)
4072 CALL grib_set_missing(gaid, 'scaledValueOfRadiusOfSphericalEarth', ierr)
4073 CALL grib_set_missing(gaid, 'scaleFactorOfEarthMajorAxis', ierr)
4074 CALL grib_set_missing(gaid, 'scaledValueOfEarthMajorAxis', ierr)
4075 CALL grib_set_missing(gaid, 'scaleFactorOfEarthMinorAxis', ierr)
4076 CALL grib_set_missing(gaid, 'scaledValueOfEarthMinorAxis', ierr)
4077
4078 SELECT CASE(ellips_type)
4079 CASE(ellips_grs80) ! iag-grs80
4080 CALL grib_set(gaid, 'shapeOfTheEarth', 4)
4081 CASE(ellips_wgs84) ! wgs84
4082 CALL grib_set(gaid, 'shapeOfTheEarth', 5)
4083 CASE default
4084 IF (f == 0.0d0) THEN ! spherical Earth
4085 IF (r1 == 6367470.0d0) THEN ! spherical
4086 CALL grib_set(gaid, 'shapeOfTheEarth', 0)
4087 ELSE IF (r1 == 6371229.0d0) THEN ! spherical
4088 CALL grib_set(gaid, 'shapeOfTheEarth', 6)
4089 ELSE ! spherical generic
4090 CALL grib_set(gaid, 'shapeOfTheEarth', 1)
4091 CALL grib_set(gaid, 'scaleFactorOfRadiusOfSphericalEarth', 2)
4092 CALL grib_set(gaid, 'scaledValueOfRadiusOfSphericalEarth', int(r1*100.0d0))
4093 ENDIF
4094 ELSE ! ellipsoidal
4095 IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
4096 CALL grib_set(gaid, 'shapeOfTheEarth', 2)
4097 ELSE ! ellipsoidal generic
4098 CALL grib_set(gaid, 'shapeOfTheEarth', 3)
4099 r2 = r1*(1.0d0 - f)
4100 CALL grib_set(gaid, 'scaleFactorOfEarthMajorAxis', 5)
4101 CALL grib_set(gaid, 'scaledValueOfEarthMajorAxis', &
4102 int(r1*100.0d0))
4103 CALL grib_set(gaid, 'scaleFactorOfEarthMinorAxis', 5)
4104 CALL grib_set(gaid, 'scaledValueOfEarthMinorAxis', &
4105 int(r2*100.0d0))
4106 ENDIF
4107 ENDIF
4108 END SELECT
4109
4110ELSE
4111
4112 IF (r1 == 6367470.0d0 .AND. f == 0.0d0) THEN ! spherical
4113 CALL grib_set(gaid, 'earthIsOblate', 0)
4114 ELSE IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
4115 CALL grib_set(gaid, 'earthIsOblate', 1)
4116 ELSE IF (f == 0.0d0) THEN ! generic spherical
4117 CALL grib_set(gaid, 'earthIsOblate', 0)
4119 ¬ supported in grib 1, coding with default radius of 6367470 m')
4120 ELSE ! generic ellipsoidal
4121 CALL grib_set(gaid, 'earthIsOblate', 1)
4123 ¬ supported in grib 1, coding with default iau65 ellipsoid')
4124 ENDIF
4125
4126ENDIF
4127
4128END SUBROUTINE griddim_export_ellipsoid
4129
4130SUBROUTINE get_unproj_first(this, iScansNegatively, jScansPositively, &
4131 loFirst, laFirst)
4132TYPE(griddim_def),INTENT(in) :: this ! griddim object
4133INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
4134DOUBLE PRECISION,INTENT(out) :: loFirst, laFirst
4135
4136! compute lon and lat of first point from projected extremes
4137IF (iscansnegatively == 0) THEN
4138 IF (jscanspositively == 0) THEN
4140 ELSE
4142 ENDIF
4143ELSE
4144 IF (jscanspositively == 0) THEN
4146 ELSE
4148 ENDIF
4149ENDIF
4150! use the values kept for personal pleasure ?
4151! loFirst = this%grid%proj%polar%lon1
4152! laFirst = this%grid%proj%polar%lat1
4153END SUBROUTINE get_unproj_first
4154
4155SUBROUTINE get_unproj_last(this, iScansNegatively, jScansPositively, &
4156 loLast, laLast)
4157TYPE(griddim_def),INTENT(in) :: this ! griddim object
4158INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
4159DOUBLE PRECISION,INTENT(out) :: loLast, laLast
4160
4161! compute lon and lat of last point from projected extremes
4162IF (iscansnegatively == 0) THEN
4163 IF (jscanspositively == 0) THEN
4165 ELSE
4167 ENDIF
4168ELSE
4169 IF (jscanspositively == 0) THEN
4171 ELSE
4173 ENDIF
4174ENDIF
4175! use the values kept for personal pleasure ?
4176! loLast = this%grid%proj%polar%lon?
4177! laLast = this%grid%proj%polar%lat?
4178END SUBROUTINE get_unproj_last
4179
4180END SUBROUTINE griddim_export_gribapi
4181#endif
4182
4183
4184#ifdef HAVE_LIBGDAL
4185! gdal driver
4186SUBROUTINE griddim_import_gdal(this, gdalid, gdal_options)
4187USE gdal
4188TYPE(griddim_def),INTENT(inout) :: this ! griddim object
4189TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal rasterband pointer
4190TYPE(gdal_file_id_options),INTENT(in) :: gdal_options
4191
4192TYPE(gdaldataseth) :: hds
4193REAL(kind=c_double) :: geotrans(6) !, invgeotrans(6), x1, y1, x2, y2
4194INTEGER(kind=c_int) :: offsetx, offsety
4195INTEGER :: ier
4196
4197hds = gdalgetbanddataset(gdalid) ! go back to dataset
4198ier = gdalgetgeotransform(hds, geotrans)
4199
4200IF (ier /= 0) THEN
4202 'griddim_import_gdal: error in accessing gdal dataset')
4203 CALL raise_error()
4204 RETURN
4205ENDIF
4206IF (geotrans(3) /= 0.0_c_double .OR. geotrans(5) /= 0.0_c_double) THEN ! dataset has not a diagonal transformation
4208 'griddim_import_gdal: dataset has a non-diagonal transformation matrix, unsupported')
4209 CALL raise_error()
4210 RETURN
4211ENDIF
4212
4213CALL gdaldatasetbbsize_f(hds, gdal_options%xmin, gdal_options%ymin, &
4214 gdal_options%xmax, gdal_options%ymax, &
4215 this%dim%nx, this%dim%ny, offsetx, offsety, &
4216 this%grid%grid%xmin, this%grid%grid%ymin, this%grid%grid%xmax, this%grid%grid%ymax)
4217
4218IF (this%dim%nx == 0 .OR. this%dim%ny == 0) THEN
4220 'griddim_import_gdal: requested bounding box '//t2c(gdal_options%xmin)//','// &
4221 t2c(gdal_options%ymin)//','//t2c(gdal_options%xmax)//','//&
4222 t2c(gdal_options%ymax))
4224 'determines an empty dataset '//t2c(this%dim%nx)//'x'//t2c(this%dim%ny))
4225ENDIF
4226
4227! get grid corners
4228!CALL gdalapplygeotransform(geotrans, 0.5_c_double, 0.5_c_double, x1, y1)
4229!CALL gdalapplygeotransform(geotrans, gdalgetrasterbandxsize(gdalid)-0.5_c_double, &
4230! gdalgetrasterbandysize(gdalid)-0.5_c_double, x2, y2)
4231
4232!IF (geotrans(3) == 0.0_c_double .AND. geotrans(5) == 0.0_c_double) THEN ! transformation is diagonal, no transposing
4233! this%dim%nx = gdalgetrasterbandxsize(gdalid)
4234! this%dim%ny = gdalgetrasterbandysize(gdalid)
4235! this%grid%grid%xmin = MIN(x1, x2)
4236! this%grid%grid%xmax = MAX(x1, x2)
4237! this%grid%grid%ymin = MIN(y1, y2)
4238! this%grid%grid%ymax = MAX(y1, y2)
4239!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
4240!
4241! this%dim%nx = gdalgetrasterbandysize(gdalid)
4242! this%dim%ny = gdalgetrasterbandxsize(gdalid)
4243! this%grid%grid%xmin = MIN(y1, y2)
4244! this%grid%grid%xmax = MAX(y1, y2)
4245! this%grid%grid%ymin = MIN(x1, x2)
4246! this%grid%grid%ymax = MAX(x1, x2)
4247!
4248!ELSE ! transformation is a rotation, not supported
4249!ENDIF
4250
4251this%grid%proj%proj_type = 'regular_ll' ! forced, only one supported (check?)
4252
4253CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
4254this%grid%grid%component_flag = 0
4255
4256END SUBROUTINE griddim_import_gdal
4257#endif
4258
4259
4261SUBROUTINE griddim_display(this)
4262TYPE(griddim_def),INTENT(in) :: this
4263
4264print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
4265
4269
4270print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
4271
4272END SUBROUTINE griddim_display
4273
4274
4275#define VOL7D_POLY_TYPE TYPE(grid_def)
4276#define VOL7D_POLY_TYPES _grid
4277#include "array_utilities_inc.F90"
4278#undef VOL7D_POLY_TYPE
4279#undef VOL7D_POLY_TYPES
4280
4281#define VOL7D_POLY_TYPE TYPE(griddim_def)
4282#define VOL7D_POLY_TYPES _griddim
4283#include "array_utilities_inc.F90"
4284#undef VOL7D_POLY_TYPE
4285#undef VOL7D_POLY_TYPES
4286
4287
4299SUBROUTINE griddim_wind_unrot(this, rot_mat)
4301TYPE(griddim_def), INTENT(in) :: this
4302DOUBLE PRECISION, POINTER :: rot_mat(:,:,:)
4303
4304DOUBLE PRECISION :: dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j
4305DOUBLE PRECISION :: lat_factor
4306INTEGER :: i, j, ip1, im1, jp1, jm1;
4307
4308IF (this%dim%nx <= 0 .OR. this%dim%ny <= 0 .OR. &
4309 .NOT. ASSOCIATED(this%dim%lon) .OR. .NOT. ASSOCIATED(this%dim%lat)) THEN
4310 CALL l4f_category_log(this%category, l4f_error, 'rotate_uv must be called after correct init of griddim object')
4311 NULLIFY(rot_mat)
4312 RETURN
4313ENDIF
4314
4315ALLOCATE(rot_mat(this%dim%nx, this%dim%ny, 4))
4316
4317DO j = 1, this%dim%ny
4318 jp1 = min(j+1, this%dim%ny)
4319 jm1 = max(j-1, 1)
4320 DO i = 1, this%dim%nx
4321 ip1 = min(i+1, this%dim%nx)
4322 im1 = max(i-1, 1)
4323
4324 dlat_i = this%dim%lat(ip1,j) - this%dim%lat(im1,j)
4325 dlat_j = this%dim%lat(i,jp1) - this%dim%lat(i,jm1)
4326
4327 dlon_i = this%dim%lon(ip1,j) - this%dim%lon(im1,j)
4328! if ( dlon_i > pi ) dlon_i -= 2*pi;
4329! if ( dlon_i < (-pi) ) dlon_i += 2*pi;
4330 dlon_j = this%dim%lon(i,jp1) - this%dim%lon(i,jm1)
4331! if ( dlon_j > pi ) dlon_j -= 2*pi;
4332! if ( dlon_j < (-pi) ) dlon_j += 2*pi;
4333
4334! check whether this is really necessary !!!!
4335 lat_factor = cos(degrad*this%dim%lat(i,j))
4336 dlon_i = dlon_i * lat_factor
4337 dlon_j = dlon_j * lat_factor
4338
4339 dist_i = sqrt(dlon_i*dlon_i + dlat_i*dlat_i)
4340 dist_j = sqrt(dlon_j*dlon_j + dlat_j*dlat_j)
4341
4342 IF (dist_i > 0.d0) THEN
4343 rot_mat(i,j,1) = dlon_i/dist_i
4344 rot_mat(i,j,2) = dlat_i/dist_i
4345 ELSE
4346 rot_mat(i,j,1) = 0.0d0
4347 rot_mat(i,j,2) = 0.0d0
4348 ENDIF
4349 IF (dist_j > 0.d0) THEN
4350 rot_mat(i,j,3) = dlon_j/dist_j
4351 rot_mat(i,j,4) = dlat_j/dist_j
4352 ELSE
4353 rot_mat(i,j,3) = 0.0d0
4354 rot_mat(i,j,4) = 0.0d0
4355 ENDIF
4356
4357 ENDDO
4358ENDDO
4359
4360END SUBROUTINE griddim_wind_unrot
4361
4362
4363! compute zoom indices from geographical zoom coordinates
4364SUBROUTINE griddim_zoom_coord(this, ilon, ilat, flon, flat, ix, iy, fx, fy)
4365TYPE(griddim_def),INTENT(in) :: this
4366DOUBLE PRECISION,INTENT(in) :: ilon, ilat, flon, flat
4367INTEGER,INTENT(out) :: ix, iy, fx, fy
4368
4369DOUBLE PRECISION :: ix1, iy1, fx1, fy1
4370
4371! compute projected coordinates of vertices of desired lonlat rectangle
4374
4375CALL griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
4376
4377END SUBROUTINE griddim_zoom_coord
4378
4379
4380! compute zoom indices from projected zoom coordinates
4381SUBROUTINE griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
4382TYPE(griddim_def),INTENT(in) :: this
4383DOUBLE PRECISION,INTENT(in) :: ix1, iy1, fx1, fy1
4384INTEGER,INTENT(out) :: ix, iy, fx, fy
4385
4386INTEGER :: lix, liy, lfx, lfy
4387
4388! compute projected indices
4389lix = nint((ix1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
4390liy = nint((iy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
4391lfx = nint((fx1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
4392lfy = nint((fy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
4393! swap projected indices, in case grid is "strongly rotated"
4394ix = min(lix, lfx)
4395fx = max(lix, lfx)
4396iy = min(liy, lfy)
4397fy = max(liy, lfy)
4398
4399END SUBROUTINE griddim_zoom_projcoord
4400
4401
4405SUBROUTINE long_reset_0_360(lon)
4406DOUBLE PRECISION,INTENT(inout) :: lon
4407
4409DO WHILE(lon < 0.0d0)
4410 lon = lon + 360.0d0
4411END DO
4412DO WHILE(lon >= 360.0d0)
4413 lon = lon - 360.0d0
4414END DO
4415
4416END SUBROUTINE long_reset_0_360
4417
4418
4422SUBROUTINE long_reset_m180_360(lon)
4423DOUBLE PRECISION,INTENT(inout) :: lon
4424
4426DO WHILE(lon < -180.0d0)
4427 lon = lon + 360.0d0
4428END DO
4429DO WHILE(lon >= 360.0d0)
4430 lon = lon - 360.0d0
4431END DO
4432
4433END SUBROUTINE long_reset_m180_360
4434
4435
4439!SUBROUTINE long_reset_m90_270(lon)
4440!DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
4441!
4442!IF (.NOT.c_e(lon)) RETURN
4443!DO WHILE(lon < -90.0D0)
4444! lon = lon + 360.0D0
4445!END DO
4446!DO WHILE(lon >= 270.0D0)
4447! lon = lon - 360.0D0
4448!END DO
4449!
4450!END SUBROUTINE long_reset_m90_270
4451
4452
4456SUBROUTINE long_reset_m180_180(lon)
4457DOUBLE PRECISION,INTENT(inout) :: lon
4458
4460DO WHILE(lon < -180.0d0)
4461 lon = lon + 360.0d0
4462END DO
4463DO WHILE(lon >= 180.0d0)
4464 lon = lon - 360.0d0
4465END DO
4466
4467END SUBROUTINE long_reset_m180_180
4468
4469
4470SUBROUTINE long_reset_to_cart_closest(lon, lonref)
4471DOUBLE PRECISION,INTENT(inout) :: lon
4472DOUBLE PRECISION,INTENT(in) :: lonref
4473
4475IF (abs(lon-lonref) < 180.0d0) RETURN ! nothing to do
4476lon = lon - nint((lon-lonref)/360.0d0)*360.0d0 ! se non e` vera e` ben trovata
4477
4478END SUBROUTINE long_reset_to_cart_closest
4479
4480
4482
Method for testing the existence of the object. Definition: geo_proj_class.F90:268 Compute forward coordinate transformation from geographical system to projected system. Definition: grid_class.F90:308 Read the object from a formatted or unformatted file. Definition: grid_class.F90:334 Compute backward coordinate transformation from projected system to geographical system. Definition: grid_class.F90:314 Write the object on a formatted or unformatted file. Definition: grid_class.F90:329 Emit log message for a category with specific priority. Definition: log4fortran.F90:457 Module for describing geographically referenced regular grids. Definition: grid_class.F90:237 Module for defining the extension and coordinates of a rectangular georeferenced grid. Definition: grid_dim_class.F90:211 This module defines an abstract interface to different drivers for access to files containing gridded... Definition: grid_id_class.F90:249 Definitions of constants and functions for working with missing values. Definition: missing_values.f90:50 Module for quickly interpreting the OPTIONAL parameters passed to a subprogram. Definition: optional_values.f90:28 This object, mainly for internal use, describes a grid on a geographical projection,... Definition: grid_class.F90:257 This object completely describes a grid on a geographic projection. Definition: grid_class.F90:270 Derived type describing the extension of a grid and the geographical coordinates of each point. Definition: grid_dim_class.F90:221 |