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