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