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