61 character (len=255),
parameter:: subcategory=
"grid_class" 71 type(geo_proj) :: proj
72 type(grid_rect) :: grid
73 integer :: category = 0
83 type(grid_def) :: grid
85 integer :: category = 0
92 INTERFACE OPERATOR (==)
93 MODULE PROCEDURE grid_eq, griddim_eq
99 INTERFACE OPERATOR (/=)
100 MODULE PROCEDURE grid_ne, griddim_ne
105 MODULE PROCEDURE griddim_init
110 MODULE PROCEDURE griddim_delete
115 MODULE PROCEDURE griddim_copy
121 MODULE PROCEDURE griddim_coord_proj
127 MODULE PROCEDURE griddim_coord_unproj, griddim_unproj
132 MODULE PROCEDURE griddim_get_val
137 MODULE PROCEDURE griddim_set_val
142 MODULE PROCEDURE griddim_write_unit
147 MODULE PROCEDURE griddim_read_unit
152 MODULE PROCEDURE griddim_import_grid_id
157 MODULE PROCEDURE griddim_export_grid_id
162 MODULE PROCEDURE griddim_display
165 #define VOL7D_POLY_TYPE TYPE(grid_def) 166 #define VOL7D_POLY_TYPES _grid 167 #include "array_utilities_pre.F90" 168 #undef VOL7D_POLY_TYPE 169 #undef VOL7D_POLY_TYPES 171 #define VOL7D_POLY_TYPE TYPE(griddim_def) 172 #define VOL7D_POLY_TYPES _griddim 173 #include "array_utilities_pre.F90" 174 #undef VOL7D_POLY_TYPE 175 #undef VOL7D_POLY_TYPES 178 MODULE PROCEDURE griddim_wind_unrot
184 PUBLIC proj,
unproj, griddim_unproj, griddim_gen_coord, &
185 griddim_zoom_coord, griddim_zoom_projcoord, &
189 PUBLIC OPERATOR(==),
OPERATOR(/=)
190 PUBLIC count_distinct, pack_distinct, count_and_pack_distinct, &
191 map_distinct, map_inv_distinct,
index 193 PUBLIC griddim_central_lon, griddim_set_central_lon
197 SUBROUTINE griddim_init(this, nx, ny, &
198 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
199 proj_type, lov, zone, xoff, yoff, &
200 longitude_south_pole, latitude_south_pole, angle_rotation, &
201 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
202 latin1, latin2, lad, projection_center_flag, &
203 ellips_smaj_axis, ellips_flatt, ellips_type, &
205 TYPE(griddim_def),
INTENT(inout) :: this
206 INTEGER,
INTENT(in),
OPTIONAL :: nx
207 INTEGER,
INTENT(in),
OPTIONAL :: ny
208 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xmin
209 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xmax
210 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ymin
211 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ymax
212 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: dx
213 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: dy
216 INTEGER,
INTENT(in),
OPTIONAL :: component_flag
217 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: proj_type
218 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lov
219 INTEGER,
INTENT(in),
OPTIONAL :: zone
220 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xoff
221 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: yoff
222 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_south_pole
223 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_south_pole
224 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: angle_rotation
225 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_stretch_pole
226 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_stretch_pole
227 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: stretch_factor
228 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin1
229 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin2
230 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lad
231 INTEGER,
INTENT(in),
OPTIONAL :: projection_center_flag
232 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_smaj_axis
233 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_flatt
234 INTEGER,
INTENT(in),
OPTIONAL :: ellips_type
235 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
237 CHARACTER(len=512) :: a_name
239 IF (
PRESENT(categoryappend))
THEN 240 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
241 trim(categoryappend))
243 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
245 this%category=l4f_category_get(a_name)
248 this%grid%proj = geo_proj_new( &
249 proj_type=proj_type, lov=lov, zone=zone, xoff=xoff, yoff=yoff, &
250 longitude_south_pole=longitude_south_pole, &
251 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
252 longitude_stretch_pole=longitude_stretch_pole, &
253 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
254 latin1=latin1, latin2=latin2, lad=lad, projection_center_flag=projection_center_flag, &
255 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, ellips_type=ellips_type)
257 this%grid%grid = grid_rect_new( &
258 xmin, xmax, ymin, ymax, dx, dy, component_flag)
260 this%dim = grid_dim_new(nx, ny)
263 call l4f_category_log(this%category,l4f_debug,
"init gtype: "//this%grid%proj%proj_type )
266 END SUBROUTINE griddim_init
270 SUBROUTINE griddim_delete(this)
271 TYPE(griddim_def),
INTENT(inout) :: this
274 CALL delete(this%grid%proj)
275 CALL delete(this%grid%grid)
277 call l4f_category_delete(this%category)
279 END SUBROUTINE griddim_delete
283 SUBROUTINE griddim_copy(this, that, categoryappend)
284 TYPE(griddim_def),
INTENT(in) :: this
285 TYPE(griddim_def),
INTENT(out) :: that
286 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
288 CHARACTER(len=512) :: a_name
292 CALL copy(this%grid%proj, that%grid%proj)
293 CALL copy(this%grid%grid, that%grid%grid)
294 CALL copy(this%dim, that%dim)
297 IF (
PRESENT(categoryappend))
THEN 298 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."// &
299 trim(categoryappend))
301 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
303 that%category=l4f_category_get(a_name)
305 END SUBROUTINE griddim_copy
310 ELEMENTAL SUBROUTINE griddim_coord_proj(this, lon, lat, x, y)
311 TYPE(griddim_def),
INTENT(in) :: this
313 DOUBLE PRECISION,
INTENT(in) :: lon, lat
315 DOUBLE PRECISION,
INTENT(out) :: x, y
317 CALL proj(this%grid%proj, lon, lat, x, y)
319 END SUBROUTINE griddim_coord_proj
324 ELEMENTAL SUBROUTINE griddim_coord_unproj(this, x, y, lon, lat)
325 TYPE(griddim_def),
INTENT(in) :: this
327 DOUBLE PRECISION,
INTENT(in) :: x, y
329 DOUBLE PRECISION,
INTENT(out) :: lon, lat
331 CALL unproj(this%grid%proj, x, y, lon, lat)
333 END SUBROUTINE griddim_coord_unproj
358 SUBROUTINE griddim_unproj(this)
359 TYPE(griddim_def),
INTENT(inout) :: this
361 IF (.NOT.
c_e(this%dim%nx) .OR. .NOT.
c_e(this%dim%ny))
RETURN 363 CALL griddim_unproj_internal(this)
365 END SUBROUTINE griddim_unproj
368 SUBROUTINE griddim_unproj_internal(this)
369 TYPE(griddim_def),
INTENT(inout) ::this
371 DOUBLE PRECISION :: x(this%dim%nx,this%dim%ny), y(this%dim%nx,this%dim%ny)
373 CALL grid_rect_coordinates(this%grid%grid, x, y)
374 CALL unproj(this, x, y, this%dim%lon, this%dim%lat)
376 END SUBROUTINE griddim_unproj_internal
380 SUBROUTINE griddim_get_val(this, nx, ny, &
381 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
382 proj, proj_type, lov, zone, xoff, yoff, &
383 longitude_south_pole, latitude_south_pole, angle_rotation, &
384 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
385 latin1, latin2, lad, projection_center_flag, &
386 ellips_smaj_axis, ellips_flatt, ellips_type)
387 TYPE(griddim_def),
INTENT(in) :: this
388 INTEGER,
INTENT(out),
OPTIONAL :: nx
389 INTEGER,
INTENT(out),
OPTIONAL :: ny
391 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: xmin, xmax, ymin, ymax
393 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: dx, dy
396 INTEGER,
INTENT(out),
OPTIONAL :: component_flag
397 TYPE(geo_proj),
INTENT(out),
OPTIONAL :: proj
398 CHARACTER(len=*),
INTENT(out),
OPTIONAL :: proj_type
399 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: lov
400 INTEGER,
INTENT(out),
OPTIONAL :: zone
401 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: xoff
402 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: yoff
403 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: longitude_south_pole
404 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: latitude_south_pole
405 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: angle_rotation
406 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: longitude_stretch_pole
407 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: latitude_stretch_pole
408 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: stretch_factor
409 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: latin1
410 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: latin2
411 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: lad
412 INTEGER,
INTENT(out),
OPTIONAL :: projection_center_flag
413 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: ellips_smaj_axis
414 DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: ellips_flatt
415 INTEGER,
INTENT(out),
OPTIONAL :: ellips_type
417 IF (
PRESENT(nx)) nx = this%dim%nx
418 IF (
PRESENT(ny)) ny = this%dim%ny
420 IF (
PRESENT(
proj))
proj = this%grid%proj
422 CALL get_val(this%grid%proj, proj_type=proj_type, lov=lov, zone=zone, &
423 xoff=xoff, yoff=yoff, &
424 longitude_south_pole=longitude_south_pole, &
425 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
426 longitude_stretch_pole=longitude_stretch_pole, &
427 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
428 latin1=latin1, latin2=latin2, lad=lad, &
429 projection_center_flag=projection_center_flag, &
430 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
431 ellips_type=ellips_type)
434 xmin, xmax, ymin, ymax, dx, dy, component_flag)
436 END SUBROUTINE griddim_get_val
440 SUBROUTINE griddim_set_val(this, nx, ny, &
441 xmin, xmax, ymin, ymax, dx, dy, component_flag, &
442 proj_type, lov, zone, xoff, yoff, &
443 longitude_south_pole, latitude_south_pole, angle_rotation, &
444 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
445 latin1, latin2, lad, projection_center_flag, &
446 ellips_smaj_axis, ellips_flatt, ellips_type)
447 TYPE(griddim_def),
INTENT(inout) :: this
448 INTEGER,
INTENT(in),
OPTIONAL :: nx
449 INTEGER,
INTENT(in),
OPTIONAL :: ny
451 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xmin, xmax, ymin, ymax
453 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: dx, dy
456 INTEGER,
INTENT(in),
OPTIONAL :: component_flag
457 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: proj_type
458 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lov
459 INTEGER,
INTENT(in),
OPTIONAL :: zone
460 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xoff
461 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: yoff
462 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_south_pole
463 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_south_pole
464 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: angle_rotation
465 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_stretch_pole
466 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_stretch_pole
467 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: stretch_factor
468 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin1
469 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin2
470 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lad
471 INTEGER,
INTENT(in),
OPTIONAL :: projection_center_flag
472 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_smaj_axis
473 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_flatt
474 INTEGER,
INTENT(in),
OPTIONAL :: ellips_type
476 IF (
PRESENT(nx)) this%dim%nx = nx
477 IF (
PRESENT(ny)) this%dim%ny = ny
479 CALL set_val(this%grid%proj, proj_type=proj_type, lov=lov, zone=zone, &
480 xoff=xoff, yoff=yoff, longitude_south_pole=longitude_south_pole, &
481 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation, &
482 longitude_stretch_pole=longitude_stretch_pole, &
483 latitude_stretch_pole=latitude_stretch_pole, stretch_factor=stretch_factor, &
484 latin1=latin1, latin2=latin2, lad=lad, &
485 projection_center_flag=projection_center_flag, &
486 ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
487 ellips_type=ellips_type)
490 xmin, xmax, ymin, ymax, dx, dy, component_flag)
492 END SUBROUTINE griddim_set_val
499 SUBROUTINE griddim_read_unit(this, unit)
500 TYPE(griddim_def),
INTENT(out) :: this
501 INTEGER,
INTENT(in) :: unit
508 END SUBROUTINE griddim_read_unit
515 SUBROUTINE griddim_write_unit(this, unit)
516 TYPE(griddim_def),
INTENT(in) :: this
517 INTEGER,
INTENT(in) :: unit
524 END SUBROUTINE griddim_write_unit
530 FUNCTION griddim_central_lon(this)
RESULT(lon)
531 TYPE(griddim_def),
INTENT(inout) :: this
533 DOUBLE PRECISION :: lon
535 CALL griddim_pistola_central_lon(this, lon)
537 END FUNCTION griddim_central_lon
543 SUBROUTINE griddim_set_central_lon(this, lonref)
544 TYPE(griddim_def),
INTENT(inout) :: this
545 DOUBLE PRECISION,
INTENT(in) :: lonref
547 DOUBLE PRECISION :: lon
549 CALL griddim_pistola_central_lon(this, lon, lonref)
551 END SUBROUTINE griddim_set_central_lon
555 SUBROUTINE griddim_pistola_central_lon(this, lon, lonref)
556 TYPE(griddim_def),
INTENT(inout) :: this
557 DOUBLE PRECISION,
INTENT(inout) :: lon
558 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lonref
561 DOUBLE PRECISION :: lonsp, latsp, londelta, lov
562 CHARACTER(len=80) :: ptype
565 CALL get_val(this%grid%proj, unit=unit)
566 IF (unit == geo_proj_unit_meter)
THEN 567 CALL get_val(this%grid%proj, lov=lon)
568 IF (
PRESENT(lonref))
THEN 569 CALL long_reset_to_cart_closest(lov, lonref)
570 CALL set_val(this%grid%proj, lov=lon)
573 ELSE IF (unit == geo_proj_unit_degree)
THEN 574 CALL get_val(this%grid%proj, proj_type=ptype, &
575 longitude_south_pole=lonsp, latitude_south_pole=latsp)
577 CASE(
'rotated_ll',
'stretched_rotated_ll')
578 IF (latsp < 0.0d0)
THEN 580 IF (
PRESENT(lonref))
THEN 581 CALL long_reset_to_cart_closest(lov, lonref)
582 CALL set_val(this%grid%proj, longitude_south_pole=lonref)
585 lon = modulo(lonsp + 180.0d0, 360.0d0)
592 IF (
c_e(this%grid%grid%xmin) .AND.
c_e(this%grid%grid%xmin))
THEN 593 lon = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
595 IF (
PRESENT(lonref))
THEN 597 CALL long_reset_to_cart_closest(londelta, lonref)
598 londelta = londelta - lon
599 this%grid%grid%xmin = this%grid%grid%xmin + londelta
600 this%grid%grid%xmax = this%grid%grid%xmax + londelta
605 END SUBROUTINE griddim_pistola_central_lon
611 SUBROUTINE griddim_gen_coord(this, x, y)
612 TYPE(griddim_def),
INTENT(in) :: this
613 DOUBLE PRECISION,
INTENT(out) :: x(:,:)
614 DOUBLE PRECISION,
INTENT(out) :: y(:,:)
617 CALL grid_rect_coordinates(this%grid%grid, x, y)
619 END SUBROUTINE griddim_gen_coord
623 SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
624 TYPE(griddim_def),
INTENT(in) :: this
625 INTEGER,
INTENT(in) :: nx
626 INTEGER,
INTENT(in) :: ny
627 DOUBLE PRECISION,
INTENT(out) :: dx
628 DOUBLE PRECISION,
INTENT(out) :: dy
630 CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
632 END SUBROUTINE griddim_steps
636 SUBROUTINE griddim_setsteps(this)
637 TYPE(griddim_def),
INTENT(inout) :: this
639 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
641 END SUBROUTINE griddim_setsteps
646 ELEMENTAL FUNCTION grid_eq(this, that)
RESULT(res)
647 TYPE(grid_def),
INTENT(IN) :: this, that
650 res = this%proj == that%proj .AND. &
651 this%grid == that%grid
656 ELEMENTAL FUNCTION griddim_eq(this, that)
RESULT(res)
660 res = this%grid == that%grid .AND. &
663 END FUNCTION griddim_eq
666 ELEMENTAL FUNCTION grid_ne(this, that)
RESULT(res)
667 TYPE(
grid_def),
INTENT(IN) :: this, that
670 res = .NOT.(this == that)
675 ELEMENTAL FUNCTION griddim_ne(this, that)
RESULT(res)
679 res = .NOT.(this == that)
681 END FUNCTION griddim_ne
689 SUBROUTINE griddim_import_grid_id(this, ingrid_id)
694 TYPE(
grid_id),
INTENT(in) :: ingrid_id
696 #ifdef HAVE_LIBGRIBAPI 700 TYPE(gdalrasterbandh) :: gdalid
704 #ifdef HAVE_LIBGRIBAPI 705 gaid = grid_id_get_gaid(ingrid_id)
706 IF (
c_e(gaid))
CALL griddim_import_gribapi(this, gaid)
709 gdalid = grid_id_get_gdalid(ingrid_id)
710 IF (gdalassociated(gdalid))
CALL griddim_import_gdal(this, gdalid, &
711 grid_id_get_gdal_options(ingrid_id))
714 END SUBROUTINE griddim_import_grid_id
721 SUBROUTINE griddim_export_grid_id(this, outgrid_id)
726 TYPE(
grid_id),
INTENT(inout) :: outgrid_id
728 #ifdef HAVE_LIBGRIBAPI 732 TYPE(gdalrasterbandh) :: gdalid
735 #ifdef HAVE_LIBGRIBAPI 736 gaid = grid_id_get_gaid(outgrid_id)
737 IF (
c_e(gaid))
CALL griddim_export_gribapi(this, gaid)
740 gdalid = grid_id_get_gdalid(outgrid_id)
745 END SUBROUTINE griddim_export_grid_id
748 #ifdef HAVE_LIBGRIBAPI 750 SUBROUTINE griddim_import_gribapi(this, gaid)
752 TYPE(griddim_def),
INTENT(inout) :: this
753 INTEGER,
INTENT(in) :: gaid
755 DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, x1, y1
756 INTEGER :: EditionNumber, iScansNegatively, jScansPositively, zone, datum, &
760 CALL grib_get(gaid,
'typeOfGrid', this%grid%proj%proj_type)
763 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type))
765 CALL grib_get(gaid,
'GRIBEditionNumber',editionnumber)
768 CALL grib_get(gaid,
'Ni', this%dim%nx)
769 CALL grib_get(gaid,
'Nj', this%dim%ny)
770 CALL griddim_import_ellipsoid(this, gaid)
772 CALL grib_get(gaid,
'iScansNegatively',iscansnegatively)
773 CALL grib_get(gaid,
'jScansPositively',jscanspositively)
774 CALL grib_get(gaid,
'uvRelativeToGrid',this%grid%grid%component_flag)
777 CALL grib_get_dmiss(gaid,
'longitudeOfSouthernPoleInDegrees', &
778 this%grid%proj%rotated%longitude_south_pole)
779 CALL grib_get_dmiss(gaid,
'latitudeOfSouthernPoleInDegrees', &
780 this%grid%proj%rotated%latitude_south_pole)
781 CALL grib_get_dmiss(gaid,
'angleOfRotationInDegrees', &
782 this%grid%proj%rotated%angle_rotation)
787 IF (editionnumber == 1)
THEN 788 CALL grib_get_dmiss(gaid,
'longitudeOfStretchingPoleInDegrees', &
789 this%grid%proj%stretched%longitude_stretch_pole)
790 CALL grib_get_dmiss(gaid,
'latitudeOfStretchingPoleInDegrees', &
791 this%grid%proj%stretched%latitude_stretch_pole)
792 CALL grib_get_dmiss(gaid,
'stretchingFactor', &
793 this%grid%proj%stretched%stretch_factor)
794 ELSE IF (editionnumber == 2)
THEN 795 CALL grib_get_dmiss(gaid,
'longitudeOfThePoleOfStretching', &
796 this%grid%proj%stretched%longitude_stretch_pole)
797 CALL grib_get_dmiss(gaid,
'latitudeOfThePoleOfStretching', &
798 this%grid%proj%stretched%latitude_stretch_pole)
799 CALL grib_get_dmiss(gaid,
'stretchingFactor', &
800 this%grid%proj%stretched%stretch_factor)
801 IF (
c_e(this%grid%proj%stretched%stretch_factor)) &
802 this%grid%proj%stretched%stretch_factor = this%grid%proj%stretched%stretch_factor*1.0d-6
806 SELECT CASE (this%grid%proj%proj_type)
809 CASE (
'regular_ll',
'rotated_ll',
'stretched_ll',
'stretched_rotated_ll')
811 CALL grib_get(gaid,
'longitudeOfFirstGridPointInDegrees',lofirst)
812 CALL grib_get(gaid,
'longitudeOfLastGridPointInDegrees',lolast)
813 CALL grib_get(gaid,
'latitudeOfFirstGridPointInDegrees',lafirst)
814 CALL grib_get(gaid,
'latitudeOfLastGridPointInDegrees',lalast)
820 CALL long_reset_0_360(lofirst)
821 CALL long_reset_0_360(lolast)
823 IF (iscansnegatively == 0)
THEN 824 this%grid%grid%xmin = lofirst
825 this%grid%grid%xmax = lolast
827 this%grid%grid%xmax = lofirst
828 this%grid%grid%xmin = lolast
830 IF (jscanspositively == 0)
THEN 831 this%grid%grid%ymax = lafirst
832 this%grid%grid%ymin = lalast
834 this%grid%grid%ymin = lafirst
835 this%grid%grid%ymax = lalast
839 IF (this%grid%grid%xmax < this%grid%grid%xmin) &
840 this%grid%grid%xmin = this%grid%grid%xmin - 360.d0
843 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
846 CASE (
'polar_stereographic',
'lambert',
'albers')
848 CALL grib_get(gaid,
'DxInMetres', this%grid%grid%dx)
849 CALL grib_get(gaid,
'DyInMetres', this%grid%grid%dy)
851 CALL grib_get_dmiss(gaid,
'Latin1InDegrees',this%grid%proj%polar%latin1)
852 CALL grib_get_dmiss(gaid,
'Latin2InDegrees',this%grid%proj%polar%latin2)
855 "griddim_import_gribapi, latin1/2 "// &
856 trim(to_char(this%grid%proj%polar%latin1))//
" "// &
857 trim(to_char(this%grid%proj%polar%latin2)))
860 CALL grib_get(gaid,
'projectionCenterFlag',&
861 this%grid%proj%polar%projection_center_flag, ierr)
862 IF (ierr /= grib_success)
THEN 863 CALL grib_get(gaid,
'projectionCentreFlag',&
864 this%grid%proj%polar%projection_center_flag)
867 IF (iand(this%grid%proj%polar%projection_center_flag,64) == 1)
THEN 869 "griddim_import_gribapi, bi-polar projections not supported")
873 CALL grib_get(gaid,
'LoVInDegrees',this%grid%proj%lov)
876 "griddim_import_gribapi, central meridian "//trim(to_char(this%grid%proj%lov)))
880 IF (editionnumber == 1)
THEN 890 this%grid%proj%polar%lad = this%grid%proj%polar%latin1
891 ELSE IF (editionnumber == 2)
THEN 892 CALL grib_get(gaid,
'LaDInDegrees',this%grid%proj%polar%lad)
896 "griddim_import_gribapi, lad "//trim(to_char(this%grid%proj%polar%lad)))
900 CALL grib_get(gaid,
'longitudeOfFirstGridPointInDegrees',lofirst)
901 CALL grib_get(gaid,
'latitudeOfFirstGridPointInDegrees',lafirst)
902 CALL long_reset_0_360(lofirst)
903 CALL long_reset_to_cart_closest(this%grid%proj%lov, lofirst)
906 "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
908 "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
912 CALL proj(this, lofirst, lafirst, x1, y1)
913 IF (iscansnegatively == 0)
THEN 914 this%grid%grid%xmin = x1
915 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
917 this%grid%grid%xmax = x1
918 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
920 IF (jscanspositively == 0)
THEN 921 this%grid%grid%ymax = y1
922 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
924 this%grid%grid%ymin = y1
925 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
928 this%grid%proj%polar%lon1 = lofirst
929 this%grid%proj%polar%lat1 = lafirst
933 CALL grib_get(gaid,
'zone',zone)
935 CALL grib_get(gaid,
'datum',datum)
937 CALL grib_get(gaid,
'referenceLongitude',reflon)
938 CALL grib_get(gaid,
'falseEasting',this%grid%proj%xoff)
939 CALL grib_get(gaid,
'falseNorthing',this%grid%proj%yoff)
940 CALL set_val(this%grid%proj, zone=zone, lov=reflon/1.0d6)
943 CALL raise_fatal_error()
946 CALL grib_get(gaid,
'eastingOfFirstGridPoint',lofirst)
947 CALL grib_get(gaid,
'eastingOfLastGridPoint',lolast)
948 CALL grib_get(gaid,
'northingOfFirstGridPoint',lafirst)
949 CALL grib_get(gaid,
'northingOfLastGridPoint',lalast)
951 IF (iscansnegatively == 0)
THEN 952 this%grid%grid%xmin = lofirst
953 this%grid%grid%xmax = lolast
955 this%grid%grid%xmax = lofirst
956 this%grid%grid%xmin = lolast
958 IF (jscanspositively == 0)
THEN 959 this%grid%grid%ymax = lafirst
960 this%grid%grid%ymin = lalast
962 this%grid%grid%ymin = lafirst
963 this%grid%grid%ymax = lalast
967 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
971 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type)//
" not supported")
978 SUBROUTINE grib_get_dmiss(gaid, key, value)
979 INTEGER,
INTENT(in) :: gaid
980 CHARACTER(len=*),
INTENT(in) :: key
981 DOUBLE PRECISION,
INTENT(out) :: value
985 CALL grib_get(gaid, key,
value, ierr)
986 IF (ierr /= grib_success)
value = dmiss
988 END SUBROUTINE grib_get_dmiss
990 SUBROUTINE grib_get_imiss(gaid, key, value)
991 INTEGER,
INTENT(in) :: gaid
992 CHARACTER(len=*),
INTENT(in) :: key
993 INTEGER,
INTENT(out) :: value
997 CALL grib_get(gaid, key,
value, ierr)
998 IF (ierr /= grib_success)
value = imiss
1000 END SUBROUTINE grib_get_imiss
1003 SUBROUTINE griddim_import_ellipsoid(this, gaid)
1004 TYPE(griddim_def),
INTENT(inout) :: this
1005 INTEGER,
INTENT(in) :: gaid
1007 INTEGER :: shapeofearth, iv, is
1008 DOUBLE PRECISION :: r1, r2
1010 IF (editionnumber == 2)
THEN 1011 CALL grib_get(gaid,
'shapeOfTheEarth', shapeofearth)
1012 SELECT CASE(shapeofearth)
1014 CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
1016 CALL grib_get(gaid,
'scaleFactorOfRadiusOfSphericalEarth', is)
1017 CALL grib_get(gaid,
'scaledValueOfRadiusOfSphericalEarth', iv)
1018 r1 = dble(iv) / 10**is
1019 CALL set_val(this, ellips_smaj_axis=r1, ellips_flatt=0.0d0)
1021 CALL set_val(this, ellips_smaj_axis=6378160.0d0, ellips_flatt=1.0d0/297.0d0)
1023 CALL grib_get(gaid,
'scaleFactorOfEarthMajorAxis', is)
1024 CALL grib_get(gaid,
'scaledValueOfEarthMajorAxis', iv)
1025 r1 = dble(iv) / 10**is
1026 CALL grib_get(gaid,
'scaleFactorOfEarthMinorAxis', is)
1027 CALL grib_get(gaid,
'scaledValueOfEarthMinorAxis', iv)
1028 r2 = dble(iv) / 10**is
1029 IF (shapeofearth == 3)
THEN 1033 IF (abs(r1) < 1.0d-6)
THEN 1035 'read from grib, going on with spherical Earth but the results may be wrong')
1036 CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
1038 CALL set_val(this, ellips_smaj_axis=r1, ellips_flatt=(r1-r2)/r1)
1041 CALL set_val(this, ellips_type=ellips_grs80)
1043 CALL set_val(this, ellips_type=ellips_wgs84)
1045 CALL set_val(this, ellips_smaj_axis=6371229.0d0, ellips_flatt=0.0d0)
1049 t2c(shapeofearth)//
' not supported in grib2')
1050 CALL raise_fatal_error()
1056 CALL grib_get(gaid,
'earthIsOblate', shapeofearth)
1057 IF (shapeofearth == 0)
THEN 1058 CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
1060 CALL set_val(this, ellips_smaj_axis=6378160.0d0, ellips_flatt=1.0d0/297.0d0)
1065 END SUBROUTINE griddim_import_ellipsoid
1068 END SUBROUTINE griddim_import_gribapi
1072 SUBROUTINE griddim_export_gribapi(this, gaid)
1074 TYPE(griddim_def),
INTENT(in) :: this
1075 INTEGER,
INTENT(inout) :: gaid
1077 INTEGER :: EditionNumber, iScansNegatively, jScansPositively, nv, pvl, zone, ierr
1078 DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, reflon
1079 DOUBLE PRECISION :: sdx, sdy, ratio, tol
1083 CALL grib_get(gaid,
'GRIBEditionNumber',editionnumber)
1084 CALL grib_set(gaid,
'typeOfGrid' ,this%grid%proj%proj_type)
1087 "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
1091 IF (editionnumber == 1)
THEN 1093 ELSE IF (editionnumber == 2)
THEN 1100 CALL griddim_export_ellipsoid(this, gaid)
1101 CALL grib_set(gaid,
'Ni', this%dim%nx)
1102 CALL grib_set(gaid,
'Nj', this%dim%ny)
1103 CALL grib_set(gaid,
'uvRelativeToGrid',this%grid%grid%component_flag)
1105 CALL grib_get(gaid,
'iScansNegatively',iscansnegatively)
1106 CALL grib_get(gaid,
'jScansPositively',jscanspositively)
1111 CALL grib_set_dmiss(gaid,
'longitudeOfSouthernPoleInDegrees', &
1112 this%grid%proj%rotated%longitude_south_pole, 0.0d0)
1113 CALL grib_set_dmiss(gaid,
'latitudeOfSouthernPoleInDegrees', &
1114 this%grid%proj%rotated%latitude_south_pole, -90.0d0)
1115 IF (editionnumber == 1)
THEN 1116 CALL grib_set_dmiss(gaid,
'angleOfRotationInDegrees', &
1117 this%grid%proj%rotated%angle_rotation, 0.0d0)
1118 ELSE IF (editionnumber == 2)
THEN 1119 CALL grib_set_dmiss(gaid,
'angleOfRotationOfProjectionInDegrees', &
1120 this%grid%proj%rotated%angle_rotation, 0.0d0)
1126 IF (editionnumber == 1)
THEN 1127 CALL grib_set_dmiss(gaid,
'longitudeOfStretchingPoleInDegrees', &
1128 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
1129 CALL grib_set_dmiss(gaid,
'latitudeOfStretchingPoleInDegrees', &
1130 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
1131 CALL grib_set_dmiss(gaid,
'stretchingFactor', &
1132 this%grid%proj%stretched%stretch_factor, 1.0d0)
1133 ELSE IF (editionnumber == 2)
THEN 1134 CALL grib_set_dmiss(gaid,
'longitudeOfThePoleOfStretching', &
1135 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
1136 CALL grib_set_dmiss(gaid,
'latitudeOfThePoleOfStretching', &
1137 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
1138 CALL grib_set_dmiss(gaid,
'stretchingFactor', &
1139 this%grid%proj%stretched%stretch_factor, 1.0d6, 1.0d6)
1143 SELECT CASE (this%grid%proj%proj_type)
1146 CASE (
'regular_ll',
'rotated_ll',
'stretched_ll',
'stretched_rotated_ll')
1148 IF (iscansnegatively == 0)
THEN 1149 lofirst = this%grid%grid%xmin
1150 lolast = this%grid%grid%xmax
1152 lofirst = this%grid%grid%xmax
1153 lolast = this%grid%grid%xmin
1155 IF (jscanspositively == 0)
THEN 1156 lafirst = this%grid%grid%ymax
1157 lalast = this%grid%grid%ymin
1159 lafirst = this%grid%grid%ymin
1160 lalast = this%grid%grid%ymax
1164 IF (editionnumber == 1)
THEN 1165 CALL long_reset_m180_360(lofirst)
1166 CALL long_reset_m180_360(lolast)
1167 ELSE IF (editionnumber == 2)
THEN 1168 CALL long_reset_0_360(lofirst)
1169 CALL long_reset_0_360(lolast)
1172 CALL grib_set(gaid,
'longitudeOfFirstGridPointInDegrees',lofirst)
1173 CALL grib_set(gaid,
'longitudeOfLastGridPointInDegrees',lolast)
1174 CALL grib_set(gaid,
'latitudeOfFirstGridPointInDegrees',lafirst)
1175 CALL grib_set(gaid,
'latitudeOfLastGridPointInDegrees',lalast)
1179 sdx = this%grid%grid%dx*ratio
1180 sdy = this%grid%grid%dy*ratio
1183 IF (abs(nint(sdx)/sdx - 1.0d0) > tol .OR. abs(nint(sdy)/sdy - 1.0d0) > tol)
THEN 1185 "griddim_export_gribapi, increments not given: inaccurate!")
1188 "griddim_export_gribapi, dlon relative error: "//&
1189 trim(to_char(abs(nint(sdx)/sdx - 1.0d0)))//
">"//trim(to_char(tol)))
1191 "griddim_export_gribapi, dlat relative error: "//&
1192 trim(to_char(abs(nint(sdy)/sdy - 1.0d0)))//
">"//trim(to_char(tol)))
1194 CALL grib_set_missing(gaid,
'Di')
1195 CALL grib_set_missing(gaid,
'Dj')
1196 CALL grib_set(gaid,
'ijDirectionIncrementGiven',0)
1200 "griddim_export_gribapi, setting increments: "// &
1201 trim(to_char(this%grid%grid%dx))//
' '//trim(to_char(this%grid%grid%dy)))
1203 CALL grib_set(gaid,
'ijDirectionIncrementGiven',1)
1204 CALL grib_set(gaid,
'iDirectionIncrement',nint(this%grid%grid%dx*ratio))
1205 CALL grib_set(gaid,
'jDirectionIncrement',nint(this%grid%grid%dy*ratio))
1212 CASE (
'polar_stereographic',
'lambert',
'albers')
1214 CALL grib_set(gaid,
'DxInMetres', this%grid%grid%dx)
1215 CALL grib_set(gaid,
'DyInMetres', this%grid%grid%dy)
1216 CALL grib_set(gaid,
'ijDirectionIncrementGiven',1)
1218 CALL grib_set_dmiss(gaid,
'Latin1InDegrees',this%grid%proj%polar%latin1)
1219 CALL grib_set_dmiss(gaid,
'Latin2InDegrees',this%grid%proj%polar%latin2)
1221 CALL grib_set(gaid,
'projectionCenterFlag',&
1222 this%grid%proj%polar%projection_center_flag, ierr)
1223 IF (ierr /= grib_success)
THEN 1224 CALL grib_set(gaid,
'projectionCentreFlag',&
1225 this%grid%proj%polar%projection_center_flag)
1230 CALL grib_set(gaid,
'LoVInDegrees',this%grid%proj%lov)
1232 IF (editionnumber == 2)
THEN 1233 CALL grib_set(gaid,
'LaDInDegrees',this%grid%proj%polar%lad)
1237 IF (iscansnegatively == 0)
THEN 1238 IF (jscanspositively == 0)
THEN 1239 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymax, lofirst, lafirst)
1241 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymin, lofirst, lafirst)
1244 IF (jscanspositively == 0)
THEN 1245 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymax, lofirst, lafirst)
1247 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymin, lofirst, lafirst)
1254 IF (editionnumber == 1)
THEN 1255 CALL long_reset_m180_360(lofirst)
1256 ELSE IF (editionnumber == 2)
THEN 1257 CALL long_reset_0_360(lofirst)
1259 CALL grib_set(gaid,
'longitudeOfFirstGridPointInDegrees',lofirst)
1260 CALL grib_set(gaid,
'latitudeOfFirstGridPointInDegrees',lafirst)
1264 CALL grib_set(gaid,
'datum',0)
1265 CALL get_val(this, zone=zone, lov=reflon)
1266 CALL grib_set(gaid,
'referenceLongitude',nint(reflon/1.0d6))
1267 CALL grib_set(gaid,
'falseEasting',this%grid%proj%xoff)
1268 CALL grib_set(gaid,
'falseNorthing',this%grid%proj%yoff)
1270 CALL grib_set(gaid,
'iDirectionIncrement',this%grid%grid%dx)
1271 CALL grib_set(gaid,
'jDirectionIncrement',this%grid%grid%dy)
1275 CALL grib_set(gaid,
'zone',zone)
1277 IF (iscansnegatively == 0)
THEN 1278 lofirst = this%grid%grid%xmin
1279 lolast = this%grid%grid%xmax
1281 lofirst = this%grid%grid%xmax
1282 lolast = this%grid%grid%xmin
1284 IF (jscanspositively == 0)
THEN 1285 lafirst = this%grid%grid%ymax
1286 lalast = this%grid%grid%ymin
1288 lafirst = this%grid%grid%ymin
1289 lalast = this%grid%grid%ymax
1292 CALL grib_set(gaid,
'eastingOfFirstGridPoint',lofirst)
1293 CALL grib_set(gaid,
'eastingOfLastGridPoint',lolast)
1294 CALL grib_set(gaid,
'northingOfFirstGridPoint',lafirst)
1295 CALL grib_set(gaid,
'northingOfLastGridPoint',lalast)
1299 "griddim_export_gribapi, grid type "//trim(this%grid%proj%proj_type)//
" not supported")
1306 IF (editionnumber == 1)
THEN 1308 CALL grib_get(gaid,
"NV",nv)
1310 CALL l4f_category_log(this%category,l4f_debug,
"griddim_export_gribapi, coding "// &
1311 trim(to_char(nv))//
" vertical coordinate parameters")
1317 SELECT CASE (this%grid%proj%proj_type)
1320 CASE (
'polar_stereographic')
1322 CASE (
'rotated_ll',
'stretched_ll',
'lambert',
'albers')
1324 CASE (
'stretched_rotated_ll')
1331 CALL grib_set(gaid,
"pvlLocation",pvl)
1333 CALL l4f_category_log(this%category,l4f_debug,
"griddim_export_gribapi, coding "// &
1334 trim(to_char(pvl))//
" as vertical coordinate parameter location")
1342 SUBROUTINE grib_set_dmiss(gaid, key, val, default, factor)
1343 INTEGER,
INTENT(in) :: gaid
1344 CHARACTER(len=*),
INTENT(in) :: key
1345 DOUBLE PRECISION,
INTENT(in) :: val
1346 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: default
1347 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: factor
1352 IF (
PRESENT(factor))
THEN 1353 CALL grib_set(gaid, key, val*factor, ierr)
1355 CALL grib_set(gaid, key, val, ierr)
1357 ELSE IF (
PRESENT(default))
THEN 1358 CALL grib_set(gaid, key, default, ierr)
1361 END SUBROUTINE grib_set_dmiss
1363 SUBROUTINE grib_set_imiss(gaid, key, value, default)
1364 INTEGER,
INTENT(in) :: gaid
1365 CHARACTER(len=*),
INTENT(in) :: key
1366 INTEGER,
INTENT(in) :: value
1367 INTEGER,
INTENT(in),
OPTIONAL :: default
1371 IF (
c_e(
value))
THEN 1372 CALL grib_set(gaid, key,
value, ierr)
1373 ELSE IF (
PRESENT(default))
THEN 1374 CALL grib_set(gaid, key, default, ierr)
1377 END SUBROUTINE grib_set_imiss
1379 SUBROUTINE griddim_export_ellipsoid(this, gaid)
1380 TYPE(griddim_def),
INTENT(in) :: this
1381 INTEGER,
INTENT(in) :: gaid
1383 INTEGER :: ellips_type
1384 DOUBLE PRECISION :: r1, r2, f
1386 CALL get_val(this, ellips_smaj_axis=r1, ellips_flatt=f, ellips_type=ellips_type)
1388 IF (editionnumber == 2)
THEN 1390 CALL grib_set_missing(gaid,
'scaleFactorOfRadiusOfSphericalEarth')
1391 CALL grib_set_missing(gaid,
'scaledValueOfRadiusOfSphericalEarth')
1392 CALL grib_set_missing(gaid,
'scaleFactorOfEarthMajorAxis')
1393 CALL grib_set_missing(gaid,
'scaledValueOfEarthMajorAxis')
1394 CALL grib_set_missing(gaid,
'scaleFactorOfEarthMinorAxis')
1395 CALL grib_set_missing(gaid,
'scaledValueOfEarthMinorAxis')
1397 SELECT CASE(ellips_type)
1399 CALL grib_set(gaid,
'shapeOfTheEarth', 4)
1401 CALL grib_set(gaid,
'shapeOfTheEarth', 5)
1403 IF (f == 0.0d0)
THEN 1404 IF (r1 == 6367470.0d0)
THEN 1405 CALL grib_set(gaid,
'shapeOfTheEarth', 0)
1406 ELSE IF (r1 == 6371229.0d0)
THEN 1407 CALL grib_set(gaid,
'shapeOfTheEarth', 6)
1409 CALL grib_set(gaid,
'shapeOfTheEarth', 1)
1410 CALL grib_set(gaid,
'scaleFactorOfRadiusOfSphericalEarth', 2)
1411 CALL grib_set(gaid,
'scaledValueOfRadiusOfSphericalEarth', int(r1*100.0d0))
1414 IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0)
THEN 1415 CALL grib_set(gaid,
'shapeOfTheEarth', 2)
1417 CALL grib_set(gaid,
'shapeOfTheEarth', 3)
1419 CALL grib_set(gaid,
'scaleFactorOfEarthMajorAxis', 5)
1420 CALL grib_set(gaid,
'scaledValueOfEarthMajorAxis', &
1422 CALL grib_set(gaid,
'scaleFactorOfEarthMinorAxis', 5)
1423 CALL grib_set(gaid,
'scaledValueOfEarthMinorAxis', &
1431 IF (r1 == 6367470.0d0 .AND. f == 0.0d0)
THEN 1432 CALL grib_set(gaid,
'earthIsOblate', 0)
1433 ELSE IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0)
THEN 1434 CALL grib_set(gaid,
'earthIsOblate', 1)
1435 ELSE IF (f == 0.0d0)
THEN 1436 CALL grib_set(gaid,
'earthIsOblate', 0)
1437 CALL l4f_category_log(this%category,l4f_warn,
'desired spherical Earth radius & 1438 ¬ supported in grib 1, coding with default radius of 6367470 m')
1440 CALL grib_set(gaid,
'earthIsOblate', 1)
1442 ¬ supported in grib 1, coding with default iau65 ellipsoid')
1447 END SUBROUTINE griddim_export_ellipsoid
1449 END SUBROUTINE griddim_export_gribapi
1455 SUBROUTINE griddim_import_gdal(this, gdalid, gdal_options)
1457 TYPE(griddim_def),
INTENT(inout) :: this
1458 TYPE(gdalrasterbandh),
INTENT(in) :: gdalid
1459 TYPE(gdal_file_id_options),
INTENT(in) :: gdal_options
1461 TYPE(gdaldataseth) :: hds
1462 REAL(kind=c_double) :: geotrans(6)
1463 INTEGER(kind=c_int) :: offsetx, offsety
1466 hds = gdalgetbanddataset(gdalid)
1467 ier = gdalgetgeotransform(hds, geotrans)
1471 'griddim_import_gdal: error in accessing gdal dataset')
1475 IF (geotrans(3) /= 0.0_c_double .OR. geotrans(5) /= 0.0_c_double)
THEN 1477 'griddim_import_gdal: dataset has a non-diagonal transformation matrix, unsupported')
1482 CALL gdaldatasetbbsize_f(hds, gdal_options%xmin, gdal_options%ymin, &
1483 gdal_options%xmax, gdal_options%ymax, &
1484 this%dim%nx, this%dim%ny, offsetx, offsety, &
1485 this%grid%grid%xmin, this%grid%grid%ymin, this%grid%grid%xmax, this%grid%grid%ymax)
1487 IF (this%dim%nx == 0 .OR. this%dim%ny == 0)
THEN 1489 'griddim_import_gdal: requested bounding box '//t2c(gdal_options%xmin)//
','// &
1490 t2c(gdal_options%ymin)//
','//t2c(gdal_options%xmax)//
','//&
1491 t2c(gdal_options%ymax))
1493 'determines an empty dataset '//t2c(this%dim%nx)//
'x'//t2c(this%dim%ny))
1520 this%grid%proj%proj_type =
'regular_ll' 1522 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
1523 this%grid%grid%component_flag = 0
1525 END SUBROUTINE griddim_import_gdal
1530 SUBROUTINE griddim_display(this)
1531 TYPE(griddim_def),
INTENT(in) :: this
1533 print*,
"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>" 1539 print*,
"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>" 1541 END SUBROUTINE griddim_display
1544 #define VOL7D_POLY_TYPE TYPE(grid_def) 1545 #define VOL7D_POLY_TYPES _grid 1546 #include "array_utilities_inc.F90" 1547 #undef VOL7D_POLY_TYPE 1548 #undef VOL7D_POLY_TYPES 1550 #define VOL7D_POLY_TYPE TYPE(griddim_def) 1551 #define VOL7D_POLY_TYPES _griddim 1552 #include "array_utilities_inc.F90" 1553 #undef VOL7D_POLY_TYPE 1554 #undef VOL7D_POLY_TYPES 1568 SUBROUTINE griddim_wind_unrot(this, rot_mat)
1570 TYPE(griddim_def),
INTENT(in) :: this
1571 DOUBLE PRECISION,
POINTER :: rot_mat(:,:,:)
1573 DOUBLE PRECISION :: dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j
1574 DOUBLE PRECISION :: lat_factor
1575 INTEGER :: i, j, ip1, im1, jp1, jm1;
1577 IF (this%dim%nx <= 0 .OR. this%dim%ny <= 0 .OR. &
1578 .NOT.
ASSOCIATED(this%dim%lon) .OR. .NOT.
ASSOCIATED(this%dim%lat))
THEN 1579 CALL l4f_category_log(this%category, l4f_error,
'rotate_uv must be called after correct init of griddim object')
1584 ALLOCATE(rot_mat(this%dim%nx, this%dim%ny, 4))
1586 DO j = 1, this%dim%ny
1587 jp1 = min(j+1, this%dim%ny)
1589 DO i = 1, this%dim%nx
1590 ip1 = min(i+1, this%dim%nx)
1593 dlat_i = this%dim%lat(ip1,j) - this%dim%lat(im1,j)
1594 dlat_j = this%dim%lat(i,jp1) - this%dim%lat(i,jm1)
1596 dlon_i = this%dim%lon(ip1,j) - this%dim%lon(im1,j)
1599 dlon_j = this%dim%lon(i,jp1) - this%dim%lon(i,jm1)
1604 lat_factor = cos(degrad*this%dim%lat(i,j))
1605 dlon_i = dlon_i * lat_factor
1606 dlon_j = dlon_j * lat_factor
1608 dist_i = sqrt(dlon_i*dlon_i + dlat_i*dlat_i)
1609 dist_j = sqrt(dlon_j*dlon_j + dlat_j*dlat_j)
1611 IF (dist_i > 0.d0)
THEN 1612 rot_mat(i,j,1) = dlon_i/dist_i
1613 rot_mat(i,j,2) = dlat_i/dist_i
1615 rot_mat(i,j,1) = 0.0d0
1616 rot_mat(i,j,2) = 0.0d0
1618 IF (dist_j > 0.d0)
THEN 1619 rot_mat(i,j,3) = dlon_j/dist_j
1620 rot_mat(i,j,4) = dlat_j/dist_j
1622 rot_mat(i,j,3) = 0.0d0
1623 rot_mat(i,j,4) = 0.0d0
1629 END SUBROUTINE griddim_wind_unrot
1633 SUBROUTINE griddim_zoom_coord(this, ilon, ilat, flon, flat, ix, iy, fx, fy)
1634 TYPE(griddim_def),
INTENT(in) :: this
1635 DOUBLE PRECISION,
INTENT(in) :: ilon, ilat, flon, flat
1636 INTEGER,
INTENT(out) :: ix, iy, fx, fy
1638 DOUBLE PRECISION :: ix1, iy1, fx1, fy1
1641 CALL proj(this, ilon, ilat, ix1, iy1)
1642 CALL proj(this, flon, flat, fx1, fy1)
1644 CALL griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
1646 END SUBROUTINE griddim_zoom_coord
1650 SUBROUTINE griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
1651 TYPE(griddim_def),
INTENT(in) :: this
1652 DOUBLE PRECISION,
INTENT(in) :: ix1, iy1, fx1, fy1
1653 INTEGER,
INTENT(out) :: ix, iy, fx, fy
1655 INTEGER :: lix, liy, lfx, lfy
1658 lix = nint((ix1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
1659 liy = nint((iy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
1660 lfx = nint((fx1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
1661 lfy = nint((fy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
1668 END SUBROUTINE griddim_zoom_projcoord
1674 SUBROUTINE long_reset_0_360(lon)
1675 DOUBLE PRECISION,
INTENT(inout) :: lon
1677 IF (.NOT.
c_e(lon))
RETURN 1678 DO WHILE(lon < 0.0d0)
1681 DO WHILE(lon >= 360.0d0)
1685 END SUBROUTINE long_reset_0_360
1691 SUBROUTINE long_reset_m180_360(lon)
1692 DOUBLE PRECISION,
INTENT(inout) :: lon
1694 IF (.NOT.
c_e(lon))
RETURN 1695 DO WHILE(lon < -180.0d0)
1698 DO WHILE(lon >= 360.0d0)
1702 END SUBROUTINE long_reset_m180_360
1725 SUBROUTINE long_reset_m180_180(lon)
1726 DOUBLE PRECISION,
INTENT(inout) :: lon
1728 IF (.NOT.
c_e(lon))
RETURN 1729 DO WHILE(lon < -180.0d0)
1732 DO WHILE(lon >= 180.0d0)
1736 END SUBROUTINE long_reset_m180_180
1739 SUBROUTINE long_reset_to_cart_closest(lon, lonref)
1740 DOUBLE PRECISION,
INTENT(inout) :: lon
1741 DOUBLE PRECISION,
INTENT(in) :: lonref
1743 IF (.NOT.
c_e(lon) .OR. .NOT.
c_e(lonref))
RETURN 1744 IF (abs(lon-lonref) < 180.0d0)
RETURN 1745 lon = lon - nint((lon-lonref)/360.0d0)*360.0d0
1747 END SUBROUTINE long_reset_to_cart_closest
Export griddim object to grid_id.
Method for setting the contents of the object.
This object, mainly for internal use, describes a grid on a geographical projection, except the grid dimensions.
Compute forward coordinate transformation from geographical system to projected system.
Constructors of the corresponding objects.
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Print a brief description on stdout.
This module defines an abstract interface to different drivers for access to files containing gridded...
Module for describing geographically referenced regular grids.
Write the object on a formatted or unformatted file.
Compute backward coordinate transformation from projected system to geographical system.
This object completely describes a grid on a geographic projection.
Method for returning the contents of the object.
Definitions of constants and functions for working with missing values.
Costanti fisiche (DOUBLEPRECISION).
Module for defining the extension and coordinates of a rectangular georeferenced grid.
classe per la gestione del logging
Method for testing the existence of the object.
Copy an object, creating a fully new instance.
Read the object from a formatted or unformatted file.
Emit log message for a category with specific priority.
Derived type describing the extension of a grid and the geographical coordinates of each point...
Destructors of the corresponding objects.
Import griddim object from grid_id.