61 character (len=255),
parameter:: subcategory=
"grid_class"
71 type(geo_proj
) ::
proj
72 type(grid_rect
) :: grid
73 integer :: category = 0
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, &
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)
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)
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)
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)
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)
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)
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)
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)
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)
501 INTEGER,
INTENT(in) :: unit
508 END SUBROUTINE griddim_read_unit
515 SUBROUTINE griddim_write_unit(this, unit)
517 INTEGER,
INTENT(in) :: unit
524 END SUBROUTINE griddim_write_unit
530 FUNCTION griddim_central_lon(this) RESULT(lon)
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)
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)
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)
688 res = .NOT.(this == that)
690 END FUNCTION griddim_ne
698 SUBROUTINE griddim_import_grid_id(this, ingrid_id)
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)
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)
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)
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)
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 "//&
1272 "griddim_export_gribapi, error on y "//&
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)
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, &
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, &
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)
1616 TYPE(gdalrasterbandh
),
INTENT(in) :: gdalid
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)
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)
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)
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)
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.
Definitions of constants and functions for working with missing values.
Functions that return a trimmed CHARACTER representation of the input variable.
Constructors of the corresponding objects.
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...
Operatore di valore assoluto di un intervallo.
Print a brief description on stdout.
This module defines an abstract interface to different drivers for access to files containing gridded...
Derived type containing driver-specific options for gdal.
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.
Module for defining the extension and coordinates of a rectangular georeferenced grid.
This object completely describes a grid on a geographic projection.
Method for returning the contents of the object.
Restituiscono il valore dell'oggetto in forma di stringa stampabile.
classe per la gestione del logging
Copy an object, creating a fully new instance.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
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...
Costanti fisiche (DOUBLEPRECISION).
Destructors of the corresponding objects.
Import griddim object from grid_id.