libsim Versione 7.2.1
|
◆ map_distinct_griddim()
map distinct Definizione alla linea 2721 del file grid_class.F90. 2722! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2723! authors:
2724! Davide Cesari <dcesari@arpa.emr.it>
2725! Paolo Patruno <ppatruno@arpa.emr.it>
2726
2727! This program is free software; you can redistribute it and/or
2728! modify it under the terms of the GNU General Public License as
2729! published by the Free Software Foundation; either version 2 of
2730! the License, or (at your option) any later version.
2731
2732! This program is distributed in the hope that it will be useful,
2733! but WITHOUT ANY WARRANTY; without even the implied warranty of
2734! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2735! GNU General Public License for more details.
2736
2737! You should have received a copy of the GNU General Public License
2738! along with this program. If not, see <http://www.gnu.org/licenses/>.
2739#include "config.h"
2740
2771use geo_proj_class
2773use grid_rect_class
2779implicit none
2780
2781
2782character (len=255),parameter:: subcategory="grid_class"
2783
2784
2791 private
2792 type(geo_proj) :: proj
2793 type(grid_rect) :: grid
2794 integer :: category = 0
2796
2797
2804 type(grid_def) :: grid
2805 type(grid_dim) :: dim
2806 integer :: category = 0
2808
2809
2813INTERFACE OPERATOR (==)
2814 MODULE PROCEDURE grid_eq, griddim_eq
2815END INTERFACE
2816
2820INTERFACE OPERATOR (/=)
2821 MODULE PROCEDURE grid_ne, griddim_ne
2822END INTERFACE
2823
2826 MODULE PROCEDURE griddim_init
2827END INTERFACE
2828
2831 MODULE PROCEDURE griddim_delete
2832END INTERFACE
2833
2836 MODULE PROCEDURE griddim_copy
2837END INTERFACE
2838
2842 MODULE PROCEDURE griddim_coord_proj!, griddim_proj
2843END INTERFACE
2844
2848 MODULE PROCEDURE griddim_coord_unproj, griddim_unproj
2849END INTERFACE
2850
2853 MODULE PROCEDURE griddim_get_val
2854END INTERFACE
2855
2858 MODULE PROCEDURE griddim_set_val
2859END INTERFACE
2860
2863 MODULE PROCEDURE griddim_write_unit
2864END INTERFACE
2865
2868 MODULE PROCEDURE griddim_read_unit
2869END INTERFACE
2870
2872INTERFACE import
2873 MODULE PROCEDURE griddim_import_grid_id
2874END INTERFACE
2875
2878 MODULE PROCEDURE griddim_export_grid_id
2879END INTERFACE
2880
2883 MODULE PROCEDURE griddim_display
2884END INTERFACE
2885
2886#define VOL7D_POLY_TYPE TYPE(grid_def)
2887#define VOL7D_POLY_TYPES _grid
2888#include "array_utilities_pre.F90"
2889#undef VOL7D_POLY_TYPE
2890#undef VOL7D_POLY_TYPES
2891
2892#define VOL7D_POLY_TYPE TYPE(griddim_def)
2893#define VOL7D_POLY_TYPES _griddim
2894#include "array_utilities_pre.F90"
2895#undef VOL7D_POLY_TYPE
2896#undef VOL7D_POLY_TYPES
2897
2898INTERFACE wind_unrot
2899 MODULE PROCEDURE griddim_wind_unrot
2900END INTERFACE
2901
2902
2903PRIVATE
2904
2906 griddim_zoom_coord, griddim_zoom_projcoord, &
2910PUBLIC OPERATOR(==),OPERATOR(/=)
2911PUBLIC count_distinct, pack_distinct, count_and_pack_distinct, &
2912 map_distinct, map_inv_distinct,index
2914PUBLIC griddim_central_lon, griddim_set_central_lon
2915CONTAINS
2916
2918SUBROUTINE griddim_init(this, nx, ny, &
2919 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
2920 proj_type, lov, zone, xoff, yoff, &
2921 longitude_south_pole, latitude_south_pole, angle_rotation, &
2922 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
2923 latin1, latin2, lad, projection_center_flag, &
2924 ellips_smaj_axis, ellips_flatt, ellips_type, &
2925 categoryappend)
2926TYPE(griddim_def),INTENT(inout) :: this
2927INTEGER,INTENT(in),OPTIONAL :: nx
2928INTEGER,INTENT(in),OPTIONAL :: ny
2929DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin
2930DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmax
2931DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymin
2932DOUBLE PRECISION,INTENT(in),OPTIONAL :: ymax
2933DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx
2934DOUBLE PRECISION,INTENT(in),OPTIONAL :: dy
2937INTEGER,INTENT(in),OPTIONAL :: component_flag
2938CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
2939DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
2940INTEGER,INTENT(in),OPTIONAL :: zone
2941DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
2942DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
2943DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
2944DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
2945DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
2946DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
2947DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
2948DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
2949DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
2950DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
2951DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
2952INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
2953DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
2954DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
2955INTEGER,INTENT(in),OPTIONAL :: ellips_type
2956CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2957
2958CHARACTER(len=512) :: a_name
2959
2960IF (PRESENT(categoryappend)) THEN
2961 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
2962 trim(categoryappend))
2963ELSE
2964 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
2965ENDIF
2966this%category=l4f_category_get(a_name)
2967
2968! geographical projection
2969this%grid%proj = geo_proj_new( &
2970 proj_type=proj_type, lov=lov, zone=zone, xoff=xoff, yoff=yoff, &
2971 longitude_south_pole=longitude_south_pole, &
2972 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
2973 longitude_stretch_pole=longitude_stretch_pole, &
2974 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
2975 latin1=latin1, latin2=latin2, lad=lad, projection_center_flag=projection_center_flag, &
2976 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, ellips_type=ellips_type)
2977! grid extension
2978this%grid%grid = grid_rect_new( &
2979 xmin, xmax, ymin, ymax, dx, dy, component_flag)
2980! grid size
2981this%dim = grid_dim_new(nx, ny)
2982
2983#ifdef DEBUG
2985#endif
2986
2987END SUBROUTINE griddim_init
2988
2989
2991SUBROUTINE griddim_delete(this)
2992TYPE(griddim_def),INTENT(inout) :: this
2993
2997
2998call l4f_category_delete(this%category)
2999
3000END SUBROUTINE griddim_delete
3001
3002
3004SUBROUTINE griddim_copy(this, that, categoryappend)
3005TYPE(griddim_def),INTENT(in) :: this
3006TYPE(griddim_def),INTENT(out) :: that
3007CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
3008
3009CHARACTER(len=512) :: a_name
3010
3012
3016
3017! new category
3018IF (PRESENT(categoryappend)) THEN
3019 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."// &
3020 trim(categoryappend))
3021ELSE
3022 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
3023ENDIF
3024that%category=l4f_category_get(a_name)
3025
3026END SUBROUTINE griddim_copy
3027
3028
3031ELEMENTAL SUBROUTINE griddim_coord_proj(this, lon, lat, x, y)
3032TYPE(griddim_def),INTENT(in) :: this
3034DOUBLE PRECISION,INTENT(in) :: lon, lat
3036DOUBLE PRECISION,INTENT(out) :: x, y
3037
3039
3040END SUBROUTINE griddim_coord_proj
3041
3042
3045ELEMENTAL SUBROUTINE griddim_coord_unproj(this, x, y, lon, lat)
3046TYPE(griddim_def),INTENT(in) :: this
3048DOUBLE PRECISION,INTENT(in) :: x, y
3050DOUBLE PRECISION,INTENT(out) :: lon, lat
3051
3053
3054END SUBROUTINE griddim_coord_unproj
3055
3056
3057! Computes and sets the grid parameters required to compute
3058! coordinates of grid points in the projected system.
3059! probably meaningless
3060!SUBROUTINE griddim_proj(this)
3061!TYPE(griddim_def),INTENT(inout) :: this !< definition of projection and grid
3062!
3063!CALL proj(this, this%dim%lon(1,1), this%dim%lat(1,1), &
3064! this%grid%grid%xmin, this%grid%grid%ymin)
3065!
3066!CALL proj(this, this%dim%lon(this%dim%nx,this%dim%ny), &
3067! this%dim%lat(this%dim%nx,this%dim%ny), &
3068! this%grid%grid%xmax, this%grid%grid%ymax)
3069!
3070!END SUBROUTINE griddim_proj
3071
3079SUBROUTINE griddim_unproj(this)
3080TYPE(griddim_def),INTENT(inout) :: this
3081
3083CALL alloc(this%dim)
3084CALL griddim_unproj_internal(this)
3085
3086END SUBROUTINE griddim_unproj
3087
3088! internal subroutine needed for allocating automatic arrays
3089SUBROUTINE griddim_unproj_internal(this)
3090TYPE(griddim_def),INTENT(inout) ::this ! definition of projection and grid
3091
3092DOUBLE PRECISION :: x(this%dim%nx,this%dim%ny), y(this%dim%nx,this%dim%ny)
3093
3094CALL grid_rect_coordinates(this%grid%grid, x, y)
3096
3097END SUBROUTINE griddim_unproj_internal
3098
3099
3101SUBROUTINE griddim_get_val(this, nx, ny, &
3102 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
3103 proj, proj_type, lov, zone, xoff, yoff, &
3104 longitude_south_pole, latitude_south_pole, angle_rotation, &
3105 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
3106 latin1, latin2, lad, projection_center_flag, &
3107 ellips_smaj_axis, ellips_flatt, ellips_type)
3108TYPE(griddim_def),INTENT(in) :: this
3109INTEGER,INTENT(out),OPTIONAL :: nx
3110INTEGER,INTENT(out),OPTIONAL :: ny
3112DOUBLE PRECISION,INTENT(out),OPTIONAL :: xmin, xmax, ymin, ymax
3114DOUBLE PRECISION,INTENT(out),OPTIONAL :: dx, dy
3117INTEGER,INTENT(out),OPTIONAL :: component_flag
3118TYPE(geo_proj),INTENT(out),OPTIONAL :: proj
3119CHARACTER(len=*),INTENT(out),OPTIONAL :: proj_type
3120DOUBLE PRECISION,INTENT(out),OPTIONAL :: lov
3121INTEGER,INTENT(out),OPTIONAL :: zone
3122DOUBLE PRECISION,INTENT(out),OPTIONAL :: xoff
3123DOUBLE PRECISION,INTENT(out),OPTIONAL :: yoff
3124DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_south_pole
3125DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_south_pole
3126DOUBLE PRECISION,INTENT(out),OPTIONAL :: angle_rotation
3127DOUBLE PRECISION,INTENT(out),OPTIONAL :: longitude_stretch_pole
3128DOUBLE PRECISION,INTENT(out),OPTIONAL :: latitude_stretch_pole
3129DOUBLE PRECISION,INTENT(out),OPTIONAL :: stretch_factor
3130DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin1
3131DOUBLE PRECISION,INTENT(out),OPTIONAL :: latin2
3132DOUBLE PRECISION,INTENT(out),OPTIONAL :: lad
3133INTEGER,INTENT(out),OPTIONAL :: projection_center_flag
3134DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_smaj_axis
3135DOUBLE PRECISION,INTENT(out),OPTIONAL :: ellips_flatt
3136INTEGER,INTENT(out),OPTIONAL :: ellips_type
3137
3138IF (PRESENT(nx)) nx = this%dim%nx
3139IF (PRESENT(ny)) ny = this%dim%ny
3140
3142
3144 xoff=xoff, yoff=yoff, &
3145 longitude_south_pole=longitude_south_pole, &
3146 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
3147 longitude_stretch_pole=longitude_stretch_pole, &
3148 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
3149 latin1=latin1, latin2=latin2, lad=lad, &
3150 projection_center_flag=projection_center_flag, &
3151 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
3152 ellips_type=ellips_type)
3153
3155 xmin, xmax, ymin, ymax, dx, dy, component_flag)
3156
3157END SUBROUTINE griddim_get_val
3158
3159
3161SUBROUTINE griddim_set_val(this, nx, ny, &
3162 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
3163 proj_type, lov, zone, xoff, yoff, &
3164 longitude_south_pole, latitude_south_pole, angle_rotation, &
3165 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
3166 latin1, latin2, lad, projection_center_flag, &
3167 ellips_smaj_axis, ellips_flatt, ellips_type)
3168TYPE(griddim_def),INTENT(inout) :: this
3169INTEGER,INTENT(in),OPTIONAL :: nx
3170INTEGER,INTENT(in),OPTIONAL :: ny
3172DOUBLE PRECISION,INTENT(in),OPTIONAL :: xmin, xmax, ymin, ymax
3174DOUBLE PRECISION,INTENT(in),OPTIONAL :: dx, dy
3177INTEGER,INTENT(in),OPTIONAL :: component_flag
3178CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
3179DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
3180INTEGER,INTENT(in),OPTIONAL :: zone
3181DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
3182DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
3183DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
3184DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
3185DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
3186DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
3187DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
3188DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
3189DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
3190DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
3191DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
3192INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
3193DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
3194DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
3195INTEGER,INTENT(in),OPTIONAL :: ellips_type
3196
3197IF (PRESENT(nx)) this%dim%nx = nx
3198IF (PRESENT(ny)) this%dim%ny = ny
3199
3201 xoff=xoff, yoff=yoff, longitude_south_pole=longitude_south_pole, &
3202 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
3203 longitude_stretch_pole=longitude_stretch_pole, &
3204 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
3205 latin1=latin1, latin2=latin2, lad=lad, &
3206 projection_center_flag=projection_center_flag, &
3207 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
3208 ellips_type=ellips_type)
3209
3211 xmin, xmax, ymin, ymax, dx, dy, component_flag)
3212
3213END SUBROUTINE griddim_set_val
3214
3215
3220SUBROUTINE griddim_read_unit(this, unit)
3221TYPE(griddim_def),INTENT(out) :: this
3222INTEGER, INTENT(in) :: unit
3223
3224
3228
3229END SUBROUTINE griddim_read_unit
3230
3231
3236SUBROUTINE griddim_write_unit(this, unit)
3237TYPE(griddim_def),INTENT(in) :: this
3238INTEGER, INTENT(in) :: unit
3239
3240
3244
3245END SUBROUTINE griddim_write_unit
3246
3247
3251FUNCTION griddim_central_lon(this) RESULT(lon)
3252TYPE(griddim_def),INTENT(inout) :: this
3253
3254DOUBLE PRECISION :: lon
3255
3256CALL griddim_pistola_central_lon(this, lon)
3257
3258END FUNCTION griddim_central_lon
3259
3260
3264SUBROUTINE griddim_set_central_lon(this, lonref)
3265TYPE(griddim_def),INTENT(inout) :: this
3266DOUBLE PRECISION,INTENT(in) :: lonref
3267
3268DOUBLE PRECISION :: lon
3269
3270CALL griddim_pistola_central_lon(this, lon, lonref)
3271
3272END SUBROUTINE griddim_set_central_lon
3273
3274
3275! internal subroutine for performing tasks common to the prevous two
3276SUBROUTINE griddim_pistola_central_lon(this, lon, lonref)
3277TYPE(griddim_def),INTENT(inout) :: this ! grid descriptor
3278DOUBLE PRECISION,INTENT(inout) :: lon ! central longitude
3279DOUBLE PRECISION,INTENT(in),OPTIONAL :: lonref ! reference longitude
3280
3281INTEGER :: unit
3282DOUBLE PRECISION :: lonsp, latsp, londelta, lov, lonrot
3283CHARACTER(len=80) :: ptype
3284
3285lon = dmiss
3287IF (unit == geo_proj_unit_meter) THEN ! it is a plane projection
3289 IF (PRESENT(lonref)) THEN
3290 CALL long_reset_to_cart_closest(lov, lonref)
3292 ENDIF
3293
3294ELSE IF (unit == geo_proj_unit_degree) THEN ! it is a spheric projection
3296 longitude_south_pole=lonsp, latitude_south_pole=latsp)
3297 SELECT CASE(ptype)
3298 CASE('rotated_ll','stretched_rotated_ll') ! use origin of rotated system
3299 IF (latsp < 0.0d0) THEN
3300 lon = lonsp
3301 IF (PRESENT(lonref)) THEN
3302 CALL long_reset_to_cart_closest(lon, lonref)
3304! now reset rotated coordinates around zero
3306 lonrot = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
3307 ENDIF
3308 londelta = lonrot
3309 CALL long_reset_to_cart_closest(londelta, 0.0d0)
3310 londelta = londelta - lonrot
3311 this%grid%grid%xmin = this%grid%grid%xmin + londelta
3312 this%grid%grid%xmax = this%grid%grid%xmax + londelta
3313 ENDIF
3314 ELSE ! this part to be checked
3315 lon = modulo(lonsp + 180.0d0, 360.0d0)
3316! IF (PRESENT(lonref)) THEN
3317! CALL long_reset_to_cart_closest(lov, lonref)
3318! CALL set_val(this%grid%proj, longitude_south_pole=lonref)
3319! ENDIF
3320 ENDIF
3321 CASE default ! use real grid limits
3323 lon = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
3324 ENDIF
3325 IF (PRESENT(lonref)) THEN
3326 londelta = lon
3327 CALL long_reset_to_cart_closest(londelta, lonref)
3328 londelta = londelta - lon
3329 this%grid%grid%xmin = this%grid%grid%xmin + londelta
3330 this%grid%grid%xmax = this%grid%grid%xmax + londelta
3331 ENDIF
3332 END SELECT
3333ENDIF
3334
3335END SUBROUTINE griddim_pistola_central_lon
3336
3337
3341SUBROUTINE griddim_gen_coord(this, x, y)
3342TYPE(griddim_def),INTENT(in) :: this
3343DOUBLE PRECISION,INTENT(out) :: x(:,:)
3344DOUBLE PRECISION,INTENT(out) :: y(:,:)
3345
3346
3347CALL grid_rect_coordinates(this%grid%grid, x, y)
3348
3349END SUBROUTINE griddim_gen_coord
3350
3351
3353SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
3354TYPE(griddim_def), INTENT(in) :: this
3355INTEGER,INTENT(in) :: nx
3356INTEGER,INTENT(in) :: ny
3357DOUBLE PRECISION,INTENT(out) :: dx
3358DOUBLE PRECISION,INTENT(out) :: dy
3359
3360CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
3361
3362END SUBROUTINE griddim_steps
3363
3364
3366SUBROUTINE griddim_setsteps(this)
3367TYPE(griddim_def), INTENT(inout) :: this
3368
3369CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3370
3371END SUBROUTINE griddim_setsteps
3372
3373
3374! TODO
3375! bisogna sviluppare gli altri operatori
3376ELEMENTAL FUNCTION grid_eq(this, that) RESULT(res)
3377TYPE(grid_def),INTENT(IN) :: this, that
3378LOGICAL :: res
3379
3380res = this%proj == that%proj .AND. &
3381 this%grid == that%grid
3382
3383END FUNCTION grid_eq
3384
3385
3386ELEMENTAL FUNCTION griddim_eq(this, that) RESULT(res)
3387TYPE(griddim_def),INTENT(IN) :: this, that
3388LOGICAL :: res
3389
3390res = this%grid == that%grid .AND. &
3391 this%dim == that%dim
3392
3393END FUNCTION griddim_eq
3394
3395
3396ELEMENTAL FUNCTION grid_ne(this, that) RESULT(res)
3397TYPE(grid_def),INTENT(IN) :: this, that
3398LOGICAL :: res
3399
3400res = .NOT.(this == that)
3401
3402END FUNCTION grid_ne
3403
3404
3405ELEMENTAL FUNCTION griddim_ne(this, that) RESULT(res)
3406TYPE(griddim_def),INTENT(IN) :: this, that
3407LOGICAL :: res
3408
3409res = .NOT.(this == that)
3410
3411END FUNCTION griddim_ne
3412
3413
3419SUBROUTINE griddim_import_grid_id(this, ingrid_id)
3420#ifdef HAVE_LIBGDAL
3421USE gdal
3422#endif
3423TYPE(griddim_def),INTENT(inout) :: this
3424TYPE(grid_id),INTENT(in) :: ingrid_id
3425
3426#ifdef HAVE_LIBGRIBAPI
3427INTEGER :: gaid
3428#endif
3429#ifdef HAVE_LIBGDAL
3430TYPE(gdalrasterbandh) :: gdalid
3431#endif
3433
3434#ifdef HAVE_LIBGRIBAPI
3435gaid = grid_id_get_gaid(ingrid_id)
3437#endif
3438#ifdef HAVE_LIBGDAL
3439gdalid = grid_id_get_gdalid(ingrid_id)
3440IF (gdalassociated(gdalid)) CALL griddim_import_gdal(this, gdalid, &
3441 grid_id_get_gdal_options(ingrid_id))
3442#endif
3443
3444END SUBROUTINE griddim_import_grid_id
3445
3446
3451SUBROUTINE griddim_export_grid_id(this, outgrid_id)
3452#ifdef HAVE_LIBGDAL
3453USE gdal
3454#endif
3455TYPE(griddim_def),INTENT(in) :: this
3456TYPE(grid_id),INTENT(inout) :: outgrid_id
3457
3458#ifdef HAVE_LIBGRIBAPI
3459INTEGER :: gaid
3460#endif
3461#ifdef HAVE_LIBGDAL
3462TYPE(gdalrasterbandh) :: gdalid
3463#endif
3464
3465#ifdef HAVE_LIBGRIBAPI
3466gaid = grid_id_get_gaid(outgrid_id)
3468#endif
3469#ifdef HAVE_LIBGDAL
3470gdalid = grid_id_get_gdalid(outgrid_id)
3471!IF (gdalassociated(gdalid)
3472! export for gdal not implemented, log?
3473#endif
3474
3475END SUBROUTINE griddim_export_grid_id
3476
3477
3478#ifdef HAVE_LIBGRIBAPI
3479! grib_api driver
3480SUBROUTINE griddim_import_gribapi(this, gaid)
3481USE grib_api
3482TYPE(griddim_def),INTENT(inout) :: this ! griddim object
3483INTEGER, INTENT(in) :: gaid ! grib_api id of the grib loaded in memory to import
3484
3485DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, x1, y1, orient
3486INTEGER :: EditionNumber, iScansNegatively, jScansPositively, zone, datum, &
3487 reflon, ierr
3488
3489! Generic keys
3490CALL grib_get(gaid, 'typeOfGrid', this%grid%proj%proj_type)
3491#ifdef DEBUG
3493 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type))
3494#endif
3495CALL grib_get(gaid,'GRIBEditionNumber',editionnumber)
3496
3497IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
3498 CALL grib_get(gaid,'numberOfDataPoints', this%dim%nx)
3499 this%dim%ny = 1
3500 this%grid%grid%component_flag = 0
3501 CALL griddim_import_ellipsoid(this, gaid)
3502 RETURN
3503ENDIF
3504
3505! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
3506CALL grib_get(gaid, 'Ni', this%dim%nx)
3507CALL grib_get(gaid, 'Nj', this%dim%ny)
3508CALL griddim_import_ellipsoid(this, gaid) ! placed here, not valid for utm datum /= 1
3509
3510CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
3511CALL grib_get(gaid,'jScansPositively',jscanspositively)
3512CALL grib_get(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
3513
3514! Keys for rotated grids (checked through missing values)
3515CALL grib_get_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
3516 this%grid%proj%rotated%longitude_south_pole)
3517CALL grib_get_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
3518 this%grid%proj%rotated%latitude_south_pole)
3519CALL grib_get_dmiss(gaid,'angleOfRotationInDegrees', &
3520 this%grid%proj%rotated%angle_rotation)
3521
3522! Keys for stretched grids (checked through missing values)
3523! units must be verified, still experimental in grib_api
3524! # TODO: Is it a float? Is it signed?
3525IF (editionnumber == 1) THEN
3526 CALL grib_get_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
3527 this%grid%proj%stretched%longitude_stretch_pole)
3528 CALL grib_get_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
3529 this%grid%proj%stretched%latitude_stretch_pole)
3530 CALL grib_get_dmiss(gaid,'stretchingFactor', &
3531 this%grid%proj%stretched%stretch_factor)
3532ELSE IF (editionnumber == 2) THEN
3533 CALL grib_get_dmiss(gaid,'longitudeOfThePoleOfStretching', &
3534 this%grid%proj%stretched%longitude_stretch_pole)
3535 CALL grib_get_dmiss(gaid,'latitudeOfThePoleOfStretching', &
3536 this%grid%proj%stretched%latitude_stretch_pole)
3537 CALL grib_get_dmiss(gaid,'stretchingFactor', &
3538 this%grid%proj%stretched%stretch_factor)
3540 this%grid%proj%stretched%stretch_factor = this%grid%proj%stretched%stretch_factor*1.0d-6
3541ENDIF
3542
3543! Projection-dependent keys
3544SELECT CASE (this%grid%proj%proj_type)
3545
3546! Keys for spherical coordinate systems
3547CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
3548
3549 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3550 CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3551 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3552 CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3553
3554! longitudes are sometimes wrongly coded even in grib2 and even by the
3555! Metoffice!
3556! longitudeOfFirstGridPointInDegrees = 354.911;
3557! longitudeOfLastGridPointInDegrees = 363.311;
3558 CALL long_reset_0_360(lofirst)
3559 CALL long_reset_0_360(lolast)
3560
3561 IF (iscansnegatively == 0) THEN
3562 this%grid%grid%xmin = lofirst
3563 this%grid%grid%xmax = lolast
3564 ELSE
3565 this%grid%grid%xmax = lofirst
3566 this%grid%grid%xmin = lolast
3567 ENDIF
3568 IF (jscanspositively == 0) THEN
3569 this%grid%grid%ymax = lafirst
3570 this%grid%grid%ymin = lalast
3571 ELSE
3572 this%grid%grid%ymin = lafirst
3573 this%grid%grid%ymax = lalast
3574 ENDIF
3575
3576! reset longitudes in order to have a Cartesian plane
3577 IF (this%grid%grid%xmax < this%grid%grid%xmin) &
3578 this%grid%grid%xmin = this%grid%grid%xmin - 360.d0
3579
3580! compute dx and dy (should we get them from grib?)
3581 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3582
3583! Keys for polar projections
3584CASE ('polar_stereographic', 'lambert', 'albers')
3585
3586 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
3587 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
3588! latin1/latin2 may be missing (e.g. stereographic)
3589 CALL grib_get_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
3590 CALL grib_get_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
3591#ifdef DEBUG
3593 "griddim_import_gribapi, latin1/2 "// &
3594 trim(to_char(this%grid%proj%polar%latin1))//" "// &
3595 trim(to_char(this%grid%proj%polar%latin2)))
3596#endif
3597! projection center flag, aka hemisphere
3598 CALL grib_get(gaid,'projectionCenterFlag',&
3599 this%grid%proj%polar%projection_center_flag, ierr)
3600 IF (ierr /= grib_success) THEN ! try center/centre
3601 CALL grib_get(gaid,'projectionCentreFlag',&
3602 this%grid%proj%polar%projection_center_flag)
3603 ENDIF
3604
3605 IF (iand(this%grid%proj%polar%projection_center_flag,64) == 1) THEN
3607 "griddim_import_gribapi, bi-polar projections not supported")
3608 CALL raise_error()
3609 ENDIF
3610! line of view, aka central meridian
3611 CALL grib_get(gaid,'LoVInDegrees',this%grid%proj%lov)
3612#ifdef DEBUG
3614 "griddim_import_gribapi, central meridian "//trim(to_char(this%grid%proj%lov)))
3615#endif
3616
3617! latitude at which dx and dy are valid
3618 IF (editionnumber == 1) THEN
3619! ECMWF (gribex/grib_api) says: Grid lengths are in metres, at the
3620! 60-degree parallel nearest to the pole on the projection plane.
3621! IF (IAND(this%projection_center_flag, 128) == 0) THEN
3622! this%grid%proj%polar%lad = 60.D0
3623! ELSE
3624! this%grid%proj%polar%lad = -60.D0
3625! ENDIF
3626! WMO says: Grid lengths are in units of metres, at the secant cone
3627! intersection parallel nearest to the pole on the projection plane.
3628 this%grid%proj%polar%lad = this%grid%proj%polar%latin1
3629 ELSE IF (editionnumber == 2) THEN
3630 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3631 ENDIF
3632#ifdef DEBUG
3634 "griddim_import_gribapi, lad "//trim(to_char(this%grid%proj%polar%lad)))
3635#endif
3636
3637! compute projected extremes from lon and lat of first point
3638 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3639 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3640 CALL long_reset_0_360(lofirst)
3641 CALL long_reset_to_cart_closest(this%grid%proj%lov, lofirst)
3642#ifdef DEBUG
3644 "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
3646 "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
3647#endif
3648
3650 IF (iscansnegatively == 0) THEN
3651 this%grid%grid%xmin = x1
3652 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
3653 ELSE
3654 this%grid%grid%xmax = x1
3655 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
3656 ENDIF
3657 IF (jscanspositively == 0) THEN
3658 this%grid%grid%ymax = y1
3659 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
3660 ELSE
3661 this%grid%grid%ymin = y1
3662 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
3663 ENDIF
3664! keep these values for personal pleasure
3665 this%grid%proj%polar%lon1 = lofirst
3666 this%grid%proj%polar%lat1 = lafirst
3667
3668! Keys for equatorial projections
3669CASE ('mercator')
3670
3671 CALL grib_get(gaid,'DxInMetres', this%grid%grid%dx)
3672 CALL grib_get(gaid,'DyInMetres', this%grid%grid%dy)
3673 CALL grib_get(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
3674 this%grid%proj%lov = 0.0d0 ! not in grib
3675
3676 CALL grib_get(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3677 CALL grib_get(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3678
3680 IF (iscansnegatively == 0) THEN
3681 this%grid%grid%xmin = x1
3682 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
3683 ELSE
3684 this%grid%grid%xmax = x1
3685 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
3686 ENDIF
3687 IF (jscanspositively == 0) THEN
3688 this%grid%grid%ymax = y1
3689 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
3690 ELSE
3691 this%grid%grid%ymin = y1
3692 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
3693 ENDIF
3694
3695 IF (editionnumber == 2) THEN
3696 CALL grib_get(gaid,'orientationOfTheGridInDegrees',orient)
3697 IF (orient /= 0.0d0) THEN
3699 "griddim_import_gribapi, Mercator grid orientation != 0 not supported")
3700 CALL raise_error()
3701 ENDIF
3702 ENDIF
3703
3704#ifdef DEBUG
3707 "griddim_import_gribapi, unprojected first point "//t2c(lofirst)//" "// &
3708 t2c(lafirst))
3709
3710 CALL grib_get(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3711 CALL grib_get(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3714 "griddim_import_gribapi, extremes from grib "//t2c(x1)//" "//t2c(y1))
3716 "griddim_import_gribapi, extremes from proj "//t2c(this%grid%grid%xmin)// &
3717 " "//t2c(this%grid%grid%ymin)//" "//t2c(this%grid%grid%xmax)//" "// &
3718 t2c(this%grid%grid%ymax))
3719#endif
3720
3721CASE ('UTM')
3722
3723 CALL grib_get(gaid,'zone',zone)
3724
3725 CALL grib_get(gaid,'datum',datum)
3726 IF (datum == 0) THEN
3727 CALL grib_get(gaid,'referenceLongitude',reflon)
3728 CALL grib_get(gaid,'falseEasting',this%grid%proj%xoff)
3729 CALL grib_get(gaid,'falseNorthing',this%grid%proj%yoff)
3731 ELSE
3733 CALL raise_fatal_error()
3734 ENDIF
3735
3736 CALL grib_get(gaid,'eastingOfFirstGridPoint',lofirst)
3737 CALL grib_get(gaid,'eastingOfLastGridPoint',lolast)
3738 CALL grib_get(gaid,'northingOfFirstGridPoint',lafirst)
3739 CALL grib_get(gaid,'northingOfLastGridPoint',lalast)
3740
3741 IF (iscansnegatively == 0) THEN
3742 this%grid%grid%xmin = lofirst
3743 this%grid%grid%xmax = lolast
3744 ELSE
3745 this%grid%grid%xmax = lofirst
3746 this%grid%grid%xmin = lolast
3747 ENDIF
3748 IF (jscanspositively == 0) THEN
3749 this%grid%grid%ymax = lafirst
3750 this%grid%grid%ymin = lalast
3751 ELSE
3752 this%grid%grid%ymin = lafirst
3753 this%grid%grid%ymax = lalast
3754 ENDIF
3755
3756! compute dx and dy (should we get them from grib?)
3757 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
3758
3759CASE default
3761 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
3762 CALL raise_error()
3763
3764END SELECT
3765
3766CONTAINS
3767! utilities routines for grib_api, is there a better place?
3768SUBROUTINE grib_get_dmiss(gaid, key, value)
3769INTEGER,INTENT(in) :: gaid
3770CHARACTER(len=*),INTENT(in) :: key
3771DOUBLE PRECISION,INTENT(out) :: value
3772
3773INTEGER :: ierr
3774
3775CALL grib_get(gaid, key, value, ierr)
3776IF (ierr /= grib_success) value = dmiss
3777
3778END SUBROUTINE grib_get_dmiss
3779
3780SUBROUTINE grib_get_imiss(gaid, key, value)
3781INTEGER,INTENT(in) :: gaid
3782CHARACTER(len=*),INTENT(in) :: key
3783INTEGER,INTENT(out) :: value
3784
3785INTEGER :: ierr
3786
3787CALL grib_get(gaid, key, value, ierr)
3788IF (ierr /= grib_success) value = imiss
3789
3790END SUBROUTINE grib_get_imiss
3791
3792
3793SUBROUTINE griddim_import_ellipsoid(this, gaid)
3794TYPE(griddim_def),INTENT(inout) :: this
3795INTEGER,INTENT(in) :: gaid
3796
3797INTEGER :: shapeofearth, iv, is
3798DOUBLE PRECISION :: r1, r2
3799
3800IF (editionnumber == 2) THEN
3801 CALL grib_get(gaid, 'shapeOfTheEarth', shapeofearth)
3802 SELECT CASE(shapeofearth)
3803 CASE(0) ! spherical
3805 CASE(1) ! spherical generic
3806 CALL grib_get(gaid, 'scaleFactorOfRadiusOfSphericalEarth', is)
3807 CALL grib_get(gaid, 'scaledValueOfRadiusOfSphericalEarth', iv)
3808 r1 = dble(iv) / 10**is
3810 CASE(2) ! iau65
3812 CASE(3,7) ! ellipsoidal generic
3813 CALL grib_get(gaid, 'scaleFactorOfEarthMajorAxis', is)
3814 CALL grib_get(gaid, 'scaledValueOfEarthMajorAxis', iv)
3815 r1 = dble(iv) / 10**is
3816 CALL grib_get(gaid, 'scaleFactorOfEarthMinorAxis', is)
3817 CALL grib_get(gaid, 'scaledValueOfEarthMinorAxis', iv)
3818 r2 = dble(iv) / 10**is
3819 IF (shapeofearth == 3) THEN ! km->m
3820 r1 = r1*1000.0d0
3821 r2 = r2*1000.0d0
3822 ENDIF
3823 IF (abs(r1) < 1.0d-6) THEN ! suspicious data read from grib
3825 'read from grib, going on with spherical Earth but the results may be wrong')
3827 ELSE
3829 ENDIF
3830 CASE(4) ! iag-grs80
3832 CASE(5) ! wgs84
3834 CASE(6) ! spherical
3836! CASE(7) ! google earth-like?
3837 CASE default
3839 t2c(shapeofearth)//' not supported in grib2')
3840 CALL raise_fatal_error()
3841
3842 END SELECT
3843
3844ELSE
3845
3846 CALL grib_get(gaid, 'earthIsOblate', shapeofearth)
3847 IF (shapeofearth == 0) THEN ! spherical
3849 ELSE ! iau65
3851 ENDIF
3852
3853ENDIF
3854
3855END SUBROUTINE griddim_import_ellipsoid
3856
3857
3858END SUBROUTINE griddim_import_gribapi
3859
3860
3861! grib_api driver
3862SUBROUTINE griddim_export_gribapi(this, gaid)
3863USE grib_api
3864TYPE(griddim_def),INTENT(in) :: this ! griddim object
3865INTEGER, INTENT(inout) :: gaid ! grib_api id of the grib loaded in memory to export
3866
3867INTEGER :: EditionNumber, iScansNegatively, jScansPositively, nv, pvl, zone, ierr
3868DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, reflon
3869DOUBLE PRECISION :: sdx, sdy, ratio, tol, ex, ey
3870
3871
3872! Generic keys
3873CALL grib_get(gaid,'GRIBEditionNumber', editionnumber)
3874! the following required since eccodes
3875IF (editionnumber == 2) CALL grib_set(gaid,'shapeOfTheEarth', 6)
3876CALL grib_set(gaid,'typeOfGrid', this%grid%proj%proj_type)
3877#ifdef DEBUG
3879 "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
3880#endif
3881
3882IF (this%grid%proj%proj_type == 'unstructured_grid') THEN
3883 CALL grib_set(gaid,'numberOfDataPoints', this%dim%nx)
3884! reenable when possible
3885! CALL griddim_export_ellipsoid(this, gaid)
3886 RETURN
3887ENDIF
3888
3889
3890! Edition dependent setup
3891IF (editionnumber == 1) THEN
3892 ratio = 1.d3
3893ELSE IF (editionnumber == 2) THEN
3894 ratio = 1.d6
3895ELSE
3896 ratio = 0.0d0 ! signal error?!
3897ENDIF
3898
3899! Keys valid for (almost?) all cases, Ni and Nj are universal aliases
3900CALL griddim_export_ellipsoid(this, gaid)
3901CALL grib_set(gaid, 'Ni', this%dim%nx)
3902CALL grib_set(gaid, 'Nj', this%dim%ny)
3903CALL grib_set(gaid,'uvRelativeToGrid',this%grid%grid%component_flag)
3904
3905CALL grib_get(gaid,'iScansNegatively',iscansnegatively)
3906CALL grib_get(gaid,'jScansPositively',jscanspositively)
3907
3908! Keys for rotated grids (checked through missing values and/or error code)
3909!SELECT CASE (this%grid%proj%proj_type)
3910!CASE ('rotated_ll', 'stretched_rotated_ll', 'polar_stereographic', 'lambert', 'albers')
3911CALL grib_set_dmiss(gaid,'longitudeOfSouthernPoleInDegrees', &
3912 this%grid%proj%rotated%longitude_south_pole, 0.0d0)
3913CALL grib_set_dmiss(gaid,'latitudeOfSouthernPoleInDegrees', &
3914 this%grid%proj%rotated%latitude_south_pole, -90.0d0)
3915IF (editionnumber == 1) THEN
3916 CALL grib_set_dmiss(gaid,'angleOfRotationInDegrees', &
3917 this%grid%proj%rotated%angle_rotation, 0.0d0)
3918ELSE IF (editionnumber == 2)THEN
3919 CALL grib_set_dmiss(gaid,'angleOfRotationOfProjectionInDegrees', &
3920 this%grid%proj%rotated%angle_rotation, 0.0d0)
3921ENDIF
3922
3923! Keys for stretched grids (checked through missing values and/or error code)
3924! units must be verified, still experimental in grib_api
3925! # TODO: Is it a float? Is it signed?
3926IF (editionnumber == 1) THEN
3927 CALL grib_set_dmiss(gaid,'longitudeOfStretchingPoleInDegrees', &
3928 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
3929 CALL grib_set_dmiss(gaid,'latitudeOfStretchingPoleInDegrees', &
3930 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
3931 CALL grib_set_dmiss(gaid,'stretchingFactor', &
3932 this%grid%proj%stretched%stretch_factor, 1.0d0)
3933ELSE IF (editionnumber == 2) THEN
3934 CALL grib_set_dmiss(gaid,'longitudeOfThePoleOfStretching', &
3935 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
3936 CALL grib_set_dmiss(gaid,'latitudeOfThePoleOfStretching', &
3937 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
3938 CALL grib_set_dmiss(gaid,'stretchingFactor', &
3939 this%grid%proj%stretched%stretch_factor, 1.0d6, 1.0d6)
3940ENDIF
3941
3942! Projection-dependent keys
3943SELECT CASE (this%grid%proj%proj_type)
3944
3945! Keys for sphaerical coordinate systems
3946CASE ('regular_ll', 'rotated_ll', 'stretched_ll', 'stretched_rotated_ll')
3947
3948 IF (iscansnegatively == 0) THEN
3949 lofirst = this%grid%grid%xmin
3950 lolast = this%grid%grid%xmax
3951 ELSE
3952 lofirst = this%grid%grid%xmax
3953 lolast = this%grid%grid%xmin
3954 ENDIF
3955 IF (jscanspositively == 0) THEN
3956 lafirst = this%grid%grid%ymax
3957 lalast = this%grid%grid%ymin
3958 ELSE
3959 lafirst = this%grid%grid%ymin
3960 lalast = this%grid%grid%ymax
3961 ENDIF
3962
3963! reset lon in standard grib 2 definition [0,360]
3964 IF (editionnumber == 1) THEN
3965 CALL long_reset_m180_360(lofirst)
3966 CALL long_reset_m180_360(lolast)
3967 ELSE IF (editionnumber == 2) THEN
3968 CALL long_reset_0_360(lofirst)
3969 CALL long_reset_0_360(lolast)
3970 ENDIF
3971
3972 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
3973 CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
3974 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
3975 CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
3976
3977! test relative coordinate truncation error with respect to tol
3978! tol should be tuned
3979 sdx = this%grid%grid%dx*ratio
3980 sdy = this%grid%grid%dy*ratio
3981! error is computed relatively to the whole grid extension
3982 ex = abs(((sdx - nint(sdx))/ratio*this%dim%nx)/ &
3983 (this%grid%grid%xmax-this%grid%grid%xmin))
3984 ey = abs(((sdy - nint(sdy))/ratio*this%dim%ny)/ &
3985 (this%grid%grid%ymax-this%grid%grid%ymin))
3986 tol = 1.0d-3
3987 IF (ex > tol .OR. ey > tol) THEN
3988#ifdef DEBUG
3990 "griddim_export_gribapi, error on x "//&
3991 trim(to_char(ex))//"/"//t2c(tol))
3993 "griddim_export_gribapi, error on y "//&
3994 trim(to_char(ey))//"/"//t2c(tol))
3995#endif
3996! previous frmula relative to a single grid step
3997! tol = 1.0d0/ratio
3998! IF (ABS(NINT(sdx)/sdx - 1.0d0) > tol .OR. ABS(NINT(sdy)/sdy - 1.0d0) > tol) THEN
3999!#ifdef DEBUG
4000! CALL l4f_category_log(this%category,L4F_DEBUG, &
4001! "griddim_export_gribapi, dlon relative error: "//&
4002! TRIM(to_char(ABS(NINT(sdx)/sdx - 1.0d0)))//">"//TRIM(to_char(tol)))
4003! CALL l4f_category_log(this%category,L4F_DEBUG, &
4004! "griddim_export_gribapi, dlat relative error: "//&
4005! TRIM(to_char(ABS(NINT(sdy)/sdy - 1.0d0)))//">"//TRIM(to_char(tol)))
4006!#endif
4008 "griddim_export_gribapi, increments not given: inaccurate!")
4009 CALL grib_set_missing(gaid,'Di')
4010 CALL grib_set_missing(gaid,'Dj')
4011 CALL grib_set(gaid,'ijDirectionIncrementGiven',0)
4012 ELSE
4013#ifdef DEBUG
4015 "griddim_export_gribapi, setting increments: "// &
4016 trim(to_char(this%grid%grid%dx))//' '//trim(to_char(this%grid%grid%dy)))
4017#endif
4018 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
4019 CALL grib_set(gaid,'iDirectionIncrement',nint(this%grid%grid%dx*ratio))
4020 CALL grib_set(gaid,'jDirectionIncrement',nint(this%grid%grid%dy*ratio))
4021! this does not work in grib_set
4022! CALL grib_set(gaid,'iDirectionIncrementInDegrees',this%grid%grid%dx)
4023! CALL grib_set(gaid,'jDirectionIncrementInDegrees',this%grid%grid%dy)
4024 ENDIF
4025
4026! Keys for polar projections
4027CASE ('polar_stereographic', 'lambert', 'albers')
4028! increments are required
4029 CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
4030 CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
4031 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
4032! latin1/latin2 may be missing (e.g. stereographic)
4033 CALL grib_set_dmiss(gaid,'Latin1InDegrees',this%grid%proj%polar%latin1)
4034 CALL grib_set_dmiss(gaid,'Latin2InDegrees',this%grid%proj%polar%latin2)
4035! projection center flag, aka hemisphere
4036 CALL grib_set(gaid,'projectionCenterFlag',&
4037 this%grid%proj%polar%projection_center_flag, ierr)
4038 IF (ierr /= grib_success) THEN ! try center/centre
4039 CALL grib_set(gaid,'projectionCentreFlag',&
4040 this%grid%proj%polar%projection_center_flag)
4041 ENDIF
4042
4043
4044! line of view, aka central meridian
4045 CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
4046! latitude at which dx and dy are valid
4047 IF (editionnumber == 2) THEN
4048 CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
4049 ENDIF
4050
4051! compute lon and lat of first point from projected extremes
4052 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
4053 lofirst, lafirst)
4054! reset lon in standard grib 2 definition [0,360]
4055 IF (editionnumber == 1) THEN
4056 CALL long_reset_m180_360(lofirst)
4057 ELSE IF (editionnumber == 2) THEN
4058 CALL long_reset_0_360(lofirst)
4059 ENDIF
4060 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
4061 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
4062
4063! Keys for equatorial projections
4064CASE ('mercator')
4065
4066! increments are required
4067 CALL grib_set(gaid,'DxInMetres', this%grid%grid%dx)
4068 CALL grib_set(gaid,'DyInMetres', this%grid%grid%dy)
4069 CALL grib_set(gaid,'ijDirectionIncrementGiven',1)
4070
4071! line of view, aka central meridian, not in grib
4072! CALL grib_set(gaid,'LoVInDegrees',this%grid%proj%lov)
4073! latitude at which dx and dy are valid
4074 CALL grib_set(gaid,'LaDInDegrees',this%grid%proj%polar%lad)
4075
4076! compute lon and lat of first and last points from projected extremes
4077 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
4078 lofirst, lafirst)
4079 CALL grib_set(gaid,'longitudeOfFirstGridPointInDegrees',lofirst)
4080 CALL grib_set(gaid,'latitudeOfFirstGridPointInDegrees',lafirst)
4081 CALL get_unproj_last(this, iscansnegatively, jscanspositively, &
4082 lolast, lalast)
4083 CALL grib_set(gaid,'longitudeOfLastGridPointInDegrees',lolast)
4084 CALL grib_set(gaid,'latitudeOfLastGridPointInDegrees',lalast)
4085
4086 IF (editionnumber == 2) THEN
4087 CALL grib_set(gaid,'orientationOfTheGridInDegrees',0.)
4088 ENDIF
4089
4090CASE ('UTM')
4091
4092 CALL grib_set(gaid,'datum',0)
4094 CALL grib_set(gaid,'referenceLongitude',nint(reflon/1.0d6))
4095 CALL grib_set(gaid,'falseEasting',this%grid%proj%xoff)
4096 CALL grib_set(gaid,'falseNorthing',this%grid%proj%yoff)
4097
4098 CALL grib_set(gaid,'iDirectionIncrement',this%grid%grid%dx)
4099 CALL grib_set(gaid,'jDirectionIncrement',this%grid%grid%dy)
4100
4101!res/scann ??
4102
4103 CALL grib_set(gaid,'zone',zone)
4104
4105 IF (iscansnegatively == 0) THEN
4106 lofirst = this%grid%grid%xmin
4107 lolast = this%grid%grid%xmax
4108 ELSE
4109 lofirst = this%grid%grid%xmax
4110 lolast = this%grid%grid%xmin
4111 ENDIF
4112 IF (jscanspositively == 0) THEN
4113 lafirst = this%grid%grid%ymax
4114 lalast = this%grid%grid%ymin
4115 ELSE
4116 lafirst = this%grid%grid%ymin
4117 lalast = this%grid%grid%ymax
4118 ENDIF
4119
4120 CALL grib_set(gaid,'eastingOfFirstGridPoint',lofirst)
4121 CALL grib_set(gaid,'eastingOfLastGridPoint',lolast)
4122 CALL grib_set(gaid,'northingOfFirstGridPoint',lafirst)
4123 CALL grib_set(gaid,'northingOfLastGridPoint',lalast)
4124
4125CASE default
4127 "griddim_export_gribapi, grid type "//trim(this%grid%proj%proj_type)//" not supported")
4128 CALL raise_error()
4129
4130END SELECT
4131
4132! hack for position of vertical coordinate parameters
4133! buggy in grib_api
4134IF (editionnumber == 1) THEN
4135! CALL grib_get(gaid,"PVPresent",pvp) ! alias, probably useless
4136 CALL grib_get(gaid,"NV",nv)
4137#ifdef DEBUG
4139 trim(to_char(nv))//" vertical coordinate parameters")
4140#endif
4141
4142 IF (nv == 0) THEN
4143 pvl = 255
4144 ELSE
4145 SELECT CASE (this%grid%proj%proj_type)
4146 CASE ('regular_ll') ! check whether "29-32 Set to zero (reserved)" are required
4147 pvl = 33
4148 CASE ('polar_stereographic')
4149 pvl = 33
4150 CASE ('rotated_ll', 'stretched_ll', 'lambert', 'albers')
4151 pvl = 43
4152 CASE ('stretched_rotated_ll')
4153 pvl = 43
4154 CASE DEFAULT
4155 pvl = 43 !?
4156 END SELECT
4157 ENDIF
4158
4159 CALL grib_set(gaid,"pvlLocation",pvl)
4160#ifdef DEBUG
4162 trim(to_char(pvl))//" as vertical coordinate parameter location")
4163#endif
4164
4165ENDIF
4166
4167
4168CONTAINS
4169! utilities routines for grib_api, is there a better place?
4170SUBROUTINE grib_set_dmiss(gaid, key, val, default, factor)
4171INTEGER,INTENT(in) :: gaid
4172CHARACTER(len=*),INTENT(in) :: key
4173DOUBLE PRECISION,INTENT(in) :: val
4174DOUBLE PRECISION,INTENT(in),OPTIONAL :: default
4175DOUBLE PRECISION,INTENT(in),OPTIONAL :: factor
4176
4177INTEGER :: ierr
4178
4180 IF (PRESENT(factor)) THEN
4181 CALL grib_set(gaid, key, val*factor, ierr)
4182 ELSE
4183 CALL grib_set(gaid, key, val, ierr)
4184 ENDIF
4185ELSE IF (PRESENT(default)) THEN
4186 CALL grib_set(gaid, key, default, ierr)
4187ENDIF
4188
4189END SUBROUTINE grib_set_dmiss
4190
4191SUBROUTINE grib_set_imiss(gaid, key, value, default)
4192INTEGER,INTENT(in) :: gaid
4193CHARACTER(len=*),INTENT(in) :: key
4194INTEGER,INTENT(in) :: value
4195INTEGER,INTENT(in),OPTIONAL :: default
4196
4197INTEGER :: ierr
4198
4200 CALL grib_set(gaid, key, value, ierr)
4201ELSE IF (PRESENT(default)) THEN
4202 CALL grib_set(gaid, key, default, ierr)
4203ENDIF
4204
4205END SUBROUTINE grib_set_imiss
4206
4207SUBROUTINE griddim_export_ellipsoid(this, gaid)
4208TYPE(griddim_def),INTENT(in) :: this
4209INTEGER,INTENT(in) :: gaid
4210
4211INTEGER :: ellips_type, ierr
4212DOUBLE PRECISION :: r1, r2, f
4213
4215
4216IF (editionnumber == 2) THEN
4217
4218! why does it produce a message even with ierr?
4219 CALL grib_set_missing(gaid, 'scaleFactorOfRadiusOfSphericalEarth', ierr)
4220 CALL grib_set_missing(gaid, 'scaledValueOfRadiusOfSphericalEarth', ierr)
4221 CALL grib_set_missing(gaid, 'scaleFactorOfEarthMajorAxis', ierr)
4222 CALL grib_set_missing(gaid, 'scaledValueOfEarthMajorAxis', ierr)
4223 CALL grib_set_missing(gaid, 'scaleFactorOfEarthMinorAxis', ierr)
4224 CALL grib_set_missing(gaid, 'scaledValueOfEarthMinorAxis', ierr)
4225
4226 SELECT CASE(ellips_type)
4227 CASE(ellips_grs80) ! iag-grs80
4228 CALL grib_set(gaid, 'shapeOfTheEarth', 4)
4229 CASE(ellips_wgs84) ! wgs84
4230 CALL grib_set(gaid, 'shapeOfTheEarth', 5)
4231 CASE default
4232 IF (f == 0.0d0) THEN ! spherical Earth
4233 IF (r1 == 6367470.0d0) THEN ! spherical
4234 CALL grib_set(gaid, 'shapeOfTheEarth', 0)
4235 ELSE IF (r1 == 6371229.0d0) THEN ! spherical
4236 CALL grib_set(gaid, 'shapeOfTheEarth', 6)
4237 ELSE ! spherical generic
4238 CALL grib_set(gaid, 'shapeOfTheEarth', 1)
4239 CALL grib_set(gaid, 'scaleFactorOfRadiusOfSphericalEarth', 2)
4240 CALL grib_set(gaid, 'scaledValueOfRadiusOfSphericalEarth', int(r1*100.0d0))
4241 ENDIF
4242 ELSE ! ellipsoidal
4243 IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
4244 CALL grib_set(gaid, 'shapeOfTheEarth', 2)
4245 ELSE ! ellipsoidal generic
4246 CALL grib_set(gaid, 'shapeOfTheEarth', 3)
4247 r2 = r1*(1.0d0 - f)
4248 CALL grib_set(gaid, 'scaleFactorOfEarthMajorAxis', 5)
4249 CALL grib_set(gaid, 'scaledValueOfEarthMajorAxis', &
4250 int(r1*100.0d0))
4251 CALL grib_set(gaid, 'scaleFactorOfEarthMinorAxis', 5)
4252 CALL grib_set(gaid, 'scaledValueOfEarthMinorAxis', &
4253 int(r2*100.0d0))
4254 ENDIF
4255 ENDIF
4256 END SELECT
4257
4258ELSE
4259
4260 IF (r1 == 6367470.0d0 .AND. f == 0.0d0) THEN ! spherical
4261 CALL grib_set(gaid, 'earthIsOblate', 0)
4262 ELSE IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0) THEN ! iau65
4263 CALL grib_set(gaid, 'earthIsOblate', 1)
4264 ELSE IF (f == 0.0d0) THEN ! generic spherical
4265 CALL grib_set(gaid, 'earthIsOblate', 0)
4267 ¬ supported in grib 1, coding with default radius of 6367470 m')
4268 ELSE ! generic ellipsoidal
4269 CALL grib_set(gaid, 'earthIsOblate', 1)
4271 ¬ supported in grib 1, coding with default iau65 ellipsoid')
4272 ENDIF
4273
4274ENDIF
4275
4276END SUBROUTINE griddim_export_ellipsoid
4277
4278SUBROUTINE get_unproj_first(this, iScansNegatively, jScansPositively, &
4279 loFirst, laFirst)
4280TYPE(griddim_def),INTENT(in) :: this ! griddim object
4281INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
4282DOUBLE PRECISION,INTENT(out) :: loFirst, laFirst
4283
4284! compute lon and lat of first point from projected extremes
4285IF (iscansnegatively == 0) THEN
4286 IF (jscanspositively == 0) THEN
4288 ELSE
4290 ENDIF
4291ELSE
4292 IF (jscanspositively == 0) THEN
4294 ELSE
4296 ENDIF
4297ENDIF
4298! use the values kept for personal pleasure ?
4299! loFirst = this%grid%proj%polar%lon1
4300! laFirst = this%grid%proj%polar%lat1
4301END SUBROUTINE get_unproj_first
4302
4303SUBROUTINE get_unproj_last(this, iScansNegatively, jScansPositively, &
4304 loLast, laLast)
4305TYPE(griddim_def),INTENT(in) :: this ! griddim object
4306INTEGER,INTENT(in) :: iScansNegatively, jScansPositively
4307DOUBLE PRECISION,INTENT(out) :: loLast, laLast
4308
4309! compute lon and lat of last point from projected extremes
4310IF (iscansnegatively == 0) THEN
4311 IF (jscanspositively == 0) THEN
4313 ELSE
4315 ENDIF
4316ELSE
4317 IF (jscanspositively == 0) THEN
4319 ELSE
4321 ENDIF
4322ENDIF
4323! use the values kept for personal pleasure ?
4324! loLast = this%grid%proj%polar%lon?
4325! laLast = this%grid%proj%polar%lat?
4326END SUBROUTINE get_unproj_last
4327
4328END SUBROUTINE griddim_export_gribapi
4329#endif
4330
4331
4332#ifdef HAVE_LIBGDAL
4333! gdal driver
4334SUBROUTINE griddim_import_gdal(this, gdalid, gdal_options)
4335USE gdal
4336TYPE(griddim_def),INTENT(inout) :: this ! griddim object
4337TYPE(gdalrasterbandh),INTENT(in) :: gdalid ! gdal rasterband pointer
4338TYPE(gdal_file_id_options),INTENT(in) :: gdal_options
4339
4340TYPE(gdaldataseth) :: hds
4341REAL(kind=c_double) :: geotrans(6) !, invgeotrans(6), x1, y1, x2, y2
4342INTEGER(kind=c_int) :: offsetx, offsety
4343INTEGER :: ier
4344
4345hds = gdalgetbanddataset(gdalid) ! go back to dataset
4346ier = gdalgetgeotransform(hds, geotrans)
4347
4348IF (ier /= 0) THEN
4350 'griddim_import_gdal: error in accessing gdal dataset')
4351 CALL raise_error()
4352 RETURN
4353ENDIF
4354IF (geotrans(3) /= 0.0_c_double .OR. geotrans(5) /= 0.0_c_double) THEN ! dataset has not a diagonal transformation
4356 'griddim_import_gdal: dataset has a non-diagonal transformation matrix, unsupported')
4357 CALL raise_error()
4358 RETURN
4359ENDIF
4360
4361CALL gdaldatasetbbsize_f(hds, gdal_options%xmin, gdal_options%ymin, &
4362 gdal_options%xmax, gdal_options%ymax, &
4363 this%dim%nx, this%dim%ny, offsetx, offsety, &
4364 this%grid%grid%xmin, this%grid%grid%ymin, this%grid%grid%xmax, this%grid%grid%ymax)
4365
4366IF (this%dim%nx == 0 .OR. this%dim%ny == 0) THEN
4368 'griddim_import_gdal: requested bounding box '//t2c(gdal_options%xmin)//','// &
4369 t2c(gdal_options%ymin)//','//t2c(gdal_options%xmax)//','//&
4370 t2c(gdal_options%ymax))
4372 'determines an empty dataset '//t2c(this%dim%nx)//'x'//t2c(this%dim%ny))
4373ENDIF
4374
4375! get grid corners
4376!CALL gdalapplygeotransform(geotrans, 0.5_c_double, 0.5_c_double, x1, y1)
4377!CALL gdalapplygeotransform(geotrans, gdalgetrasterbandxsize(gdalid)-0.5_c_double, &
4378! gdalgetrasterbandysize(gdalid)-0.5_c_double, x2, y2)
4379
4380!IF (geotrans(3) == 0.0_c_double .AND. geotrans(5) == 0.0_c_double) THEN ! transformation is diagonal, no transposing
4381! this%dim%nx = gdalgetrasterbandxsize(gdalid)
4382! this%dim%ny = gdalgetrasterbandysize(gdalid)
4383! this%grid%grid%xmin = MIN(x1, x2)
4384! this%grid%grid%xmax = MAX(x1, x2)
4385! this%grid%grid%ymin = MIN(y1, y2)
4386! this%grid%grid%ymax = MAX(y1, y2)
4387!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
4388!
4389! this%dim%nx = gdalgetrasterbandysize(gdalid)
4390! this%dim%ny = gdalgetrasterbandxsize(gdalid)
4391! this%grid%grid%xmin = MIN(y1, y2)
4392! this%grid%grid%xmax = MAX(y1, y2)
4393! this%grid%grid%ymin = MIN(x1, x2)
4394! this%grid%grid%ymax = MAX(x1, x2)
4395!
4396!ELSE ! transformation is a rotation, not supported
4397!ENDIF
4398
4399this%grid%proj%proj_type = 'regular_ll' ! forced, only one supported (check?)
4400
4401CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
4402this%grid%grid%component_flag = 0
4403
4404END SUBROUTINE griddim_import_gdal
4405#endif
4406
4407
4409SUBROUTINE griddim_display(this)
4410TYPE(griddim_def),INTENT(in) :: this
4411
4412print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
4413
4417
4418print*,"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
4419
4420END SUBROUTINE griddim_display
4421
4422
4423#define VOL7D_POLY_TYPE TYPE(grid_def)
4424#define VOL7D_POLY_TYPES _grid
4425#include "array_utilities_inc.F90"
4426#undef VOL7D_POLY_TYPE
4427#undef VOL7D_POLY_TYPES
4428
4429#define VOL7D_POLY_TYPE TYPE(griddim_def)
4430#define VOL7D_POLY_TYPES _griddim
4431#include "array_utilities_inc.F90"
4432#undef VOL7D_POLY_TYPE
4433#undef VOL7D_POLY_TYPES
4434
4435
4447SUBROUTINE griddim_wind_unrot(this, rot_mat)
4449TYPE(griddim_def), INTENT(in) :: this
4450DOUBLE PRECISION, POINTER :: rot_mat(:,:,:)
4451
4452DOUBLE PRECISION :: dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j
4453DOUBLE PRECISION :: lat_factor
4454INTEGER :: i, j, ip1, im1, jp1, jm1;
4455
4456IF (this%dim%nx <= 0 .OR. this%dim%ny <= 0 .OR. &
4457 .NOT. ASSOCIATED(this%dim%lon) .OR. .NOT. ASSOCIATED(this%dim%lat)) THEN
4458 CALL l4f_category_log(this%category, l4f_error, 'rotate_uv must be called after correct init of griddim object')
4459 NULLIFY(rot_mat)
4460 RETURN
4461ENDIF
4462
4463ALLOCATE(rot_mat(this%dim%nx, this%dim%ny, 4))
4464
4465DO j = 1, this%dim%ny
4466 jp1 = min(j+1, this%dim%ny)
4467 jm1 = max(j-1, 1)
4468 DO i = 1, this%dim%nx
4469 ip1 = min(i+1, this%dim%nx)
4470 im1 = max(i-1, 1)
4471
4472 dlat_i = this%dim%lat(ip1,j) - this%dim%lat(im1,j)
4473 dlat_j = this%dim%lat(i,jp1) - this%dim%lat(i,jm1)
4474
4475 dlon_i = this%dim%lon(ip1,j) - this%dim%lon(im1,j)
4476! if ( dlon_i > pi ) dlon_i -= 2*pi;
4477! if ( dlon_i < (-pi) ) dlon_i += 2*pi;
4478 dlon_j = this%dim%lon(i,jp1) - this%dim%lon(i,jm1)
4479! if ( dlon_j > pi ) dlon_j -= 2*pi;
4480! if ( dlon_j < (-pi) ) dlon_j += 2*pi;
4481
4482! check whether this is really necessary !!!!
4483 lat_factor = cos(degrad*this%dim%lat(i,j))
4484 dlon_i = dlon_i * lat_factor
4485 dlon_j = dlon_j * lat_factor
4486
4487 dist_i = sqrt(dlon_i*dlon_i + dlat_i*dlat_i)
4488 dist_j = sqrt(dlon_j*dlon_j + dlat_j*dlat_j)
4489
4490 IF (dist_i > 0.d0) THEN
4491 rot_mat(i,j,1) = dlon_i/dist_i
4492 rot_mat(i,j,2) = dlat_i/dist_i
4493 ELSE
4494 rot_mat(i,j,1) = 0.0d0
4495 rot_mat(i,j,2) = 0.0d0
4496 ENDIF
4497 IF (dist_j > 0.d0) THEN
4498 rot_mat(i,j,3) = dlon_j/dist_j
4499 rot_mat(i,j,4) = dlat_j/dist_j
4500 ELSE
4501 rot_mat(i,j,3) = 0.0d0
4502 rot_mat(i,j,4) = 0.0d0
4503 ENDIF
4504
4505 ENDDO
4506ENDDO
4507
4508END SUBROUTINE griddim_wind_unrot
4509
4510
4511! compute zoom indices from geographical zoom coordinates
4512SUBROUTINE griddim_zoom_coord(this, ilon, ilat, flon, flat, ix, iy, fx, fy)
4513TYPE(griddim_def),INTENT(in) :: this
4514DOUBLE PRECISION,INTENT(in) :: ilon, ilat, flon, flat
4515INTEGER,INTENT(out) :: ix, iy, fx, fy
4516
4517DOUBLE PRECISION :: ix1, iy1, fx1, fy1
4518
4519! compute projected coordinates of vertices of desired lonlat rectangle
4522
4523CALL griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
4524
4525END SUBROUTINE griddim_zoom_coord
4526
4527
4528! compute zoom indices from projected zoom coordinates
4529SUBROUTINE griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
4530TYPE(griddim_def),INTENT(in) :: this
4531DOUBLE PRECISION,INTENT(in) :: ix1, iy1, fx1, fy1
4532INTEGER,INTENT(out) :: ix, iy, fx, fy
4533
4534INTEGER :: lix, liy, lfx, lfy
4535
4536! compute projected indices
4537lix = nint((ix1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
4538liy = nint((iy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
4539lfx = nint((fx1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
4540lfy = nint((fy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
4541! swap projected indices, in case grid is "strongly rotated"
4542ix = min(lix, lfx)
4543fx = max(lix, lfx)
4544iy = min(liy, lfy)
4545fy = max(liy, lfy)
4546
4547END SUBROUTINE griddim_zoom_projcoord
4548
4549
4553SUBROUTINE long_reset_0_360(lon)
4554DOUBLE PRECISION,INTENT(inout) :: lon
4555
4557DO WHILE(lon < 0.0d0)
4558 lon = lon + 360.0d0
4559END DO
4560DO WHILE(lon >= 360.0d0)
4561 lon = lon - 360.0d0
4562END DO
4563
4564END SUBROUTINE long_reset_0_360
4565
4566
4570SUBROUTINE long_reset_m180_360(lon)
4571DOUBLE PRECISION,INTENT(inout) :: lon
4572
4574DO WHILE(lon < -180.0d0)
4575 lon = lon + 360.0d0
4576END DO
4577DO WHILE(lon >= 360.0d0)
4578 lon = lon - 360.0d0
4579END DO
4580
4581END SUBROUTINE long_reset_m180_360
4582
4583
4587!SUBROUTINE long_reset_m90_270(lon)
4588!DOUBLE PRECISION,INTENT(inout) :: lon !< the longitude to reset
4589!
4590!IF (.NOT.c_e(lon)) RETURN
4591!DO WHILE(lon < -90.0D0)
4592! lon = lon + 360.0D0
4593!END DO
4594!DO WHILE(lon >= 270.0D0)
4595! lon = lon - 360.0D0
4596!END DO
4597!
4598!END SUBROUTINE long_reset_m90_270
4599
4600
4604SUBROUTINE long_reset_m180_180(lon)
4605DOUBLE PRECISION,INTENT(inout) :: lon
4606
4608DO WHILE(lon < -180.0d0)
4609 lon = lon + 360.0d0
4610END DO
4611DO WHILE(lon >= 180.0d0)
4612 lon = lon - 360.0d0
4613END DO
4614
4615END SUBROUTINE long_reset_m180_180
4616
4617
4618SUBROUTINE long_reset_to_cart_closest(lon, lonref)
4619DOUBLE PRECISION,INTENT(inout) :: lon
4620DOUBLE PRECISION,INTENT(in) :: lonref
4621
4623IF (abs(lon-lonref) < 180.0d0) RETURN ! nothing to do
4624lon = lon - nint((lon-lonref)/360.0d0)*360.0d0 ! se non e` vera e` ben trovata
4625
4626END SUBROUTINE long_reset_to_cart_closest
4627
4628
4630
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 |