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)
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)
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)
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
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)
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)
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)
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)
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)
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)
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: "//&
1191 "griddim_export_gribapi, dlat relative error: "//&
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)
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)
1458 TYPE(gdalrasterbandh
),
INTENT(in) :: gdalid
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)
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)
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)
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)
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.
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.