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