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)
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, lonrot
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(lon, lonref)
582 CALL set_val(this%grid%proj, longitude_south_pole=lonref)
584 IF (
c_e(this%grid%grid%xmin) .AND.
c_e(this%grid%grid%xmax))
THEN 585 lonrot = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
588 CALL long_reset_to_cart_closest(londelta, 0.0d0)
589 londelta = londelta - lonrot
590 this%grid%grid%xmin = this%grid%grid%xmin + londelta
591 this%grid%grid%xmax = this%grid%grid%xmax + londelta
594 lon = modulo(lonsp + 180.0d0, 360.0d0)
601 IF (
c_e(this%grid%grid%xmin) .AND.
c_e(this%grid%grid%xmax))
THEN 602 lon = 0.5d0*(this%grid%grid%xmin + this%grid%grid%xmax)
604 IF (
PRESENT(lonref))
THEN 606 CALL long_reset_to_cart_closest(londelta, lonref)
607 londelta = londelta - lon
608 this%grid%grid%xmin = this%grid%grid%xmin + londelta
609 this%grid%grid%xmax = this%grid%grid%xmax + londelta
614 END SUBROUTINE griddim_pistola_central_lon
620 SUBROUTINE griddim_gen_coord(this, x, y)
622 DOUBLE PRECISION,
INTENT(out) :: x(:,:)
623 DOUBLE PRECISION,
INTENT(out) :: y(:,:)
626 CALL grid_rect_coordinates(this%grid%grid, x, y)
628 END SUBROUTINE griddim_gen_coord
632 SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
634 INTEGER,
INTENT(in) :: nx
635 INTEGER,
INTENT(in) :: ny
636 DOUBLE PRECISION,
INTENT(out) :: dx
637 DOUBLE PRECISION,
INTENT(out) :: dy
639 CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
641 END SUBROUTINE griddim_steps
645 SUBROUTINE griddim_setsteps(this)
648 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
650 END SUBROUTINE griddim_setsteps
655 ELEMENTAL FUNCTION grid_eq(this, that)
RESULT(res)
656 TYPE(
grid_def),
INTENT(IN) :: this, that
659 res = this%proj == that%proj .AND. &
660 this%grid == that%grid
665 ELEMENTAL FUNCTION griddim_eq(this, that)
RESULT(res)
669 res = this%grid == that%grid .AND. &
672 END FUNCTION griddim_eq
675 ELEMENTAL FUNCTION grid_ne(this, that)
RESULT(res)
676 TYPE(grid_def),
INTENT(IN) :: this, that
679 res = .NOT.(this == that)
684 ELEMENTAL FUNCTION griddim_ne(this, that)
RESULT(res)
685 TYPE(griddim_def),
INTENT(IN) :: this, that
688 res = .NOT.(this == that)
690 END FUNCTION griddim_ne
698 SUBROUTINE griddim_import_grid_id(this, ingrid_id)
702 TYPE(griddim_def),
INTENT(inout) :: this
703 TYPE(grid_id),
INTENT(in) :: ingrid_id
705 #ifdef HAVE_LIBGRIBAPI 709 TYPE(gdalrasterbandh) :: gdalid
713 #ifdef HAVE_LIBGRIBAPI 714 gaid = grid_id_get_gaid(ingrid_id)
715 IF (
c_e(gaid))
CALL griddim_import_gribapi(this, gaid)
718 gdalid = grid_id_get_gdalid(ingrid_id)
719 IF (gdalassociated(gdalid))
CALL griddim_import_gdal(this, gdalid, &
720 grid_id_get_gdal_options(ingrid_id))
723 END SUBROUTINE griddim_import_grid_id
730 SUBROUTINE griddim_export_grid_id(this, outgrid_id)
734 TYPE(griddim_def),
INTENT(in) :: this
735 TYPE(grid_id),
INTENT(inout) :: outgrid_id
737 #ifdef HAVE_LIBGRIBAPI 741 TYPE(gdalrasterbandh) :: gdalid
744 #ifdef HAVE_LIBGRIBAPI 745 gaid = grid_id_get_gaid(outgrid_id)
746 IF (
c_e(gaid))
CALL griddim_export_gribapi(this, gaid)
749 gdalid = grid_id_get_gdalid(outgrid_id)
754 END SUBROUTINE griddim_export_grid_id
757 #ifdef HAVE_LIBGRIBAPI 759 SUBROUTINE griddim_import_gribapi(this, gaid)
761 TYPE(griddim_def),
INTENT(inout) :: this
762 INTEGER,
INTENT(in) :: gaid
764 DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, x1, y1, orient
765 INTEGER :: EditionNumber, iScansNegatively, jScansPositively, zone, datum, &
769 CALL grib_get(gaid,
'typeOfGrid', this%grid%proj%proj_type)
772 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type))
774 CALL grib_get(gaid,
'GRIBEditionNumber',editionnumber)
776 IF (this%grid%proj%proj_type ==
'unstructured_grid')
THEN 777 CALL grib_get(gaid,
'numberOfDataPoints', this%dim%nx)
779 this%grid%grid%component_flag = 0
780 CALL griddim_import_ellipsoid(this, gaid)
785 CALL grib_get(gaid,
'Ni', this%dim%nx)
786 CALL grib_get(gaid,
'Nj', this%dim%ny)
787 CALL griddim_import_ellipsoid(this, gaid)
789 CALL grib_get(gaid,
'iScansNegatively',iscansnegatively)
790 CALL grib_get(gaid,
'jScansPositively',jscanspositively)
791 CALL grib_get(gaid,
'uvRelativeToGrid',this%grid%grid%component_flag)
794 CALL grib_get_dmiss(gaid,
'longitudeOfSouthernPoleInDegrees', &
795 this%grid%proj%rotated%longitude_south_pole)
796 CALL grib_get_dmiss(gaid,
'latitudeOfSouthernPoleInDegrees', &
797 this%grid%proj%rotated%latitude_south_pole)
798 CALL grib_get_dmiss(gaid,
'angleOfRotationInDegrees', &
799 this%grid%proj%rotated%angle_rotation)
804 IF (editionnumber == 1)
THEN 805 CALL grib_get_dmiss(gaid,
'longitudeOfStretchingPoleInDegrees', &
806 this%grid%proj%stretched%longitude_stretch_pole)
807 CALL grib_get_dmiss(gaid,
'latitudeOfStretchingPoleInDegrees', &
808 this%grid%proj%stretched%latitude_stretch_pole)
809 CALL grib_get_dmiss(gaid,
'stretchingFactor', &
810 this%grid%proj%stretched%stretch_factor)
811 ELSE IF (editionnumber == 2)
THEN 812 CALL grib_get_dmiss(gaid,
'longitudeOfThePoleOfStretching', &
813 this%grid%proj%stretched%longitude_stretch_pole)
814 CALL grib_get_dmiss(gaid,
'latitudeOfThePoleOfStretching', &
815 this%grid%proj%stretched%latitude_stretch_pole)
816 CALL grib_get_dmiss(gaid,
'stretchingFactor', &
817 this%grid%proj%stretched%stretch_factor)
818 IF (
c_e(this%grid%proj%stretched%stretch_factor)) &
819 this%grid%proj%stretched%stretch_factor = this%grid%proj%stretched%stretch_factor*1.0d-6
823 SELECT CASE (this%grid%proj%proj_type)
826 CASE (
'regular_ll',
'rotated_ll',
'stretched_ll',
'stretched_rotated_ll')
828 CALL grib_get(gaid,
'longitudeOfFirstGridPointInDegrees',lofirst)
829 CALL grib_get(gaid,
'longitudeOfLastGridPointInDegrees',lolast)
830 CALL grib_get(gaid,
'latitudeOfFirstGridPointInDegrees',lafirst)
831 CALL grib_get(gaid,
'latitudeOfLastGridPointInDegrees',lalast)
837 CALL long_reset_0_360(lofirst)
838 CALL long_reset_0_360(lolast)
840 IF (iscansnegatively == 0)
THEN 841 this%grid%grid%xmin = lofirst
842 this%grid%grid%xmax = lolast
844 this%grid%grid%xmax = lofirst
845 this%grid%grid%xmin = lolast
847 IF (jscanspositively == 0)
THEN 848 this%grid%grid%ymax = lafirst
849 this%grid%grid%ymin = lalast
851 this%grid%grid%ymin = lafirst
852 this%grid%grid%ymax = lalast
856 IF (this%grid%grid%xmax < this%grid%grid%xmin) &
857 this%grid%grid%xmin = this%grid%grid%xmin - 360.d0
860 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
863 CASE (
'polar_stereographic',
'lambert',
'albers')
865 CALL grib_get(gaid,
'DxInMetres', this%grid%grid%dx)
866 CALL grib_get(gaid,
'DyInMetres', this%grid%grid%dy)
868 CALL grib_get_dmiss(gaid,
'Latin1InDegrees',this%grid%proj%polar%latin1)
869 CALL grib_get_dmiss(gaid,
'Latin2InDegrees',this%grid%proj%polar%latin2)
872 "griddim_import_gribapi, latin1/2 "// &
873 trim(to_char(this%grid%proj%polar%latin1))//
" "// &
874 trim(to_char(this%grid%proj%polar%latin2)))
877 CALL grib_get(gaid,
'projectionCenterFlag',&
878 this%grid%proj%polar%projection_center_flag, ierr)
879 IF (ierr /= grib_success)
THEN 880 CALL grib_get(gaid,
'projectionCentreFlag',&
881 this%grid%proj%polar%projection_center_flag)
884 IF (iand(this%grid%proj%polar%projection_center_flag,64) == 1)
THEN 886 "griddim_import_gribapi, bi-polar projections not supported")
890 CALL grib_get(gaid,
'LoVInDegrees',this%grid%proj%lov)
893 "griddim_import_gribapi, central meridian "//trim(to_char(this%grid%proj%lov)))
897 IF (editionnumber == 1)
THEN 907 this%grid%proj%polar%lad = this%grid%proj%polar%latin1
908 ELSE IF (editionnumber == 2)
THEN 909 CALL grib_get(gaid,
'LaDInDegrees',this%grid%proj%polar%lad)
913 "griddim_import_gribapi, lad "//trim(to_char(this%grid%proj%polar%lad)))
917 CALL grib_get(gaid,
'longitudeOfFirstGridPointInDegrees',lofirst)
918 CALL grib_get(gaid,
'latitudeOfFirstGridPointInDegrees',lafirst)
919 CALL long_reset_0_360(lofirst)
920 CALL long_reset_to_cart_closest(this%grid%proj%lov, lofirst)
923 "griddim_import_gribapi, longitude of first point "//trim(to_char(lofirst)))
925 "griddim_import_gribapi, central meridian reset "//trim(to_char(this%grid%proj%lov)))
928 CALL proj(this, lofirst, lafirst, x1, y1)
929 IF (iscansnegatively == 0)
THEN 930 this%grid%grid%xmin = x1
931 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
933 this%grid%grid%xmax = x1
934 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
936 IF (jscanspositively == 0)
THEN 937 this%grid%grid%ymax = y1
938 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
940 this%grid%grid%ymin = y1
941 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
944 this%grid%proj%polar%lon1 = lofirst
945 this%grid%proj%polar%lat1 = lafirst
950 CALL grib_get(gaid,
'DxInMetres', this%grid%grid%dx)
951 CALL grib_get(gaid,
'DyInMetres', this%grid%grid%dy)
952 CALL grib_get(gaid,
'LaDInDegrees',this%grid%proj%polar%lad)
953 this%grid%proj%lov = 0.0d0
955 CALL grib_get(gaid,
'longitudeOfFirstGridPointInDegrees',lofirst)
956 CALL grib_get(gaid,
'latitudeOfFirstGridPointInDegrees',lafirst)
958 CALL proj(this, lofirst, lafirst, x1, y1)
959 IF (iscansnegatively == 0)
THEN 960 this%grid%grid%xmin = x1
961 this%grid%grid%xmax = x1 + this%grid%grid%dx*dble(this%dim%nx - 1)
963 this%grid%grid%xmax = x1
964 this%grid%grid%xmin = x1 - this%grid%grid%dx*dble(this%dim%nx - 1)
966 IF (jscanspositively == 0)
THEN 967 this%grid%grid%ymax = y1
968 this%grid%grid%ymin = y1 - this%grid%grid%dy*dble(this%dim%ny - 1)
970 this%grid%grid%ymin = y1
971 this%grid%grid%ymax = y1 + this%grid%grid%dy*dble(this%dim%ny - 1)
974 IF (editionnumber == 2)
THEN 975 CALL grib_get(gaid,
'orientationOfTheGridInDegrees',orient)
976 IF (orient /= 0.0d0)
THEN 978 "griddim_import_gribapi, Mercator grid orientation != 0 not supported")
984 CALL unproj(this, x1, y1, lofirst, lafirst)
986 "griddim_import_gribapi, unprojected first point "//t2c(lofirst)//
" "// &
989 CALL grib_get(gaid,
'longitudeOfLastGridPointInDegrees',lolast)
990 CALL grib_get(gaid,
'latitudeOfLastGridPointInDegrees',lalast)
991 CALL proj(this, lolast, lalast, x1, y1)
993 "griddim_import_gribapi, extremes from grib "//t2c(x1)//
" "//t2c(y1))
995 "griddim_import_gribapi, extremes from proj "//t2c(this%grid%grid%xmin)// &
996 " "//t2c(this%grid%grid%ymin)//
" "//t2c(this%grid%grid%xmax)//
" "// &
997 t2c(this%grid%grid%ymax))
1002 CALL grib_get(gaid,
'zone',zone)
1004 CALL grib_get(gaid,
'datum',datum)
1005 IF (datum == 0)
THEN 1006 CALL grib_get(gaid,
'referenceLongitude',reflon)
1007 CALL grib_get(gaid,
'falseEasting',this%grid%proj%xoff)
1008 CALL grib_get(gaid,
'falseNorthing',this%grid%proj%yoff)
1009 CALL set_val(this%grid%proj, zone=zone, lov=reflon/1.0d6)
1012 CALL raise_fatal_error()
1015 CALL grib_get(gaid,
'eastingOfFirstGridPoint',lofirst)
1016 CALL grib_get(gaid,
'eastingOfLastGridPoint',lolast)
1017 CALL grib_get(gaid,
'northingOfFirstGridPoint',lafirst)
1018 CALL grib_get(gaid,
'northingOfLastGridPoint',lalast)
1020 IF (iscansnegatively == 0)
THEN 1021 this%grid%grid%xmin = lofirst
1022 this%grid%grid%xmax = lolast
1024 this%grid%grid%xmax = lofirst
1025 this%grid%grid%xmin = lolast
1027 IF (jscanspositively == 0)
THEN 1028 this%grid%grid%ymax = lafirst
1029 this%grid%grid%ymin = lalast
1031 this%grid%grid%ymin = lafirst
1032 this%grid%grid%ymax = lalast
1036 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
1040 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type)//
" not supported")
1047 SUBROUTINE grib_get_dmiss(gaid, key, value)
1048 INTEGER,
INTENT(in) :: gaid
1049 CHARACTER(len=*),
INTENT(in) :: key
1050 DOUBLE PRECISION,
INTENT(out) :: value
1054 CALL grib_get(gaid, key,
value, ierr)
1055 IF (ierr /= grib_success)
value = dmiss
1057 END SUBROUTINE grib_get_dmiss
1059 SUBROUTINE grib_get_imiss(gaid, key, value)
1060 INTEGER,
INTENT(in) :: gaid
1061 CHARACTER(len=*),
INTENT(in) :: key
1062 INTEGER,
INTENT(out) :: value
1066 CALL grib_get(gaid, key,
value, ierr)
1067 IF (ierr /= grib_success)
value = imiss
1069 END SUBROUTINE grib_get_imiss
1072 SUBROUTINE griddim_import_ellipsoid(this, gaid)
1073 TYPE(griddim_def),
INTENT(inout) :: this
1074 INTEGER,
INTENT(in) :: gaid
1076 INTEGER :: shapeofearth, iv, is
1077 DOUBLE PRECISION :: r1, r2
1079 IF (editionnumber == 2)
THEN 1080 CALL grib_get(gaid,
'shapeOfTheEarth', shapeofearth)
1081 SELECT CASE(shapeofearth)
1083 CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
1085 CALL grib_get(gaid,
'scaleFactorOfRadiusOfSphericalEarth', is)
1086 CALL grib_get(gaid,
'scaledValueOfRadiusOfSphericalEarth', iv)
1087 r1 = dble(iv) / 10**is
1088 CALL set_val(this, ellips_smaj_axis=r1, ellips_flatt=0.0d0)
1090 CALL set_val(this, ellips_smaj_axis=6378160.0d0, ellips_flatt=1.0d0/297.0d0)
1092 CALL grib_get(gaid,
'scaleFactorOfEarthMajorAxis', is)
1093 CALL grib_get(gaid,
'scaledValueOfEarthMajorAxis', iv)
1094 r1 = dble(iv) / 10**is
1095 CALL grib_get(gaid,
'scaleFactorOfEarthMinorAxis', is)
1096 CALL grib_get(gaid,
'scaledValueOfEarthMinorAxis', iv)
1097 r2 = dble(iv) / 10**is
1098 IF (shapeofearth == 3)
THEN 1102 IF (abs(r1) < 1.0d-6)
THEN 1104 'read from grib, going on with spherical Earth but the results may be wrong')
1105 CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
1107 CALL set_val(this, ellips_smaj_axis=r1, ellips_flatt=(r1-r2)/r1)
1110 CALL set_val(this, ellips_type=ellips_grs80)
1112 CALL set_val(this, ellips_type=ellips_wgs84)
1114 CALL set_val(this, ellips_smaj_axis=6371229.0d0, ellips_flatt=0.0d0)
1118 t2c(shapeofearth)//
' not supported in grib2')
1119 CALL raise_fatal_error()
1125 CALL grib_get(gaid,
'earthIsOblate', shapeofearth)
1126 IF (shapeofearth == 0)
THEN 1127 CALL set_val(this, ellips_smaj_axis=6367470.0d0, ellips_flatt=0.0d0)
1129 CALL set_val(this, ellips_smaj_axis=6378160.0d0, ellips_flatt=1.0d0/297.0d0)
1134 END SUBROUTINE griddim_import_ellipsoid
1137 END SUBROUTINE griddim_import_gribapi
1141 SUBROUTINE griddim_export_gribapi(this, gaid)
1143 TYPE(griddim_def),
INTENT(in) :: this
1144 INTEGER,
INTENT(inout) :: gaid
1146 INTEGER :: EditionNumber, iScansNegatively, jScansPositively, nv, pvl, zone, ierr
1147 DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, reflon
1148 DOUBLE PRECISION :: sdx, sdy, ratio, tol, ex, ey
1152 CALL grib_get(gaid,
'GRIBEditionNumber', editionnumber)
1154 IF (editionnumber == 2)
CALL grib_set(gaid,
'shapeOfTheEarth', 6)
1155 CALL grib_set(gaid,
'typeOfGrid', this%grid%proj%proj_type)
1158 "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
1161 IF (this%grid%proj%proj_type ==
'unstructured_grid')
THEN 1162 CALL grib_set(gaid,
'numberOfDataPoints', this%dim%nx)
1170 IF (editionnumber == 1)
THEN 1172 ELSE IF (editionnumber == 2)
THEN 1179 CALL griddim_export_ellipsoid(this, gaid)
1180 CALL grib_set(gaid,
'Ni', this%dim%nx)
1181 CALL grib_set(gaid,
'Nj', this%dim%ny)
1182 CALL grib_set(gaid,
'uvRelativeToGrid',this%grid%grid%component_flag)
1184 CALL grib_get(gaid,
'iScansNegatively',iscansnegatively)
1185 CALL grib_get(gaid,
'jScansPositively',jscanspositively)
1190 CALL grib_set_dmiss(gaid,
'longitudeOfSouthernPoleInDegrees', &
1191 this%grid%proj%rotated%longitude_south_pole, 0.0d0)
1192 CALL grib_set_dmiss(gaid,
'latitudeOfSouthernPoleInDegrees', &
1193 this%grid%proj%rotated%latitude_south_pole, -90.0d0)
1194 IF (editionnumber == 1)
THEN 1195 CALL grib_set_dmiss(gaid,
'angleOfRotationInDegrees', &
1196 this%grid%proj%rotated%angle_rotation, 0.0d0)
1197 ELSE IF (editionnumber == 2)
THEN 1198 CALL grib_set_dmiss(gaid,
'angleOfRotationOfProjectionInDegrees', &
1199 this%grid%proj%rotated%angle_rotation, 0.0d0)
1205 IF (editionnumber == 1)
THEN 1206 CALL grib_set_dmiss(gaid,
'longitudeOfStretchingPoleInDegrees', &
1207 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
1208 CALL grib_set_dmiss(gaid,
'latitudeOfStretchingPoleInDegrees', &
1209 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
1210 CALL grib_set_dmiss(gaid,
'stretchingFactor', &
1211 this%grid%proj%stretched%stretch_factor, 1.0d0)
1212 ELSE IF (editionnumber == 2)
THEN 1213 CALL grib_set_dmiss(gaid,
'longitudeOfThePoleOfStretching', &
1214 this%grid%proj%stretched%longitude_stretch_pole, 0.0d0)
1215 CALL grib_set_dmiss(gaid,
'latitudeOfThePoleOfStretching', &
1216 this%grid%proj%stretched%latitude_stretch_pole, -90.0d0)
1217 CALL grib_set_dmiss(gaid,
'stretchingFactor', &
1218 this%grid%proj%stretched%stretch_factor, 1.0d6, 1.0d6)
1222 SELECT CASE (this%grid%proj%proj_type)
1225 CASE (
'regular_ll',
'rotated_ll',
'stretched_ll',
'stretched_rotated_ll')
1227 IF (iscansnegatively == 0)
THEN 1228 lofirst = this%grid%grid%xmin
1229 lolast = this%grid%grid%xmax
1231 lofirst = this%grid%grid%xmax
1232 lolast = this%grid%grid%xmin
1234 IF (jscanspositively == 0)
THEN 1235 lafirst = this%grid%grid%ymax
1236 lalast = this%grid%grid%ymin
1238 lafirst = this%grid%grid%ymin
1239 lalast = this%grid%grid%ymax
1243 IF (editionnumber == 1)
THEN 1244 CALL long_reset_m180_360(lofirst)
1245 CALL long_reset_m180_360(lolast)
1246 ELSE IF (editionnumber == 2)
THEN 1247 CALL long_reset_0_360(lofirst)
1248 CALL long_reset_0_360(lolast)
1251 CALL grib_set(gaid,
'longitudeOfFirstGridPointInDegrees',lofirst)
1252 CALL grib_set(gaid,
'longitudeOfLastGridPointInDegrees',lolast)
1253 CALL grib_set(gaid,
'latitudeOfFirstGridPointInDegrees',lafirst)
1254 CALL grib_set(gaid,
'latitudeOfLastGridPointInDegrees',lalast)
1258 sdx = this%grid%grid%dx*ratio
1259 sdy = this%grid%grid%dy*ratio
1261 ex = abs(((sdx - nint(sdx))/ratio*this%dim%nx)/ &
1262 (this%grid%grid%xmax-this%grid%grid%xmin))
1263 ey = abs(((sdy - nint(sdy))/ratio*this%dim%ny)/ &
1264 (this%grid%grid%ymax-this%grid%grid%ymin))
1266 IF (ex > tol .OR. ey > tol)
THEN 1269 "griddim_export_gribapi, error on x "//&
1270 trim(to_char(ex))//
"/"//t2c(tol))
1272 "griddim_export_gribapi, error on y "//&
1273 trim(to_char(ey))//
"/"//t2c(tol))
1287 "griddim_export_gribapi, increments not given: inaccurate!")
1288 CALL grib_set_missing(gaid,
'Di')
1289 CALL grib_set_missing(gaid,
'Dj')
1290 CALL grib_set(gaid,
'ijDirectionIncrementGiven',0)
1294 "griddim_export_gribapi, setting increments: "// &
1295 trim(to_char(this%grid%grid%dx))//
' '//trim(to_char(this%grid%grid%dy)))
1297 CALL grib_set(gaid,
'ijDirectionIncrementGiven',1)
1298 CALL grib_set(gaid,
'iDirectionIncrement',nint(this%grid%grid%dx*ratio))
1299 CALL grib_set(gaid,
'jDirectionIncrement',nint(this%grid%grid%dy*ratio))
1306 CASE (
'polar_stereographic',
'lambert',
'albers')
1308 CALL grib_set(gaid,
'DxInMetres', this%grid%grid%dx)
1309 CALL grib_set(gaid,
'DyInMetres', this%grid%grid%dy)
1310 CALL grib_set(gaid,
'ijDirectionIncrementGiven',1)
1312 CALL grib_set_dmiss(gaid,
'Latin1InDegrees',this%grid%proj%polar%latin1)
1313 CALL grib_set_dmiss(gaid,
'Latin2InDegrees',this%grid%proj%polar%latin2)
1315 CALL grib_set(gaid,
'projectionCenterFlag',&
1316 this%grid%proj%polar%projection_center_flag, ierr)
1317 IF (ierr /= grib_success)
THEN 1318 CALL grib_set(gaid,
'projectionCentreFlag',&
1319 this%grid%proj%polar%projection_center_flag)
1324 CALL grib_set(gaid,
'LoVInDegrees',this%grid%proj%lov)
1326 IF (editionnumber == 2)
THEN 1327 CALL grib_set(gaid,
'LaDInDegrees',this%grid%proj%polar%lad)
1331 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
1334 IF (editionnumber == 1)
THEN 1335 CALL long_reset_m180_360(lofirst)
1336 ELSE IF (editionnumber == 2)
THEN 1337 CALL long_reset_0_360(lofirst)
1339 CALL grib_set(gaid,
'longitudeOfFirstGridPointInDegrees',lofirst)
1340 CALL grib_set(gaid,
'latitudeOfFirstGridPointInDegrees',lafirst)
1346 CALL grib_set(gaid,
'DxInMetres', this%grid%grid%dx)
1347 CALL grib_set(gaid,
'DyInMetres', this%grid%grid%dy)
1348 CALL grib_set(gaid,
'ijDirectionIncrementGiven',1)
1353 CALL grib_set(gaid,
'LaDInDegrees',this%grid%proj%polar%lad)
1356 CALL get_unproj_first(this, iscansnegatively, jscanspositively, &
1358 CALL grib_set(gaid,
'longitudeOfFirstGridPointInDegrees',lofirst)
1359 CALL grib_set(gaid,
'latitudeOfFirstGridPointInDegrees',lafirst)
1360 CALL get_unproj_last(this, iscansnegatively, jscanspositively, &
1362 CALL grib_set(gaid,
'longitudeOfLastGridPointInDegrees',lolast)
1363 CALL grib_set(gaid,
'latitudeOfLastGridPointInDegrees',lalast)
1365 IF (editionnumber == 2)
THEN 1366 CALL grib_set(gaid,
'orientationOfTheGridInDegrees',0.)
1371 CALL grib_set(gaid,
'datum',0)
1372 CALL get_val(this, zone=zone, lov=reflon)
1373 CALL grib_set(gaid,
'referenceLongitude',nint(reflon/1.0d6))
1374 CALL grib_set(gaid,
'falseEasting',this%grid%proj%xoff)
1375 CALL grib_set(gaid,
'falseNorthing',this%grid%proj%yoff)
1377 CALL grib_set(gaid,
'iDirectionIncrement',this%grid%grid%dx)
1378 CALL grib_set(gaid,
'jDirectionIncrement',this%grid%grid%dy)
1382 CALL grib_set(gaid,
'zone',zone)
1384 IF (iscansnegatively == 0)
THEN 1385 lofirst = this%grid%grid%xmin
1386 lolast = this%grid%grid%xmax
1388 lofirst = this%grid%grid%xmax
1389 lolast = this%grid%grid%xmin
1391 IF (jscanspositively == 0)
THEN 1392 lafirst = this%grid%grid%ymax
1393 lalast = this%grid%grid%ymin
1395 lafirst = this%grid%grid%ymin
1396 lalast = this%grid%grid%ymax
1399 CALL grib_set(gaid,
'eastingOfFirstGridPoint',lofirst)
1400 CALL grib_set(gaid,
'eastingOfLastGridPoint',lolast)
1401 CALL grib_set(gaid,
'northingOfFirstGridPoint',lafirst)
1402 CALL grib_set(gaid,
'northingOfLastGridPoint',lalast)
1406 "griddim_export_gribapi, grid type "//trim(this%grid%proj%proj_type)//
" not supported")
1413 IF (editionnumber == 1)
THEN 1415 CALL grib_get(gaid,
"NV",nv)
1417 CALL l4f_category_log(this%category,l4f_debug,
"griddim_export_gribapi, coding "// &
1418 trim(to_char(nv))//
" vertical coordinate parameters")
1424 SELECT CASE (this%grid%proj%proj_type)
1427 CASE (
'polar_stereographic')
1429 CASE (
'rotated_ll',
'stretched_ll',
'lambert',
'albers')
1431 CASE (
'stretched_rotated_ll')
1438 CALL grib_set(gaid,
"pvlLocation",pvl)
1440 CALL l4f_category_log(this%category,l4f_debug,
"griddim_export_gribapi, coding "// &
1441 trim(to_char(pvl))//
" as vertical coordinate parameter location")
1449 SUBROUTINE grib_set_dmiss(gaid, key, val, default, factor)
1450 INTEGER,
INTENT(in) :: gaid
1451 CHARACTER(len=*),
INTENT(in) :: key
1452 DOUBLE PRECISION,
INTENT(in) :: val
1453 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: default
1454 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: factor
1459 IF (
PRESENT(factor))
THEN 1460 CALL grib_set(gaid, key, val*factor, ierr)
1462 CALL grib_set(gaid, key, val, ierr)
1464 ELSE IF (
PRESENT(default))
THEN 1465 CALL grib_set(gaid, key, default, ierr)
1468 END SUBROUTINE grib_set_dmiss
1470 SUBROUTINE grib_set_imiss(gaid, key, value, default)
1471 INTEGER,
INTENT(in) :: gaid
1472 CHARACTER(len=*),
INTENT(in) :: key
1473 INTEGER,
INTENT(in) :: value
1474 INTEGER,
INTENT(in),
OPTIONAL :: default
1478 IF (
c_e(
value))
THEN 1479 CALL grib_set(gaid, key,
value, ierr)
1480 ELSE IF (
PRESENT(default))
THEN 1481 CALL grib_set(gaid, key, default, ierr)
1484 END SUBROUTINE grib_set_imiss
1486 SUBROUTINE griddim_export_ellipsoid(this, gaid)
1487 TYPE(griddim_def),
INTENT(in) :: this
1488 INTEGER,
INTENT(in) :: gaid
1490 INTEGER :: ellips_type, ierr
1491 DOUBLE PRECISION :: r1, r2, f
1493 CALL get_val(this, ellips_smaj_axis=r1, ellips_flatt=f, ellips_type=ellips_type)
1495 IF (editionnumber == 2)
THEN 1498 CALL grib_set_missing(gaid,
'scaleFactorOfRadiusOfSphericalEarth', ierr)
1499 CALL grib_set_missing(gaid,
'scaledValueOfRadiusOfSphericalEarth', ierr)
1500 CALL grib_set_missing(gaid,
'scaleFactorOfEarthMajorAxis', ierr)
1501 CALL grib_set_missing(gaid,
'scaledValueOfEarthMajorAxis', ierr)
1502 CALL grib_set_missing(gaid,
'scaleFactorOfEarthMinorAxis', ierr)
1503 CALL grib_set_missing(gaid,
'scaledValueOfEarthMinorAxis', ierr)
1505 SELECT CASE(ellips_type)
1507 CALL grib_set(gaid,
'shapeOfTheEarth', 4)
1509 CALL grib_set(gaid,
'shapeOfTheEarth', 5)
1511 IF (f == 0.0d0)
THEN 1512 IF (r1 == 6367470.0d0)
THEN 1513 CALL grib_set(gaid,
'shapeOfTheEarth', 0)
1514 ELSE IF (r1 == 6371229.0d0)
THEN 1515 CALL grib_set(gaid,
'shapeOfTheEarth', 6)
1517 CALL grib_set(gaid,
'shapeOfTheEarth', 1)
1518 CALL grib_set(gaid,
'scaleFactorOfRadiusOfSphericalEarth', 2)
1519 CALL grib_set(gaid,
'scaledValueOfRadiusOfSphericalEarth', int(r1*100.0d0))
1522 IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0)
THEN 1523 CALL grib_set(gaid,
'shapeOfTheEarth', 2)
1525 CALL grib_set(gaid,
'shapeOfTheEarth', 3)
1527 CALL grib_set(gaid,
'scaleFactorOfEarthMajorAxis', 5)
1528 CALL grib_set(gaid,
'scaledValueOfEarthMajorAxis', &
1530 CALL grib_set(gaid,
'scaleFactorOfEarthMinorAxis', 5)
1531 CALL grib_set(gaid,
'scaledValueOfEarthMinorAxis', &
1539 IF (r1 == 6367470.0d0 .AND. f == 0.0d0)
THEN 1540 CALL grib_set(gaid,
'earthIsOblate', 0)
1541 ELSE IF (r1 == 6378160.0d0 .AND. f == 1.0d0/297.0d0)
THEN 1542 CALL grib_set(gaid,
'earthIsOblate', 1)
1543 ELSE IF (f == 0.0d0)
THEN 1544 CALL grib_set(gaid,
'earthIsOblate', 0)
1545 CALL l4f_category_log(this%category,l4f_warn,
'desired spherical Earth radius & 1546 ¬ supported in grib 1, coding with default radius of 6367470 m')
1548 CALL grib_set(gaid,
'earthIsOblate', 1)
1550 ¬ supported in grib 1, coding with default iau65 ellipsoid')
1555 END SUBROUTINE griddim_export_ellipsoid
1557 SUBROUTINE get_unproj_first(this, iScansNegatively, jScansPositively, &
1559 TYPE(griddim_def),
INTENT(in) :: this
1560 INTEGER,
INTENT(in) :: iScansNegatively, jScansPositively
1561 DOUBLE PRECISION,
INTENT(out) :: loFirst, laFirst
1564 IF (iscansnegatively == 0)
THEN 1565 IF (jscanspositively == 0)
THEN 1566 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymax, lofirst, lafirst)
1568 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymin, lofirst, lafirst)
1571 IF (jscanspositively == 0)
THEN 1572 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymax, lofirst, lafirst)
1574 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymin, lofirst, lafirst)
1580 END SUBROUTINE get_unproj_first
1582 SUBROUTINE get_unproj_last(this, iScansNegatively, jScansPositively, &
1584 TYPE(griddim_def),
INTENT(in) :: this
1585 INTEGER,
INTENT(in) :: iScansNegatively, jScansPositively
1586 DOUBLE PRECISION,
INTENT(out) :: loLast, laLast
1589 IF (iscansnegatively == 0)
THEN 1590 IF (jscanspositively == 0)
THEN 1591 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymin, lolast, lalast)
1593 CALL unproj(this, this%grid%grid%xmax, this%grid%grid%ymax, lolast, lalast)
1596 IF (jscanspositively == 0)
THEN 1597 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymin, lolast, lalast)
1599 CALL unproj(this, this%grid%grid%xmin, this%grid%grid%ymax, lolast, lalast)
1605 END SUBROUTINE get_unproj_last
1607 END SUBROUTINE griddim_export_gribapi
1613 SUBROUTINE griddim_import_gdal(this, gdalid, gdal_options)
1615 TYPE(griddim_def),
INTENT(inout) :: this
1616 TYPE(gdalrasterbandh),
INTENT(in) :: gdalid
1617 TYPE(gdal_file_id_options),
INTENT(in) :: gdal_options
1619 TYPE(gdaldataseth) :: hds
1620 REAL(kind=c_double) :: geotrans(6)
1621 INTEGER(kind=c_int) :: offsetx, offsety
1624 hds = gdalgetbanddataset(gdalid)
1625 ier = gdalgetgeotransform(hds, geotrans)
1629 'griddim_import_gdal: error in accessing gdal dataset')
1633 IF (geotrans(3) /= 0.0_c_double .OR. geotrans(5) /= 0.0_c_double)
THEN 1635 'griddim_import_gdal: dataset has a non-diagonal transformation matrix, unsupported')
1640 CALL gdaldatasetbbsize_f(hds, gdal_options%xmin, gdal_options%ymin, &
1641 gdal_options%xmax, gdal_options%ymax, &
1642 this%dim%nx, this%dim%ny, offsetx, offsety, &
1643 this%grid%grid%xmin, this%grid%grid%ymin, this%grid%grid%xmax, this%grid%grid%ymax)
1645 IF (this%dim%nx == 0 .OR. this%dim%ny == 0)
THEN 1647 'griddim_import_gdal: requested bounding box '//t2c(gdal_options%xmin)//
','// &
1648 t2c(gdal_options%ymin)//
','//t2c(gdal_options%xmax)//
','//&
1649 t2c(gdal_options%ymax))
1651 'determines an empty dataset '//t2c(this%dim%nx)//
'x'//t2c(this%dim%ny))
1678 this%grid%proj%proj_type =
'regular_ll' 1680 CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
1681 this%grid%grid%component_flag = 0
1683 END SUBROUTINE griddim_import_gdal
1688 SUBROUTINE griddim_display(this)
1689 TYPE(griddim_def),
INTENT(in) :: this
1691 print*,
"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>" 1697 print*,
"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>" 1699 END SUBROUTINE griddim_display
1702 #define VOL7D_POLY_TYPE TYPE(grid_def) 1703 #define VOL7D_POLY_TYPES _grid 1704 #include "array_utilities_inc.F90" 1705 #undef VOL7D_POLY_TYPE 1706 #undef VOL7D_POLY_TYPES 1708 #define VOL7D_POLY_TYPE TYPE(griddim_def) 1709 #define VOL7D_POLY_TYPES _griddim 1710 #include "array_utilities_inc.F90" 1711 #undef VOL7D_POLY_TYPE 1712 #undef VOL7D_POLY_TYPES 1726 SUBROUTINE griddim_wind_unrot(this, rot_mat)
1728 TYPE(griddim_def),
INTENT(in) :: this
1729 DOUBLE PRECISION,
POINTER :: rot_mat(:,:,:)
1731 DOUBLE PRECISION :: dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j
1732 DOUBLE PRECISION :: lat_factor
1733 INTEGER :: i, j, ip1, im1, jp1, jm1;
1735 IF (this%dim%nx <= 0 .OR. this%dim%ny <= 0 .OR. &
1736 .NOT.
ASSOCIATED(this%dim%lon) .OR. .NOT.
ASSOCIATED(this%dim%lat))
THEN 1737 CALL l4f_category_log(this%category, l4f_error,
'rotate_uv must be called after correct init of griddim object')
1742 ALLOCATE(rot_mat(this%dim%nx, this%dim%ny, 4))
1744 DO j = 1, this%dim%ny
1745 jp1 = min(j+1, this%dim%ny)
1747 DO i = 1, this%dim%nx
1748 ip1 = min(i+1, this%dim%nx)
1751 dlat_i = this%dim%lat(ip1,j) - this%dim%lat(im1,j)
1752 dlat_j = this%dim%lat(i,jp1) - this%dim%lat(i,jm1)
1754 dlon_i = this%dim%lon(ip1,j) - this%dim%lon(im1,j)
1757 dlon_j = this%dim%lon(i,jp1) - this%dim%lon(i,jm1)
1762 lat_factor = cos(degrad*this%dim%lat(i,j))
1763 dlon_i = dlon_i * lat_factor
1764 dlon_j = dlon_j * lat_factor
1766 dist_i = sqrt(dlon_i*dlon_i + dlat_i*dlat_i)
1767 dist_j = sqrt(dlon_j*dlon_j + dlat_j*dlat_j)
1769 IF (dist_i > 0.d0)
THEN 1770 rot_mat(i,j,1) = dlon_i/dist_i
1771 rot_mat(i,j,2) = dlat_i/dist_i
1773 rot_mat(i,j,1) = 0.0d0
1774 rot_mat(i,j,2) = 0.0d0
1776 IF (dist_j > 0.d0)
THEN 1777 rot_mat(i,j,3) = dlon_j/dist_j
1778 rot_mat(i,j,4) = dlat_j/dist_j
1780 rot_mat(i,j,3) = 0.0d0
1781 rot_mat(i,j,4) = 0.0d0
1787 END SUBROUTINE griddim_wind_unrot
1791 SUBROUTINE griddim_zoom_coord(this, ilon, ilat, flon, flat, ix, iy, fx, fy)
1792 TYPE(griddim_def),
INTENT(in) :: this
1793 DOUBLE PRECISION,
INTENT(in) :: ilon, ilat, flon, flat
1794 INTEGER,
INTENT(out) :: ix, iy, fx, fy
1796 DOUBLE PRECISION :: ix1, iy1, fx1, fy1
1799 CALL proj(this, ilon, ilat, ix1, iy1)
1800 CALL proj(this, flon, flat, fx1, fy1)
1802 CALL griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
1804 END SUBROUTINE griddim_zoom_coord
1808 SUBROUTINE griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
1809 TYPE(griddim_def),
INTENT(in) :: this
1810 DOUBLE PRECISION,
INTENT(in) :: ix1, iy1, fx1, fy1
1811 INTEGER,
INTENT(out) :: ix, iy, fx, fy
1813 INTEGER :: lix, liy, lfx, lfy
1816 lix = nint((ix1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
1817 liy = nint((iy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
1818 lfx = nint((fx1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
1819 lfy = nint((fy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
1826 END SUBROUTINE griddim_zoom_projcoord
1832 SUBROUTINE long_reset_0_360(lon)
1833 DOUBLE PRECISION,
INTENT(inout) :: lon
1835 IF (.NOT.
c_e(lon))
RETURN 1836 DO WHILE(lon < 0.0d0)
1839 DO WHILE(lon >= 360.0d0)
1843 END SUBROUTINE long_reset_0_360
1849 SUBROUTINE long_reset_m180_360(lon)
1850 DOUBLE PRECISION,
INTENT(inout) :: lon
1852 IF (.NOT.
c_e(lon))
RETURN 1853 DO WHILE(lon < -180.0d0)
1856 DO WHILE(lon >= 360.0d0)
1860 END SUBROUTINE long_reset_m180_360
1883 SUBROUTINE long_reset_m180_180(lon)
1884 DOUBLE PRECISION,
INTENT(inout) :: lon
1886 IF (.NOT.
c_e(lon))
RETURN 1887 DO WHILE(lon < -180.0d0)
1890 DO WHILE(lon >= 180.0d0)
1894 END SUBROUTINE long_reset_m180_180
1897 SUBROUTINE long_reset_to_cart_closest(lon, lonref)
1898 DOUBLE PRECISION,
INTENT(inout) :: lon
1899 DOUBLE PRECISION,
INTENT(in) :: lonref
1901 IF (.NOT.
c_e(lon) .OR. .NOT.
c_e(lonref))
RETURN 1902 IF (abs(lon-lonref) < 180.0d0)
RETURN 1903 lon = lon - nint((lon-lonref)/360.0d0)*360.0d0
1905 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.
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.