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