61character (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
92INTERFACE OPERATOR (==)
93 MODULE PROCEDURE grid_eq, griddim_eq
99INTERFACE 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
184PUBLIC proj,
unproj, griddim_unproj, griddim_gen_coord, &
185 griddim_zoom_coord, griddim_zoom_projcoord, &
189PUBLIC OPERATOR(==),
OPERATOR(/=)
190PUBLIC count_distinct, pack_distinct, count_and_pack_distinct, &
191 map_distinct, map_inv_distinct,
index
193PUBLIC griddim_central_lon, griddim_set_central_lon
197SUBROUTINE 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, &
205TYPE(griddim_def),
INTENT(inout) :: this
206INTEGER,
INTENT(in),
OPTIONAL :: nx
207INTEGER,
INTENT(in),
OPTIONAL :: ny
208DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xmin
209DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xmax
210DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ymin
211DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ymax
212DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: dx
213DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: dy
216INTEGER,
INTENT(in),
OPTIONAL :: component_flag
217CHARACTER(len=*),
INTENT(in),
OPTIONAL :: proj_type
218DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lov
219INTEGER,
INTENT(in),
OPTIONAL :: zone
220DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xoff
221DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: yoff
222DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_south_pole
223DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_south_pole
224DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: angle_rotation
225DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_stretch_pole
226DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_stretch_pole
227DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: stretch_factor
228DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin1
229DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin2
230DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lad
231INTEGER,
INTENT(in),
OPTIONAL :: projection_center_flag
232DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_smaj_axis
233DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_flatt
234INTEGER,
INTENT(in),
OPTIONAL :: ellips_type
235CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
237CHARACTER(len=512) :: a_name
239IF (
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))
245this%category=l4f_category_get(a_name)
248this%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)
257this%grid%grid = grid_rect_new( &
258 xmin, xmax, ymin, ymax, dx, dy, component_flag)
260this%dim = grid_dim_new(nx, ny)
266END SUBROUTINE griddim_init
270SUBROUTINE griddim_delete(this)
271TYPE(griddim_def),
INTENT(inout) :: this
274CALL delete(this%grid%proj)
275CALL delete(this%grid%grid)
277call l4f_category_delete(this%category)
279END SUBROUTINE griddim_delete
283SUBROUTINE griddim_copy(this, that, categoryappend)
286CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
288CHARACTER(len=512) :: a_name
292CALL copy(this%grid%proj, that%grid%proj)
293CALL copy(this%grid%grid, that%grid%grid)
294CALL copy(this%dim, that%dim)
297IF (
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))
303that%category=l4f_category_get(a_name)
305END SUBROUTINE griddim_copy
310ELEMENTAL SUBROUTINE griddim_coord_proj(this, lon, lat, x, y)
311TYPE(griddim_def),
INTENT(in) :: this
313DOUBLE PRECISION,
INTENT(in) :: lon, lat
315DOUBLE PRECISION,
INTENT(out) :: x, y
317CALL proj(this%grid%proj, lon, lat, x, y)
319END SUBROUTINE griddim_coord_proj
324ELEMENTAL SUBROUTINE griddim_coord_unproj(this, x, y, lon, lat)
325TYPE(griddim_def),
INTENT(in) :: this
327DOUBLE PRECISION,
INTENT(in) :: x, y
329DOUBLE PRECISION,
INTENT(out) :: lon, lat
331CALL unproj(this%grid%proj, x, y, lon, lat)
333END SUBROUTINE griddim_coord_unproj
358SUBROUTINE griddim_unproj(this)
359TYPE(griddim_def),
INTENT(inout) :: this
361IF (.NOT.
c_e(this%dim%nx) .OR. .NOT.
c_e(this%dim%ny))
RETURN
363CALL griddim_unproj_internal(this)
365END SUBROUTINE griddim_unproj
368SUBROUTINE griddim_unproj_internal(this)
369TYPE(griddim_def),
INTENT(inout) ::this
371DOUBLE PRECISION :: x(this%dim%nx,this%dim%ny), y(this%dim%nx,this%dim%ny)
373CALL grid_rect_coordinates(this%grid%grid, x, y)
374CALL unproj(this, x, y, this%dim%lon, this%dim%lat)
376END SUBROUTINE griddim_unproj_internal
380SUBROUTINE 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)
387TYPE(griddim_def),
INTENT(in) :: this
388INTEGER,
INTENT(out),
OPTIONAL :: nx
389INTEGER,
INTENT(out),
OPTIONAL :: ny
391DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: xmin, xmax, ymin, ymax
393DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: dx, dy
396INTEGER,
INTENT(out),
OPTIONAL :: component_flag
397TYPE(geo_proj),
INTENT(out),
OPTIONAL :: proj
398CHARACTER(len=*),
INTENT(out),
OPTIONAL :: proj_type
399DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: lov
400INTEGER,
INTENT(out),
OPTIONAL :: zone
401DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: xoff
402DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: yoff
403DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: longitude_south_pole
404DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: latitude_south_pole
405DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: angle_rotation
406DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: longitude_stretch_pole
407DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: latitude_stretch_pole
408DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: stretch_factor
409DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: latin1
410DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: latin2
411DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: lad
412INTEGER,
INTENT(out),
OPTIONAL :: projection_center_flag
413DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: ellips_smaj_axis
414DOUBLE PRECISION,
INTENT(out),
OPTIONAL :: ellips_flatt
415INTEGER,
INTENT(out),
OPTIONAL :: ellips_type
417IF (
PRESENT(nx)) nx = this%dim%nx
418IF (
PRESENT(ny)) ny = this%dim%ny
420IF (
PRESENT(
proj))
proj = this%grid%proj
422CALL 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)
436END SUBROUTINE griddim_get_val
440SUBROUTINE 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)
447TYPE(griddim_def),
INTENT(inout) :: this
448INTEGER,
INTENT(in),
OPTIONAL :: nx
449INTEGER,
INTENT(in),
OPTIONAL :: ny
451DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xmin, xmax, ymin, ymax
453DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: dx, dy
456INTEGER,
INTENT(in),
OPTIONAL :: component_flag
457CHARACTER(len=*),
INTENT(in),
OPTIONAL :: proj_type
458DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lov
459INTEGER,
INTENT(in),
OPTIONAL :: zone
460DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xoff
461DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: yoff
462DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_south_pole
463DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_south_pole
464DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: angle_rotation
465DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_stretch_pole
466DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_stretch_pole
467DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: stretch_factor
468DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin1
469DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin2
470DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lad
471INTEGER,
INTENT(in),
OPTIONAL :: projection_center_flag
472DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_smaj_axis
473DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_flatt
474INTEGER,
INTENT(in),
OPTIONAL :: ellips_type
476IF (
PRESENT(nx)) this%dim%nx = nx
477IF (
PRESENT(ny)) this%dim%ny = ny
479CALL 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)
492END SUBROUTINE griddim_set_val
499SUBROUTINE griddim_read_unit(this, unit)
501INTEGER,
INTENT(in) :: unit
508END SUBROUTINE griddim_read_unit
515SUBROUTINE griddim_write_unit(this, unit)
517INTEGER,
INTENT(in) :: unit
524END SUBROUTINE griddim_write_unit
530FUNCTION griddim_central_lon(this)
RESULT(lon)
533DOUBLE PRECISION :: lon
535CALL griddim_pistola_central_lon(this, lon)
537END FUNCTION griddim_central_lon
543SUBROUTINE griddim_set_central_lon(this, lonref)
545DOUBLE PRECISION,
INTENT(in) :: lonref
547DOUBLE PRECISION :: lon
549CALL griddim_pistola_central_lon(this, lon, lonref)
551END SUBROUTINE griddim_set_central_lon
555SUBROUTINE griddim_pistola_central_lon(this, lon, lonref)
557DOUBLE PRECISION,
INTENT(inout) :: lon
558DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lonref
561DOUBLE PRECISION :: lonsp, latsp, londelta, lov, lonrot
562CHARACTER(len=80) :: ptype
565CALL get_val(this%grid%proj, unit=unit)
566IF (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)
573ELSE 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
614END SUBROUTINE griddim_pistola_central_lon
620SUBROUTINE griddim_gen_coord(this, x, y)
622DOUBLE PRECISION,
INTENT(out) :: x(:,:)
623DOUBLE PRECISION,
INTENT(out) :: y(:,:)
626CALL grid_rect_coordinates(this%grid%grid, x, y)
628END SUBROUTINE griddim_gen_coord
632SUBROUTINE griddim_steps(this, nx, ny, dx, dy)
634INTEGER,
INTENT(in) :: nx
635INTEGER,
INTENT(in) :: ny
636DOUBLE PRECISION,
INTENT(out) :: dx
637DOUBLE PRECISION,
INTENT(out) :: dy
639CALL grid_rect_steps(this%grid%grid, nx, ny, dx, dy)
641END SUBROUTINE griddim_steps
645SUBROUTINE griddim_setsteps(this)
648CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
650END SUBROUTINE griddim_setsteps
655ELEMENTAL FUNCTION grid_eq(this, that)
RESULT(res)
656TYPE(
grid_def),
INTENT(IN) :: this, that
659res = this%proj == that%proj .AND. &
660 this%grid == that%grid
665ELEMENTAL FUNCTION griddim_eq(this, that)
RESULT(res)
669res = this%grid == that%grid .AND. &
672END FUNCTION griddim_eq
675ELEMENTAL FUNCTION grid_ne(this, that)
RESULT(res)
676TYPE(
grid_def),
INTENT(IN) :: this, that
679res = .NOT.(this == that)
684ELEMENTAL FUNCTION griddim_ne(this, that)
RESULT(res)
688res = .NOT.(this == that)
690END FUNCTION griddim_ne
698SUBROUTINE griddim_import_grid_id(this, ingrid_id)
703TYPE(
grid_id),
INTENT(in) :: ingrid_id
705#ifdef HAVE_LIBGRIBAPI
709TYPE(gdalrasterbandh) :: gdalid
713#ifdef HAVE_LIBGRIBAPI
714gaid = grid_id_get_gaid(ingrid_id)
715IF (
c_e(gaid))
CALL griddim_import_gribapi(this, gaid)
718gdalid = grid_id_get_gdalid(ingrid_id)
719IF (gdalassociated(gdalid))
CALL griddim_import_gdal(this, gdalid, &
720 grid_id_get_gdal_options(ingrid_id))
723END SUBROUTINE griddim_import_grid_id
730SUBROUTINE griddim_export_grid_id(this, outgrid_id)
735TYPE(
grid_id),
INTENT(inout) :: outgrid_id
737#ifdef HAVE_LIBGRIBAPI
741TYPE(gdalrasterbandh) :: gdalid
744#ifdef HAVE_LIBGRIBAPI
745gaid = grid_id_get_gaid(outgrid_id)
746IF (
c_e(gaid))
CALL griddim_export_gribapi(this, gaid)
749gdalid = grid_id_get_gdalid(outgrid_id)
754END SUBROUTINE griddim_export_grid_id
757#ifdef HAVE_LIBGRIBAPI
759SUBROUTINE griddim_import_gribapi(this, gaid)
762INTEGER,
INTENT(in) :: gaid
764DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, x1, y1, orient
765INTEGER :: EditionNumber, iScansNegatively, jScansPositively, zone, datum, &
769CALL grib_get(gaid,
'typeOfGrid', this%grid%proj%proj_type)
772 "griddim_import_gribapi, grid type "//trim(this%grid%proj%proj_type))
774CALL grib_get(gaid,
'GRIBEditionNumber',editionnumber)
776IF (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)
785CALL grib_get(gaid,
'Ni', this%dim%nx)
786CALL grib_get(gaid,
'Nj', this%dim%ny)
787CALL griddim_import_ellipsoid(this, gaid)
789CALL grib_get(gaid,
'iScansNegatively',iscansnegatively)
790CALL grib_get(gaid,
'jScansPositively',jscanspositively)
791CALL grib_get(gaid,
'uvRelativeToGrid',this%grid%grid%component_flag)
794CALL grib_get_dmiss(gaid,
'longitudeOfSouthernPoleInDegrees', &
795 this%grid%proj%rotated%longitude_south_pole)
796CALL grib_get_dmiss(gaid,
'latitudeOfSouthernPoleInDegrees', &
797 this%grid%proj%rotated%latitude_south_pole)
798CALL grib_get_dmiss(gaid,
'angleOfRotationInDegrees', &
799 this%grid%proj%rotated%angle_rotation)
804IF (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)
811ELSE 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
823SELECT CASE (this%grid%proj%proj_type)
826CASE (
'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)
863CASE (
'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")
1047SUBROUTINE grib_get_dmiss(gaid, key, value)
1048INTEGER,
INTENT(in) :: gaid
1049CHARACTER(len=*),
INTENT(in) :: key
1050DOUBLE PRECISION,
INTENT(out) :: value
1054CALL grib_get(gaid, key,
value, ierr)
1055IF (ierr /= grib_success)
value = dmiss
1057END SUBROUTINE grib_get_dmiss
1059SUBROUTINE grib_get_imiss(gaid, key, value)
1060INTEGER,
INTENT(in) :: gaid
1061CHARACTER(len=*),
INTENT(in) :: key
1062INTEGER,
INTENT(out) :: value
1066CALL grib_get(gaid, key,
value, ierr)
1067IF (ierr /= grib_success)
value = imiss
1069END SUBROUTINE grib_get_imiss
1072SUBROUTINE griddim_import_ellipsoid(this, gaid)
1074INTEGER,
INTENT(in) :: gaid
1076INTEGER :: shapeofearth, iv, is
1077DOUBLE PRECISION :: r1, r2
1079IF (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)
1134END SUBROUTINE griddim_import_ellipsoid
1137END SUBROUTINE griddim_import_gribapi
1141SUBROUTINE griddim_export_gribapi(this, gaid)
1144INTEGER,
INTENT(inout) :: gaid
1146INTEGER :: EditionNumber, iScansNegatively, jScansPositively, nv, pvl, zone, ierr
1147DOUBLE PRECISION :: loFirst, loLast, laFirst, laLast, reflon
1148DOUBLE PRECISION :: sdx, sdy, ratio, tol, ex, ey
1152CALL grib_get(gaid,
'GRIBEditionNumber', editionnumber)
1154IF (editionnumber == 2)
CALL grib_set(gaid,
'shapeOfTheEarth', 6)
1155CALL grib_set(gaid,
'typeOfGrid', this%grid%proj%proj_type)
1158 "griddim_export_gribapi, grid type "//this%grid%proj%proj_type)
1161IF (this%grid%proj%proj_type ==
'unstructured_grid')
THEN
1162 CALL grib_set(gaid,
'numberOfDataPoints', this%dim%nx)
1170IF (editionnumber == 1)
THEN
1172ELSE IF (editionnumber == 2)
THEN
1179CALL griddim_export_ellipsoid(this, gaid)
1180CALL grib_set(gaid,
'Ni', this%dim%nx)
1181CALL grib_set(gaid,
'Nj', this%dim%ny)
1182CALL grib_set(gaid,
'uvRelativeToGrid',this%grid%grid%component_flag)
1184CALL grib_get(gaid,
'iScansNegatively',iscansnegatively)
1185CALL grib_get(gaid,
'jScansPositively',jscanspositively)
1190CALL grib_set_dmiss(gaid,
'longitudeOfSouthernPoleInDegrees', &
1191 this%grid%proj%rotated%longitude_south_pole, 0.0d0)
1192CALL grib_set_dmiss(gaid,
'latitudeOfSouthernPoleInDegrees', &
1193 this%grid%proj%rotated%latitude_south_pole, -90.0d0)
1194IF (editionnumber == 1)
THEN
1195 CALL grib_set_dmiss(gaid,
'angleOfRotationInDegrees', &
1196 this%grid%proj%rotated%angle_rotation, 0.0d0)
1197ELSE IF (editionnumber == 2)
THEN
1198 CALL grib_set_dmiss(gaid,
'angleOfRotationOfProjectionInDegrees', &
1199 this%grid%proj%rotated%angle_rotation, 0.0d0)
1205IF (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)
1212ELSE 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)
1222SELECT CASE (this%grid%proj%proj_type)
1225CASE (
'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))
1306CASE (
'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")
1413IF (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")
1449SUBROUTINE grib_set_dmiss(gaid, key, val, default, factor)
1450INTEGER,
INTENT(in) :: gaid
1451CHARACTER(len=*),
INTENT(in) :: key
1452DOUBLE PRECISION,
INTENT(in) :: val
1453DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: default
1454DOUBLE 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)
1464ELSE IF (
PRESENT(default))
THEN
1465 CALL grib_set(gaid, key, default, ierr)
1468END SUBROUTINE grib_set_dmiss
1470SUBROUTINE grib_set_imiss(gaid, key, value, default)
1471INTEGER,
INTENT(in) :: gaid
1472CHARACTER(len=*),
INTENT(in) :: key
1473INTEGER,
INTENT(in) :: value
1474INTEGER,
INTENT(in),
OPTIONAL :: default
1479 CALL grib_set(gaid, key,
value, ierr)
1480ELSE IF (
PRESENT(default))
THEN
1481 CALL grib_set(gaid, key, default, ierr)
1484END SUBROUTINE grib_set_imiss
1486SUBROUTINE griddim_export_ellipsoid(this, gaid)
1488INTEGER,
INTENT(in) :: gaid
1490INTEGER :: ellips_type, ierr
1491DOUBLE PRECISION :: r1, r2, f
1493CALL get_val(this, ellips_smaj_axis=r1, ellips_flatt=f, ellips_type=ellips_type)
1495IF (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')
1555END SUBROUTINE griddim_export_ellipsoid
1557SUBROUTINE get_unproj_first(this, iScansNegatively, jScansPositively, &
1560INTEGER,
INTENT(in) :: iScansNegatively, jScansPositively
1561DOUBLE PRECISION,
INTENT(out) :: loFirst, laFirst
1564IF (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)
1580END SUBROUTINE get_unproj_first
1582SUBROUTINE get_unproj_last(this, iScansNegatively, jScansPositively, &
1585INTEGER,
INTENT(in) :: iScansNegatively, jScansPositively
1586DOUBLE PRECISION,
INTENT(out) :: loLast, laLast
1589IF (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)
1605END SUBROUTINE get_unproj_last
1607END SUBROUTINE griddim_export_gribapi
1613SUBROUTINE griddim_import_gdal(this, gdalid, gdal_options)
1616TYPE(gdalrasterbandh),
INTENT(in) :: gdalid
1619TYPE(gdaldataseth) :: hds
1620REAL(kind=c_double) :: geotrans(6)
1621INTEGER(kind=c_int) :: offsetx, offsety
1624hds = gdalgetbanddataset(gdalid)
1625ier = gdalgetgeotransform(hds, geotrans)
1629 'griddim_import_gdal: error in accessing gdal dataset')
1633IF (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')
1640CALL 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)
1645IF (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))
1678this%grid%proj%proj_type =
'regular_ll'
1680CALL grid_rect_setsteps(this%grid%grid, this%dim%nx, this%dim%ny)
1681this%grid%grid%component_flag = 0
1683END SUBROUTINE griddim_import_gdal
1688SUBROUTINE griddim_display(this)
1691print*,
"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
1697print*,
"<<<<<<<<<<<<<<< griddim >>>>>>>>>>>>>>>>"
1699END 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
1726SUBROUTINE griddim_wind_unrot(this, rot_mat)
1729DOUBLE PRECISION,
POINTER :: rot_mat(:,:,:)
1731DOUBLE PRECISION :: dlat_i, dlat_j,dlon_i,dlon_j,dist_i,dist_j
1732DOUBLE PRECISION :: lat_factor
1733INTEGER :: i, j, ip1, im1, jp1, jm1;
1735IF (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')
1742ALLOCATE(rot_mat(this%dim%nx, this%dim%ny, 4))
1744DO 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
1787END SUBROUTINE griddim_wind_unrot
1791SUBROUTINE griddim_zoom_coord(this, ilon, ilat, flon, flat, ix, iy, fx, fy)
1793DOUBLE PRECISION,
INTENT(in) :: ilon, ilat, flon, flat
1794INTEGER,
INTENT(out) :: ix, iy, fx, fy
1796DOUBLE PRECISION :: ix1, iy1, fx1, fy1
1799CALL proj(this, ilon, ilat, ix1, iy1)
1800CALL proj(this, flon, flat, fx1, fy1)
1802CALL griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
1804END SUBROUTINE griddim_zoom_coord
1808SUBROUTINE griddim_zoom_projcoord(this, ix1, iy1, fx1, fy1, ix, iy, fx, fy)
1810DOUBLE PRECISION,
INTENT(in) :: ix1, iy1, fx1, fy1
1811INTEGER,
INTENT(out) :: ix, iy, fx, fy
1813INTEGER :: lix, liy, lfx, lfy
1816lix = nint((ix1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
1817liy = nint((iy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
1818lfx = nint((fx1-this%grid%grid%xmin)/this%grid%grid%dx) + 1
1819lfy = nint((fy1-this%grid%grid%ymin)/this%grid%grid%dy) + 1
1826END SUBROUTINE griddim_zoom_projcoord
1832SUBROUTINE long_reset_0_360(lon)
1833DOUBLE PRECISION,
INTENT(inout) :: lon
1835IF (.NOT.
c_e(lon))
RETURN
1836DO WHILE(lon < 0.0d0)
1839DO WHILE(lon >= 360.0d0)
1843END SUBROUTINE long_reset_0_360
1849SUBROUTINE long_reset_m180_360(lon)
1850DOUBLE PRECISION,
INTENT(inout) :: lon
1852IF (.NOT.
c_e(lon))
RETURN
1853DO WHILE(lon < -180.0d0)
1856DO WHILE(lon >= 360.0d0)
1860END SUBROUTINE long_reset_m180_360
1883SUBROUTINE long_reset_m180_180(lon)
1884DOUBLE PRECISION,
INTENT(inout) :: lon
1886IF (.NOT.
c_e(lon))
RETURN
1887DO WHILE(lon < -180.0d0)
1890DO WHILE(lon >= 180.0d0)
1894END SUBROUTINE long_reset_m180_180
1897SUBROUTINE long_reset_to_cart_closest(lon, lonref)
1898DOUBLE PRECISION,
INTENT(inout) :: lon
1899DOUBLE PRECISION,
INTENT(in) :: lonref
1901IF (.NOT.
c_e(lon) .OR. .NOT.
c_e(lonref))
RETURN
1902IF (abs(lon-lonref) < 180.0d0)
RETURN
1903lon = lon - nint((lon-lonref)/360.0d0)*360.0d0
1905END SUBROUTINE long_reset_to_cart_closest
Method for testing the existence of the object.
Copy an object, creating a fully new instance.
Destructors of the corresponding objects.
Print a brief description on stdout.
Export griddim object to grid_id.
Method for returning the contents of the object.
Import griddim object from grid_id.
Constructors of the corresponding objects.
Compute forward coordinate transformation from geographical system to projected system.
Read the object from a formatted or unformatted file.
Method for setting the contents of the object.
Compute backward coordinate transformation from projected system to geographical system.
Write the object on a formatted or unformatted file.
Emit log message for a category with specific priority.
Costanti fisiche (DOUBLEPRECISION).
Module for describing geographically referenced regular grids.
Module for defining the extension and coordinates of a rectangular georeferenced grid.
This module defines an abstract interface to different drivers for access to files containing gridded...
classe per la gestione del logging
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
This object, mainly for internal use, describes a grid on a geographical projection,...
This object completely describes a grid on a geographic projection.
Derived type describing the extension of a grid and the geographical coordinates of each point.
Derived type containing driver-specific options for gdal.
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...