27 DOUBLE PRECISION :: latin1, latin2
28 DOUBLE PRECISION :: lad
29 DOUBLE PRECISION :: lon1, lat1
30 INTEGER :: projection_center_flag
31END TYPE geo_proj_polar
34 DOUBLE PRECISION :: longitude_south_pole, latitude_south_pole, angle_rotation
35END TYPE geo_proj_rotated
37TYPE geo_proj_stretched
38 DOUBLE PRECISION :: latitude_stretch_pole, longitude_stretch_pole, stretch_factor
39END TYPE geo_proj_stretched
43 DOUBLE PRECISION :: rf, a
44 DOUBLE PRECISION :: f, e2, e1, ep2, e11, e12, e13, e14, e4, e6, ef0, ef1, ef2, ef3, k0
45END TYPE geo_proj_ellips
48 CHARACTER(len=80) :: proj_type=cmiss
49 DOUBLE PRECISION :: xoff, yoff
50 DOUBLE PRECISION :: lov
51 TYPE(geo_proj_rotated) :: rotated
52 TYPE(geo_proj_stretched) :: stretched
53 TYPE(geo_proj_polar) :: polar
55 TYPE(geo_proj_ellips) :: ellips
61 MODULE PROCEDURE geo_proj_delete
66 MODULE PROCEDURE geo_proj_copy
71 MODULE PROCEDURE geo_proj_get_val
76 MODULE PROCEDURE geo_proj_set_val
81 MODULE PROCEDURE geo_proj_c_e
86 MODULE PROCEDURE geo_proj_write_unit
91 MODULE PROCEDURE geo_proj_read_unit
96 MODULE PROCEDURE geo_proj_display
102 MODULE PROCEDURE geo_proj_proj
108 MODULE PROCEDURE geo_proj_unproj
114INTERFACE OPERATOR (==)
115 MODULE PROCEDURE geo_proj_eq
121INTERFACE OPERATOR (/=)
122 MODULE PROCEDURE geo_proj_ne
126INTEGER,
PARAMETER,
PUBLIC :: &
127 geo_proj_unit_degree = 0, & !< coordinate unit is degrees (longitude periodic over 360 degrees)
128 geo_proj_unit_meter = 1
131INTEGER,
PARAMETER :: nellips = 41
136INTEGER,
PARAMETER,
PUBLIC :: &
137ellips_merit = 1, & !< constants for predefined ellipsoids MERIT 1983
144ellips_mod_airy = 8, &
146ellips_aust_sa = 10, &
149ellips_bess_nam = 13, &
154ellips_engelis = 18, &
155ellips_evrst30 = 19, &
156ellips_evrst48 = 20, &
157ellips_evrst56 = 21, &
158ellips_evrst69 = 22, &
159ellips_evrstss = 23, &
160ellips_fschr60 = 24, &
161ellips_fschr60m = 25, &
162ellips_fschr68 = 26, &
163ellips_helmert = 27, &
170ellips_new_intl = 34, &
171ellips_plessis = 35, &
173ellips_walbeck = 37, &
179DOUBLE PRECISION,
PARAMETER,
PRIVATE :: &
222DOUBLE PRECISION,
PARAMETER,
PRIVATE :: &
266DOUBLE PRECISION,
PARAMETER,
PRIVATE :: k0=0.9996d0
269PUBLIC geo_proj, geo_proj_rotated, geo_proj_stretched, geo_proj_polar, &
283FUNCTION geo_proj_new(proj_type, lov, zone, xoff, yoff, &
284 longitude_south_pole, latitude_south_pole, angle_rotation, &
285 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
286 latin1, latin2, lad, projection_center_flag, &
287 ellips_smaj_axis, ellips_flatt, ellips_type)
RESULT(this)
288CHARACTER(len=*),
INTENT(in),
OPTIONAL :: proj_type
289DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lov
290INTEGER,
INTENT(in),
OPTIONAL :: zone
291DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xoff
292DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: yoff
293DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_south_pole
294DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_south_pole
295DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: angle_rotation
296DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_stretch_pole
297DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_stretch_pole
298DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: stretch_factor
299DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin1
300DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin2
301DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lad
302INTEGER,
INTENT(in),
OPTIONAL :: projection_center_flag
303DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_smaj_axis
304DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_flatt
305INTEGER,
INTENT(in),
OPTIONAL :: ellips_type
307TYPE(geo_proj) :: this
309this%proj_type = optio_c(proj_type, len(this%proj_type))
312this%lov = optio_d(lov)
313CALL set_val(this, zone=zone, lov=lov)
316this%xoff = optio_d(xoff)
317IF (.NOT.
c_e(this%xoff)) this%xoff = 0.0d0
318this%yoff = optio_d(yoff)
319IF (.NOT.
c_e(this%yoff)) this%yoff = 0.0d0
321this%rotated%longitude_south_pole = optio_d(longitude_south_pole)
322this%rotated%latitude_south_pole = optio_d(latitude_south_pole)
323this%rotated%angle_rotation = optio_d(angle_rotation)
324this%stretched%longitude_stretch_pole = optio_d(longitude_stretch_pole)
325this%stretched%latitude_stretch_pole = optio_d(latitude_stretch_pole)
326this%stretched%stretch_factor = optio_d(stretch_factor)
327this%polar%latin1 = optio_d(latin1)
328this%polar%latin2 = optio_d(latin2)
329this%polar%lad = optio_d(lad)
330this%polar%projection_center_flag = optio_l(projection_center_flag)
333CALL ellips_compute(this%ellips)
334CALL set_val(this, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
335 ellips_type=ellips_type)
337END FUNCTION geo_proj_new
342SUBROUTINE ellips_compute(this, a, f)
343TYPE(geo_proj_ellips),
INTENT(inout) :: this
344DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: a, f
346IF (
PRESENT(a) .AND.
PRESENT(f))
THEN
349ELSE IF (
PRESENT(a))
THEN
357this%e2 = 2.0d0*this%f - this%f*this%f
358this%e1 = this%f/(2.0d0 - this%f)
359this%ep2 = this%e2/(1.0d0 - this%e2)
360this%e11 = 3.0d0*this%e1/2.0d0 - 27.0d0*this%e1*this%e1*this%e1/32.0d0
361this%e12 = 21.0d0*this%e1*this%e1/16.0d0 &
362 - 55.0d0*this%e1*this%e1*this%e1*this%e1/32.0d0
363this%e13 = 151.0d0*this%e1*this%e1*this%e1/96.0d0
364this%e14 = 1097.0d0*this%e1*this%e1*this%e1*this%e1/512.0d0
365this%e4 = this%e2*this%e2
366this%e6 = this%e2*this%e4
367this%ef0 = this%a*(1.0d0 - 0.25d0*this%e2&
368 *(1.0d0+this%e2/16.0d0*(3.0d0 + 1.25d0*this%e2)))
369this%ef1 = this%a*(0.375d0*this%e2 &
370 *(1.0d0 + 0.25d0*this%e2*(1.0d0 + 0.46875d0*this%e2)))
371this%ef2 = this%a*(0.05859375d0*this%e2*this%e2*(1.0d0 + 0.75d0*this%e2))
372this%ef3 = this%a*this%e2*this%e2*this%e2*35.0d0/3072.0d0
374END SUBROUTINE ellips_compute
378SUBROUTINE geo_proj_delete(this)
379TYPE(geo_proj),
INTENT(inout) :: this
381this%proj_type = cmiss
385this%rotated%longitude_south_pole = dmiss
386this%rotated%latitude_south_pole = dmiss
387this%rotated%angle_rotation = dmiss
388this%stretched%longitude_stretch_pole = dmiss
389this%stretched%latitude_stretch_pole = dmiss
390this%stretched%stretch_factor = dmiss
391this%polar%latin1 = dmiss
392this%polar%latin2 = dmiss
393this%polar%lad = dmiss
394this%polar%projection_center_flag = imiss
396END SUBROUTINE geo_proj_delete
400SUBROUTINE geo_proj_copy(this, that)
401TYPE(geo_proj),
INTENT(in) :: this
402TYPE(geo_proj),
INTENT(out) :: that
406END SUBROUTINE geo_proj_copy
409SUBROUTINE geo_proj_get_val(this, &
410 proj_type, lov, zone, xoff, yoff, &
411 longitude_south_pole, latitude_south_pole, angle_rotation, &
412 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
413 latin1, latin2, lad, projection_center_flag, &
414 ellips_smaj_axis, ellips_flatt, ellips_type, unit)
415TYPE(geo_proj),
INTENT(in) :: this
416CHARACTER(len=*),
OPTIONAL :: proj_type
417DOUBLE PRECISION,
OPTIONAL :: lov
418INTEGER,
OPTIONAL :: zone
419DOUBLE PRECISION,
OPTIONAL :: xoff
420DOUBLE PRECISION,
OPTIONAL :: yoff
421DOUBLE PRECISION,
OPTIONAL :: longitude_south_pole
422DOUBLE PRECISION,
OPTIONAL :: latitude_south_pole
423DOUBLE PRECISION,
OPTIONAL :: angle_rotation
424DOUBLE PRECISION,
OPTIONAL :: longitude_stretch_pole
425DOUBLE PRECISION,
OPTIONAL :: latitude_stretch_pole
426DOUBLE PRECISION,
OPTIONAL :: stretch_factor
427DOUBLE PRECISION,
OPTIONAL :: latin1
428DOUBLE PRECISION,
OPTIONAL :: latin2
429DOUBLE PRECISION,
OPTIONAL :: lad
430INTEGER,
OPTIONAL :: projection_center_flag
431DOUBLE PRECISION,
OPTIONAL :: ellips_smaj_axis
432DOUBLE PRECISION,
OPTIONAL :: ellips_flatt
433INTEGER,
OPTIONAL :: ellips_type
434INTEGER,
OPTIONAL :: unit
438IF (
PRESENT(proj_type)) proj_type = this%proj_type
440IF (
PRESENT(lov) .AND.
PRESENT(zone))
THEN
441 zone = nint((this%lov + 183.0d0)/6.0d0)
442 lov = this%lov - zone*6.0d0 - 183.0d0
443 IF (abs(lov) < 1.0d-6) lov = 0.0d-6
444ELSE IF (
PRESENT(lov))
THEN
446ELSE IF (
PRESENT(zone))
THEN
447 zone = nint((this%lov + 183.0d0)/6.0d0)
450IF (
PRESENT(xoff)) xoff = this%xoff
451IF (
PRESENT(yoff)) yoff = this%yoff
452IF (
PRESENT(longitude_south_pole)) longitude_south_pole = this%rotated%longitude_south_pole
453IF (
PRESENT(latitude_south_pole)) latitude_south_pole = this%rotated%latitude_south_pole
454IF (
PRESENT(angle_rotation)) angle_rotation = this%rotated%angle_rotation
455IF (
PRESENT(longitude_stretch_pole)) longitude_stretch_pole = this%stretched%longitude_stretch_pole
456IF (
PRESENT(latitude_stretch_pole)) latitude_stretch_pole = this%stretched%latitude_stretch_pole
457IF (
PRESENT(stretch_factor)) stretch_factor = this%stretched%stretch_factor
458IF (
PRESENT(latin1)) latin1 = this%polar%latin1
459IF (
PRESENT(latin2)) latin2 = this%polar%latin2
460IF (
PRESENT(lad)) lad = this%polar%lad
461IF (
PRESENT(projection_center_flag)) projection_center_flag = this%polar%projection_center_flag
463IF (
PRESENT(ellips_smaj_axis)) ellips_smaj_axis = this%ellips%a
464IF (
PRESENT(ellips_flatt)) ellips_flatt = this%ellips%f
465IF (
PRESENT(ellips_type))
THEN
468 IF (this%ellips%f == 1.0d0/rf(i) .AND. this%ellips%a == a(i))
THEN
475IF (
PRESENT(unit))
THEN
476 SELECT CASE(this%proj_type)
477 CASE(
"regular_ll",
"rotated_ll")
478 unit = geo_proj_unit_degree
479 CASE(
"lambert",
"polar_stereographic",
"UTM")
480 unit = geo_proj_unit_meter
487END SUBROUTINE geo_proj_get_val
490SUBROUTINE geo_proj_set_val(this, &
491 proj_type, lov, zone, xoff, yoff, &
492 longitude_south_pole, latitude_south_pole, angle_rotation, &
493 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
494 latin1, latin2, lad, projection_center_flag, &
495 ellips_smaj_axis, ellips_flatt, ellips_type)
496TYPE(geo_proj),
INTENT(inout) :: this
497CHARACTER(len=*),
OPTIONAL :: proj_type
498DOUBLE PRECISION,
OPTIONAL :: lov
499INTEGER,
INTENT(in),
OPTIONAL :: zone
500DOUBLE PRECISION,
OPTIONAL :: xoff
501DOUBLE PRECISION,
OPTIONAL :: yoff
502DOUBLE PRECISION,
OPTIONAL :: longitude_south_pole
503DOUBLE PRECISION,
OPTIONAL :: latitude_south_pole
504DOUBLE PRECISION,
OPTIONAL :: angle_rotation
505DOUBLE PRECISION,
OPTIONAL :: longitude_stretch_pole
506DOUBLE PRECISION,
OPTIONAL :: latitude_stretch_pole
507DOUBLE PRECISION,
OPTIONAL :: stretch_factor
508DOUBLE PRECISION,
OPTIONAL :: latin1
509DOUBLE PRECISION,
OPTIONAL :: latin2
510DOUBLE PRECISION,
OPTIONAL :: lad
511INTEGER,
OPTIONAL :: projection_center_flag
512DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_smaj_axis
513DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_flatt
514INTEGER,
INTENT(in),
OPTIONAL :: ellips_type
517DOUBLE PRECISION :: llov
523IF (
c_e(llov) .AND.
c_e(lzone))
THEN
524 this%lov = llov + zone*6.0d0 - 183.0d0
525ELSE IF (
c_e(llov))
THEN
527ELSE IF (
c_e(lzone))
THEN
528 this%lov = lzone*6.0d0 - 183.0d0
532IF (
PRESENT(ellips_smaj_axis))
THEN
534 CALL ellips_compute(this%ellips, ellips_smaj_axis, ellips_flatt)
535ELSE IF (
PRESENT(ellips_type))
THEN
536 IF (ellips_type > 0 .AND. ellips_type <= nellips)
THEN
538 CALL ellips_compute(this%ellips, a(ellips_type), 1.0d0/rf(ellips_type))
540 CALL ellips_compute(this%ellips)
546END SUBROUTINE geo_proj_set_val
549FUNCTION geo_proj_c_e(this)
RESULT(c_e)
550TYPE(geo_proj),
INTENT(in) :: this
553c_e = this%proj_type /= cmiss
555END FUNCTION geo_proj_c_e
562SUBROUTINE geo_proj_read_unit(this, unit)
563TYPE(geo_proj),
INTENT(out) :: this
564INTEGER,
INTENT(in) :: unit
566CHARACTER(len=40) :: form
568INQUIRE(unit, form=form)
569IF (form ==
'FORMATTED')
THEN
575END SUBROUTINE geo_proj_read_unit
582SUBROUTINE geo_proj_write_unit(this, unit)
583TYPE(geo_proj),
INTENT(in) :: this
584INTEGER,
INTENT(in) :: unit
586CHARACTER(len=40) :: form
588INQUIRE(unit, form=form)
589IF (form ==
'FORMATTED')
THEN
595END SUBROUTINE geo_proj_write_unit
598SUBROUTINE geo_proj_display(this)
599TYPE(geo_proj),
INTENT(in) :: this
601print*,
"<<<<<<<<<<<<<<< ",
t2c(this%proj_type,
"undefined projection"), &
604IF (
c_e(this%xoff) .AND. this%xoff /= 0.0d0)
THEN
605 print*,
"False easting",this%xoff
607IF (
c_e(this%yoff) .AND. this%yoff /= 0.0d0)
THEN
608 print*,
"False northing",this%yoff
611IF (
c_e(this%lov))
THEN
612 print*,
"Central Meridian",this%lov
615IF (this%proj_type ==
'rotated_ll' .OR. this%proj_type ==
'stretched_rotated_ll')
THEN
616 print*,
"Rotated projection:"
617 print*,
"lon of south pole",this%rotated%longitude_south_pole
618 print*,
"lat of south pole",this%rotated%latitude_south_pole
619 print*,
"angle of rotation",this%rotated%angle_rotation
622IF (this%proj_type ==
'stretched_ll' .OR. this%proj_type ==
'stretched_rotated_ll')
THEN
623 print*,
"Stretched projection:"
624 print*,
"lon of stretch pole",this%stretched%longitude_stretch_pole
625 print*,
"lat of stretch pole",this%stretched%latitude_stretch_pole
626 print*,
"stretching factor",this%stretched%stretch_factor
629IF (this%proj_type ==
'lambert' .OR. this%proj_type ==
'polar_stereographic')
THEN
630 print*,
"Polar projection:"
631 IF (
c_e(this%polar%latin1) .OR.
c_e(this%polar%latin2))
THEN
632 print*,
"lat of intersections",this%polar%latin1,this%polar%latin2
634 IF (
c_e(this%polar%lad))
THEN
635 print*,
"isometric latitude",this%polar%lad
637 IF (iand(this%polar%projection_center_flag, 128) == 0)
THEN
644IF (this%proj_type ==
'mercator')
THEN
645 IF (
c_e(this%polar%lad))
THEN
646 print*,
"isometric latitude",this%polar%lad
650IF (this%ellips%f == 0.0d0)
THEN
651 print*,
"Spherical Earth:"
652 print*,
"Radius (m)",this%ellips%a
655 print*,
"Flattening",this%ellips%f
656 print*,
"Reverse of flattening",1.0d0/this%ellips%f
657 print*,
"Semi-major axis (m)",this%ellips%a
661END SUBROUTINE geo_proj_display
666ELEMENTAL SUBROUTINE geo_proj_proj(this, lon, lat, x, y)
667TYPE(geo_proj),
INTENT(in) :: this
669DOUBLE PRECISION,
INTENT(in) :: lon, lat
671DOUBLE PRECISION,
INTENT(out) :: x, y
673SELECT CASE(this%proj_type)
676 CALL proj_regular_ll(lon, lat, x, y)
679 CALL proj_rotated_ll(lon, lat, x, y, this%rotated%longitude_south_pole, &
680 this%rotated%latitude_south_pole, this%rotated%angle_rotation)
683 CALL proj_lambert(lon, lat, x, y, this%polar%latin1, &
684 this%polar%latin2, this%lov, this%polar%lad, &
685 this%polar%projection_center_flag)
687CASE(
"polar_stereographic")
688 CALL proj_polar_stereographic(lon, lat, x, y, this%lov, &
689 this%polar%lad, this%polar%projection_center_flag)
692 CALL proj_mercator(lon, lat, x, y, this%lov, this%polar%lad)
695 CALL proj_utm(lon, lat, x, y, this%lov, this%xoff, this%yoff, this%ellips)
703END SUBROUTINE geo_proj_proj
708ELEMENTAL SUBROUTINE geo_proj_unproj(this, x, y, lon, lat)
709TYPE(geo_proj),
INTENT(in) :: this
711DOUBLE PRECISION,
INTENT(in) :: x, y
713DOUBLE PRECISION,
INTENT(out) :: lon, lat
715SELECT CASE(this%proj_type)
718 CALL unproj_regular_ll(x, y, lon, lat)
721 CALL unproj_rotated_ll(x, y, lon, lat, this%rotated%longitude_south_pole, &
722 this%rotated%latitude_south_pole, this%rotated%angle_rotation)
725 CALL unproj_lambert(x, y, lon, lat, this%polar%latin1, &
726 this%polar%latin2, this%lov, this%polar%lad, &
727 this%polar%projection_center_flag)
729CASE(
"polar_stereographic")
730 CALL unproj_polar_stereographic(x, y, lon, lat, this%lov, &
731 this%polar%lad, this%polar%projection_center_flag)
734 CALL unproj_mercator(x, y, lon, lat, this%lov, this%polar%lad)
737 CALL unproj_utm(x, y, lon, lat, this%lov, this%xoff, this%yoff, this%ellips)
745END SUBROUTINE geo_proj_unproj
748ELEMENTAL FUNCTION geo_proj_eq(this, that)
RESULT(eq)
749TYPE(geo_proj),
INTENT(in) :: this, that
752eq = this%proj_type == that%proj_type .AND. this%xoff == that%xoff .AND. &
753 this%yoff == that%yoff .AND. geo_lon_eq(this%lov, that%lov) .AND. &
754 geo_lon_eq(this%rotated%longitude_south_pole, that%rotated%longitude_south_pole) .AND. &
755 this%rotated%latitude_south_pole == that%rotated%latitude_south_pole .AND. &
756 this%rotated%angle_rotation == that%rotated%angle_rotation .AND. &
757 this%stretched%latitude_stretch_pole == that%stretched%latitude_stretch_pole .AND. &
758 geo_lon_eq(this%stretched%longitude_stretch_pole, that%stretched%longitude_stretch_pole) .AND. &
759 this%stretched%stretch_factor == that%stretched%stretch_factor .AND. &
760 this%polar%latin1 == that%polar%latin1 .AND. &
761 this%polar%latin2 == that%polar%latin2 .AND. &
762 this%polar%lad == that%polar%lad .AND. &
763 this%polar%projection_center_flag == that%polar%projection_center_flag .AND. &
764 this%ellips%f == that%ellips%f .AND. this%ellips%a == that%ellips%a
766END FUNCTION geo_proj_eq
769ELEMENTAL FUNCTION geo_proj_ne(this, that)
RESULT(ne)
770TYPE(geo_proj),
INTENT(in) :: this, that
773ne = .NOT. (this == that)
775END FUNCTION geo_proj_ne
780ELEMENTAL FUNCTION geo_lon_eq(l1, l2)
RESULT(eq)
781DOUBLE PRECISION,
INTENT(in) :: l1
782DOUBLE PRECISION,
INTENT(in) :: l2
787eq = modulo(l2-l1, 360.0d0) == 0.0d0
789END FUNCTION geo_lon_eq
795ELEMENTAL SUBROUTINE proj_regular_ll(lon,lat,x,y)
796DOUBLE PRECISION,
INTENT(in) :: lon,lat
797DOUBLE PRECISION,
INTENT(out) :: x,y
802END SUBROUTINE proj_regular_ll
804ELEMENTAL SUBROUTINE unproj_regular_ll(x,y,lon,lat)
805DOUBLE PRECISION,
INTENT(in) :: x,y
806DOUBLE PRECISION,
INTENT(out) :: lon,lat
811END SUBROUTINE unproj_regular_ll
814ELEMENTAL SUBROUTINE proj_rotated_ll(lon,lat,x,y, &
815 longitude_south_pole, latitude_south_pole, angle_rotation)
816DOUBLE PRECISION,
INTENT(in) :: lon,lat
817DOUBLE PRECISION,
INTENT(out) :: x,y
818DOUBLE PRECISION,
INTENT(in) :: longitude_south_pole, latitude_south_pole, &
821DOUBLE PRECISION :: cy0,sy0,rx,srx,crx,sy,cy,l_south_pole
825l_south_pole = (latitude_south_pole+90.)*degrad
827rx = degrad*(lon - longitude_south_pole)
831sy0 = sin(l_south_pole)
832cy0 = cos(l_south_pole)
837x = raddeg*atan2(cy*srx, cy0*cy*crx+sy0*sy)
838y = raddeg*asin(cy0*sy - sy0*cy*crx)
839x = x + angle_rotation
841END SUBROUTINE proj_rotated_ll
843ELEMENTAL SUBROUTINE unproj_rotated_ll(x,y,lon,lat,&
844 longitude_south_pole, latitude_south_pole, angle_rotation)
845DOUBLE PRECISION,
INTENT(in) :: x,y
846DOUBLE PRECISION,
INTENT(out) :: lon,lat
847DOUBLE PRECISION,
INTENT(in) :: longitude_south_pole, latitude_south_pole, &
850DOUBLE PRECISION :: cy0, sy0, l_south_pole, xr
852xr = (x - angle_rotation)*degrad
855l_south_pole = (latitude_south_pole+90.)*degrad
857cy0 = cos(l_south_pole)
858sy0 = sin(l_south_pole)
860lat = raddeg*asin(sy0*cos(degrad*y)*cos(xr)+cy0*sin(degrad*y))
861lon = longitude_south_pole + &
862 raddeg*asin(sin(xr)*cos(degrad*y)/cos(degrad*lat))
864END SUBROUTINE unproj_rotated_ll
867ELEMENTAL SUBROUTINE proj_stretched_ll(lon,lat,x,y, &
868 longitude_stretch_pole, latitude_stretch_pole, stretch_factor)
869DOUBLE PRECISION,
INTENT(in) :: lon,lat
870DOUBLE PRECISION,
INTENT(out) :: x,y
871DOUBLE PRECISION,
INTENT(in) :: longitude_stretch_pole, latitude_stretch_pole, &
874DOUBLE PRECISION :: csq
876csq = stretch_factor**2
878y = raddeg*asin((1.0d0 - csq + (1.0d0 + csq)*sin(degrad*lat)) / &
879 (1.0d0 + csq + (1.0d0 - csq)*sin(degrad*lat)))
881END SUBROUTINE proj_stretched_ll
883ELEMENTAL SUBROUTINE unproj_stretched_ll(x,y,lon,lat,&
884 longitude_stretch_pole, latitude_stretch_pole, stretch_factor)
885DOUBLE PRECISION,
INTENT(in) :: x,y
886DOUBLE PRECISION,
INTENT(out) :: lon,lat
887DOUBLE PRECISION,
INTENT(in) :: longitude_stretch_pole, latitude_stretch_pole, &
890DOUBLE PRECISION :: csq
892csq = stretch_factor**2
895lat = raddeg*asin((csq - 1.0d0 + (csq + 1.0d0)*sin(degrad*y)) / &
896 (csq + 1.0d0 + (csq - 1.0d0)*sin(degrad*y)))
898END SUBROUTINE unproj_stretched_ll
909ELEMENTAL SUBROUTINE proj_lambert(lon,lat,x,y, &
910 latin1, latin2, lov, lad, projection_center_flag)
911DOUBLE PRECISION,
INTENT(in) :: lon,lat
912DOUBLE PRECISION,
INTENT(out) :: x,y
913DOUBLE PRECISION,
INTENT(in) :: latin1, latin2, lov, lad
914INTEGER,
INTENT(in) :: projection_center_flag
916DOUBLE PRECISION :: n, f, ro0, ro, cs1, cs2, cs3, pollat, angle, cot
917DOUBLE PRECISION,
PARAMETER :: epsy = 1.0d-100
919IF (iand(projection_center_flag, 128) == 0)
THEN
920 pollat = 90.d0*degrad
922 pollat = -90.d0*degrad
924cs1 = cos(degrad*latin1)
925cs2 = tan(pi*.25d0 + degrad*latin1*.5d0)
927IF (latin1 == latin2)
THEN
928 n = sin(degrad*latin1)
930 n = log(cs1/cos(degrad*latin2)) / &
931 log(tan(pi*.25d0 + degrad*latin2*.5d0) / cs2)
933f = cs1*cs2**n/n*rearth
934angle = pi*.25d0 + pollat*.5d0
935cot = cos(angle)/sin(angle)
942angle = pi*.25d0 + degrad*lat*.5d0
943cot = cos(angle)/sin(angle)
950cs3 = degrad*n*(lon - lov)
955END SUBROUTINE proj_lambert
957ELEMENTAL SUBROUTINE unproj_lambert(x,y,lon,lat, &
958 latin1, latin2, lov, lad, projection_center_flag)
959DOUBLE PRECISION,
INTENT(in) :: x,y
960DOUBLE PRECISION,
INTENT(out) :: lon,lat
961DOUBLE PRECISION,
INTENT(in) :: latin1, latin2, lov, lad
962INTEGER,
INTENT(in) :: projection_center_flag
964DOUBLE PRECISION :: n, f, ro0, ro, theta, cs1, cs2, pollat, angle, cot
965DOUBLE PRECISION,
PARAMETER :: epsy = 1.0d-100
969IF (iand(projection_center_flag, 128) == 0)
THEN
970 pollat = 90.d0*degrad
972 pollat = -90.d0*degrad
974cs1 = cos(degrad*latin1)
975cs2 = tan(pi*.25d0 + degrad*latin1*.5d0)
977IF (latin1 == latin2)
THEN
978 n = sin(degrad*latin1)
980 n = log(cs1/cos(degrad*latin2)) / &
981 log(tan(pi*.25d0 + degrad*latin2*.5d0) / cs2)
983f = cs1*cs2**n/n*rearth
984angle = pi*.25d0 + pollat*.5d0
985cot = cos(angle)/sin(angle)
992ro = sign(sqrt(x*x + (ro0-y)*(ro0-y)), n)
993theta = raddeg*atan2(x, ro0-y)
996lat = raddeg*(2.d0*atan((f/ro)**(1.d0/n)) - pi*.5d0)
998END SUBROUTINE unproj_lambert
1002ELEMENTAL SUBROUTINE proj_polar_stereographic(lon,lat,x,y, &
1003 lov, lad, projection_center_flag)
1004DOUBLE PRECISION,
INTENT(in) :: lon,lat
1005DOUBLE PRECISION,
INTENT(out) :: x,y
1006DOUBLE PRECISION,
INTENT(in) :: lov, lad
1007INTEGER,
INTENT(in) :: projection_center_flag
1009DOUBLE PRECISION :: k, pollat
1011IF (iand(projection_center_flag, 128) == 0)
THEN
1012 pollat = 90.d0*degrad
1014 pollat = -90.d0*degrad
1017k = 2.0d0*rearth/(1.0d0 + sin(pollat)*sin(degrad*lat) + &
1018 cos(pollat)*cos(degrad*lat)*cos(degrad*(lon - lov)))
1019x = k*cos(degrad*lat)*sin(degrad*(lon - lov))
1020y = k*(cos(pollat)*sin(degrad*lat) - &
1021 sin(pollat)*cos(degrad*lat)*cos(degrad*(lon - lov)))
1023END SUBROUTINE proj_polar_stereographic
1025ELEMENTAL SUBROUTINE unproj_polar_stereographic(x,y,lon,lat, &
1026 lov, lad, projection_center_flag)
1027DOUBLE PRECISION,
INTENT(in) :: x,y
1028DOUBLE PRECISION,
INTENT(out) :: lon,lat
1029DOUBLE PRECISION,
INTENT(in) :: lov, lad
1030INTEGER,
INTENT(in) :: projection_center_flag
1032DOUBLE PRECISION :: ro, c, pollat
1034IF (iand(projection_center_flag, 128) == 0)
THEN
1035 pollat = 90.d0*degrad
1037 pollat = -90.d0*degrad
1040ro = sqrt(x**2 + y**2)
1041c = 2.0d0*atan(ro/(2.0d0*rearth))
1042lat = raddeg*asin(cos(c)*sin(pollat)+y*sin(c)*cos(pollat)/ro)
1043lon = lov + raddeg*atan2(x*sin(c), &
1044 (ro*cos(pollat)*cos(c)-y*sin(pollat)*sin(c)))
1046END SUBROUTINE unproj_polar_stereographic
1050ELEMENTAL SUBROUTINE proj_mercator(lon, lat, x, y, lov, lad)
1051DOUBLE PRECISION,
INTENT(in) :: lon,lat
1052DOUBLE PRECISION,
INTENT(out) :: x,y
1053DOUBLE PRECISION,
INTENT(in) :: lov
1054DOUBLE PRECISION,
INTENT(in) :: lad
1056DOUBLE PRECISION :: scx, scy
1058scy = cos(degrad*lad)*rearth
1061y = log(tan(0.25d0*pi + 0.5d0*lat*degrad))*scy
1063END SUBROUTINE proj_mercator
1065ELEMENTAL SUBROUTINE unproj_mercator(x, y, lon, lat, lov, lad)
1066DOUBLE PRECISION,
INTENT(in) :: x,y
1067DOUBLE PRECISION,
INTENT(out) :: lon,lat
1068DOUBLE PRECISION,
INTENT(in) :: lov
1069DOUBLE PRECISION,
INTENT(in) :: lad
1071DOUBLE PRECISION :: scx, scy
1073scy = cos(degrad*lad)*rearth
1076lat = 2.0d0*atan(exp(y/scy))*raddeg-90.0d0
1078END SUBROUTINE unproj_mercator
1081ELEMENTAL SUBROUTINE proj_utm(lon, lat, x, y, lov, false_e, false_n, ellips)
1082DOUBLE PRECISION,
INTENT(in) :: lon,lat
1083DOUBLE PRECISION,
INTENT(out) :: x,y
1084DOUBLE PRECISION,
INTENT(in) :: lov
1085DOUBLE PRECISION,
INTENT(in) :: false_e, false_n
1086TYPE(geo_proj_ellips),
INTENT(in) :: ellips
1088DOUBLE PRECISION :: deltalon, p
1089DOUBLE PRECISION :: n, t, t2, c, m, a1, a2, a3, a4, a5, a6, sinp, cosp, tanp
1092deltalon = degrad*(lon - lov)
1100n = ellips%a/sqrt(1.0d0 - ellips%e2*sinp*sinp)
1102c = ellips%ep2*cosp*cosp
1107m = ellips%ef0*p - ellips%ef1*sin(2.0d0*p) + ellips%ef2*sin(4.0d0*p) - &
1108 ellips%ef3*sin(6.0d0*p)
1118x = k0*n*(a1 + (1.0d0 - t + c)*a3/6.0d0 &
1119 + (5.0d0 - 18.0d0*t + t2 + 72.0d0*c - 58.0d0*ellips%ep2) &
1120 *a5/120.0d0) + false_e
1122 *(a2/2.0d0 + (5.0d0 - t + 9.0d0*c + &
1123 4.0d0*c*c)*a4/24.0d0 + (61.0d0 - 58.0d0*t + t2 + &
1124 600.0d0*c - 330.0d0*ellips%ep2)*a6/720.0d0))
1127END SUBROUTINE proj_utm
1129ELEMENTAL SUBROUTINE unproj_utm(x, y, lon, lat, lov, false_e, false_n, ellips)
1130DOUBLE PRECISION,
INTENT(IN) :: x, y
1131DOUBLE PRECISION,
INTENT(OUT) :: lon, lat
1132DOUBLE PRECISION,
INTENT(in) :: lov
1133DOUBLE PRECISION,
INTENT(in) :: false_e, false_n
1134TYPE(geo_proj_ellips),
INTENT(in) :: ellips
1136DOUBLE PRECISION :: xm, ym, m, u, p1, c1, c2, t1, t2, n1, &
1137 sinp1, cosp1, tanp1, sin2p1, r0, r1, d, d2, d3, d4, d5, d6, p, l
1148u = m/(ellips%a*(1.0d0-ellips%e2/4.0d0 - &
1149 3.0d0*ellips%e4/64.0d0 - 5.0d0*ellips%e6/256.0d0))
1150p1 = u + ellips%e11*sin(2.0d0*u) + ellips%e12*sin(4.0d0*u) + &
1151 ellips%e13*sin(6.0d0*u) + ellips%e14*sin(8.0d0*u)
1155c1 = ellips%ep2*cosp1**2
1160r0 = 1.0d0-ellips%e2*sin2p1
1161n1 = ellips%a/sqrt(r0)
1163r1 = r0/(1.0d0-ellips%e2)
1177p = p1 - (r1*tanp1) * (d2/2.0d0 &
1178 - (5.0d0 + 3.0d0*t1 + 10.0d0*c1 - 4.0d0*c2 &
1179 - 9.0d0*ellips%ep2)*d4/24.0d0 &
1180 + (61.0d0 + 90.0d0*t1 + 298.0d0*c1 + 45.0d0*t2 &
1181 - 252d0*ellips%ep2 - 3.0d0*c2)*d6/720.0d0)
1183l = (d - (1.0d0 + 2.0d0*t1 + c1)*d3/6.0d0 &
1184 + (5.0d0 - 2.0d0*c1 + 28.0d0*t1 - 3.0d0*c2 &
1185 + 8.0d0*ellips%ep2 + 24.0d0*t2)*d5/120.0d0)/cosp1
1188END SUBROUTINE unproj_utm
1190END MODULE geo_proj_class
Set of functions that return a trimmed CHARACTER representation of the input variable.
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.
Method for returning the contents of the object.
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.
Utilities for CHARACTER variables.
Costanti fisiche (DOUBLEPRECISION).
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.