libsim  Versione7.1.6
geo_proj_class.F90
1 ! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2 ! authors:
3 ! Davide Cesari <dcesari@arpa.emr.it>
4 ! Paolo Patruno <ppatruno@arpa.emr.it>
5 
6 ! This program is free software; you can redistribute it and/or
7 ! modify it under the terms of the GNU General Public License as
8 ! published by the Free Software Foundation; either version 2 of
9 ! the License, or (at your option) any later version.
10 
11 ! This program is distributed in the hope that it will be useful,
12 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
13 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 ! GNU General Public License for more details.
15 
16 ! You should have received a copy of the GNU General Public License
17 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
18 #include "config.h"
19 MODULE geo_proj_class
24 IMPLICIT NONE
25 
26 TYPE geo_proj_polar
27  DOUBLE PRECISION :: latin1, latin2 ! latitudes at which the projection plane intesects the sphere
28  DOUBLE PRECISION :: lad ! latitudine at which dx and dy (in m) are specified
29  DOUBLE PRECISION :: lon1, lat1 ! stored in order not to forget them, from grib
30  INTEGER :: projection_center_flag ! 0 = northern hemisphere, 128 = southern hemisphere
31 END TYPE geo_proj_polar
32 
33 TYPE geo_proj_rotated
34  DOUBLE PRECISION :: longitude_south_pole, latitude_south_pole, angle_rotation
35 END TYPE geo_proj_rotated
36 
37 TYPE geo_proj_stretched
38  DOUBLE PRECISION :: latitude_stretch_pole, longitude_stretch_pole, stretch_factor
39 END TYPE geo_proj_stretched
40 
41 TYPE geo_proj_ellips
42  PRIVATE
43  DOUBLE PRECISION :: rf, a ! inverse of flattening and semi-major axis
44  DOUBLE PRECISION :: f, e2, e1, ep2, e11, e12, e13, e14, e4, e6, ef0, ef1, ef2, ef3, k0 ! computed parameters
45 END TYPE geo_proj_ellips
46 
47 TYPE geo_proj
48  CHARACTER(len=80) :: proj_type=cmiss ! the projection type
49  DOUBLE PRECISION :: xoff, yoff ! offsets in x and y wrt origin, aka false easting and northing resp.
50  DOUBLE PRECISION :: lov ! line of view (or central meridian, reference longitude, orientation of the grid)
51  TYPE(geo_proj_rotated) :: rotated
52  TYPE(geo_proj_stretched) :: stretched
53  TYPE(geo_proj_polar) :: polar
54 !TYPE (geo_proj_equatorial) :: equatorial ! For projections like Mercator?
55  TYPE(geo_proj_ellips) :: ellips
56 END TYPE geo_proj
57 
58 
60 INTERFACE delete
61  MODULE PROCEDURE geo_proj_delete
62 END INTERFACE
63 
65 INTERFACE copy
66  MODULE PROCEDURE geo_proj_copy
67 END INTERFACE
68 
70 INTERFACE get_val
71  MODULE PROCEDURE geo_proj_get_val
72 END INTERFACE
73 
75 INTERFACE set_val
76  MODULE PROCEDURE geo_proj_set_val
77 END INTERFACE
78 
80 INTERFACE c_e
81  MODULE PROCEDURE geo_proj_c_e
82 END INTERFACE
83 
85 INTERFACE write_unit
86  MODULE PROCEDURE geo_proj_write_unit
87 END INTERFACE
88 
90 INTERFACE read_unit
91  MODULE PROCEDURE geo_proj_read_unit
92 END INTERFACE
93 
95 INTERFACE display
96  MODULE PROCEDURE geo_proj_display
97 END INTERFACE
98 
101 INTERFACE proj
102  MODULE PROCEDURE geo_proj_proj
103 END INTERFACE
104 
107 INTERFACE unproj
108  MODULE PROCEDURE geo_proj_unproj
109 END INTERFACE
110 
114 INTERFACE operator (==)
115  MODULE PROCEDURE geo_proj_eq
116 END INTERFACE
117 
121 INTERFACE operator (/=)
122  MODULE PROCEDURE geo_proj_ne
123 END INTERFACE
124 
125 
126 INTEGER,PARAMETER,PUBLIC :: &
127  geo_proj_unit_degree = 0, & !< coordinate unit is degrees (longitude periodic over 360 degrees)
128  geo_proj_unit_meter = 1
129 
130 
131 INTEGER,PARAMETER :: nellips = 41
132 
133 ! queste costanti vanno usate per specificare l'ellissoide da usare per
134 ! interpretare i dati UTM, gli ellissoidi sono tratti dal pacchetto
135 ! software "proj" (vedi l'uscita a video del comando \c proj \c -le ).
136 INTEGER,PARAMETER,PUBLIC :: &
137 ellips_merit = 1, & !< constants for predefined ellipsoids MERIT 1983
138 ellips_sgs85 = 2, & !< Soviet Geodetic System 85
139 ellips_grs80 = 3, & !< GRS 1980(IUGG, 1980)
140 ellips_iau76 = 4, & !< IAU 1976
141 ellips_airy = 5, & !< Airy 1830
142 ellips_apl4_9 = 6, & !< Appl. Physics. 1965
143 ellips_nwl9d = 7, & !< Naval Weapons Lab., 1965
144 ellips_mod_airy = 8, & !< Modified Airy
145 ellips_andrae = 9, & !< Andrae 1876 (Den., Iclnd.)
146 ellips_aust_sa = 10, & !< Australian Natl & S. Amer. 1969
147 ellips_grs67 = 11, & !< GRS 67(IUGG 1967)
148 ellips_bessel = 12, & !< Bessel 1841
149 ellips_bess_nam = 13, & !< Bessel 1841 (Namibia)
150 ellips_clrk66 = 14, & !< Clarke 1866
151 ellips_clrk80 = 15, & !< Clarke 1880 mod.
152 ellips_cpm = 16, & !< Comm. des Poids et Mesures 1799
153 ellips_delmbr = 17, & !< Delambre 1810 (Belgium)
154 ellips_engelis = 18, & !< Engelis 1985
155 ellips_evrst30 = 19, & !< Everest 1830
156 ellips_evrst48 = 20, & !< Everest 1948
157 ellips_evrst56 = 21, & !< Everest 1956
158 ellips_evrst69 = 22, & !< Everest 1969
159 ellips_evrstss = 23, & !< Everest (Sabah & Sarawak)
160 ellips_fschr60 = 24, & !< Fischer (Mercury Datum) 1960
161 ellips_fschr60m = 25, & !< Modified Fischer 1960
162 ellips_fschr68 = 26, & !< Fischer 1968
163 ellips_helmert = 27, & !< Helmert 1906
164 ellips_hough = 28, & !< Hough
165 ellips_intl = 29, & !< International 1909 (Hayford)
166 ellips_krass = 30, & !< Krassovsky, 1942
167 ellips_kaula = 31, & !< Kaula 1961
168 ellips_lerch = 32, & !< Lerch 1979
169 ellips_mprts = 33, & !< Maupertius 1738
170 ellips_new_intl = 34, & !< New International 1967
171 ellips_plessis = 35, & !< Plessis 1817 (France)
172 ellips_seasia = 36, & !< Southeast Asia
173 ellips_walbeck = 37, & !< Walbeck
174 ellips_wgs60 = 38, & !< WGS 60
175 ellips_wgs66 = 39, & !< WGS 66
176 ellips_wgs72 = 40, & !< WGS 72
177 ellips_wgs84 = 41
178 
179 DOUBLE PRECISION, PARAMETER, PRIVATE :: &
180  rf(nellips)=(/ & ! inverse of flattening for each ellipsoid
181  298.257d0, &
182  298.257d0, &
183  298.257222101d0, &
184  298.257d0, &
185  299.325d0, &
186  298.25d0, &
187  298.25d0, &
188  299.328d0, &
189  300.0d0, &
190  298.25d0, &
191  298.2471674270d0, &
192  299.1528128d0, &
193  299.1528128d0, &
194  294.98d0, &
195  293.4663d0, &
196  334.29d0, &
197  311.5d0, &
198  298.2566d0, &
199  300.8017d0, &
200  300.8017d0, &
201  300.8017d0, &
202  300.8017d0, &
203  300.8017d0, &
204  298.3d0, &
205  298.3d0, &
206  298.3d0, &
207  298.3d0, &
208  297.d0, &
209  297.d0, &
210  298.3d0, &
211  298.24d0, &
212  298.257d0, &
213  191.d0, &
214  298.247d0, &
215  308.641d0, &
216  298.302d0, &
217  302.782d0, &
218  298.3d0, &
219  298.25d0, &
220  298.26d0, &
221  298.257223563d0 /)
222 DOUBLE PRECISION, PARAMETER, PRIVATE :: &
223  a(nellips)=(/ & ! semi-major axis for each ellipsoid
224  6378137.0d0, &
225  6378136.0d0, &
226  6378137.0d0, &
227  6378140.0d0, &
228  6377563.396d0, &
229  6378137.0d0, &
230  6378145.0d0, &
231  6377340.189d0, &
232  6377104.43d0, &
233  6378160.0d0, &
234  6378160.0d0, &
235  6377397.155d0, &
236  6377483.865d0, &
237  6378206.4d0, &
238  6378249.145d0, &
239  6375738.7d0, &
240  6376428.d0, &
241  6378136.05d0, &
242  6377276.345d0, &
243  6377304.063d0, &
244  6377301.243d0, &
245  6377295.664d0, &
246  6377298.556d0, &
247  6378166.d0, &
248  6378155.d0, &
249  6378150.d0, &
250  6378200.d0, &
251  6378270.0d0, &
252  6378388.0d0, &
253  6378245.0d0, &
254  6378163.d0, &
255  6378139.d0, &
256  6397300.d0, &
257  6378157.5d0, &
258  6376523.d0, &
259  6378155.0d0, &
260  6376896.0d0, &
261  6378165.0d0, &
262  6378145.0d0, &
263  6378135.0d0, &
264  6378137.0d0 /)
266 DOUBLE PRECISION,PARAMETER,PRIVATE :: k0=0.9996d0 ! scale factor at central meridian (check whether this is correct and constant)
267 
268 PRIVATE
269 PUBLIC geo_proj, geo_proj_rotated, geo_proj_stretched, geo_proj_polar, &
270  geo_proj_ellips, &
271  geo_proj_new, delete, copy, get_val, set_val, c_e, &
272  write_unit, read_unit, display, proj, unproj, operator(==), operator(/=)
273 
274 
275 CONTAINS
276 
283 FUNCTION 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)
288 CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
289 DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
290 INTEGER,INTENT(in),OPTIONAL :: zone
291 DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
292 DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
293 DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
294 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
295 DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
296 DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
297 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
298 DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
299 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
300 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
301 DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
302 INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
303 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
304 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
305 INTEGER,INTENT(in),OPTIONAL :: ellips_type
307 TYPE(geo_proj) :: this
308 
309 this%proj_type = optio_c(proj_type, len(this%proj_type))
310 
311 ! line of view / central meridian, use set_val
312 this%lov = optio_d(lov)
313 CALL set_val(this, zone=zone, lov=lov)
314 
315 ! offset / false *ing
316 this%xoff = optio_d(xoff)
317 IF (.NOT.c_e(this%xoff)) this%xoff = 0.0d0
318 this%yoff = optio_d(yoff)
319 IF (.NOT.c_e(this%yoff)) this%yoff = 0.0d0
320 
321 this%rotated%longitude_south_pole = optio_d(longitude_south_pole)
322 this%rotated%latitude_south_pole = optio_d(latitude_south_pole)
323 this%rotated%angle_rotation = optio_d(angle_rotation)
324 this%stretched%longitude_stretch_pole = optio_d(longitude_stretch_pole)
325 this%stretched%latitude_stretch_pole = optio_d(latitude_stretch_pole)
326 this%stretched%stretch_factor = optio_d(stretch_factor)
327 this%polar%latin1 = optio_d(latin1)
328 this%polar%latin2 = optio_d(latin2)
329 this%polar%lad = optio_d(lad)
330 this%polar%projection_center_flag = optio_l(projection_center_flag)
331 
332 ! ellipsoid, start from sphere, then use set_val
333 CALL ellips_compute(this%ellips)
334 CALL set_val(this, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
335  ellips_type=ellips_type)
336 
337 END FUNCTION geo_proj_new
338 
339 
340 ! compute constants related to the desired ellipsoid as a function of
341 ! semi-major axis and inverse of flattening
342 SUBROUTINE ellips_compute(this, a, f)
343 TYPE(geo_proj_ellips),INTENT(inout) :: this
344 DOUBLE PRECISION,INTENT(in),OPTIONAL :: a, f
345 
346 IF (present(a) .AND. present(f)) THEN
347  this%f = f
348  this%a = a
349 ELSE IF (present(a)) THEN ! parameters for a spherical Earth with given radius
350  this%f = 0.0d0
351  this%a = a
352 ELSE ! parameters for a standard spherical Earth
353  this%f = 0.0d0
354  this%a = rearth
355 ENDIF
356 
357 this%e2 = 2.0d0*this%f - this%f*this%f ! Eccentricity
358 this%e1 = this%f/(2.0d0 - this%f)
359 this%ep2 = this%e2/(1.0d0 - this%e2)
360 this%e11 = 3.0d0*this%e1/2.0d0 - 27.0d0*this%e1*this%e1*this%e1/32.0d0
361 this%e12 = 21.0d0*this%e1*this%e1/16.0d0 &
362  - 55.0d0*this%e1*this%e1*this%e1*this%e1/32.0d0
363 this%e13 = 151.0d0*this%e1*this%e1*this%e1/96.0d0
364 this%e14 = 1097.0d0*this%e1*this%e1*this%e1*this%e1/512.0d0
365 this%e4 = this%e2*this%e2
366 this%e6 = this%e2*this%e4
367 this%ef0 = this%a*(1.0d0 - 0.25d0*this%e2&
368  *(1.0d0+this%e2/16.0d0*(3.0d0 + 1.25d0*this%e2)))
369 this%ef1 = this%a*(0.375d0*this%e2 &
370  *(1.0d0 + 0.25d0*this%e2*(1.0d0 + 0.46875d0*this%e2)))
371 this%ef2 = this%a*(0.05859375d0*this%e2*this%e2*(1.0d0 + 0.75d0*this%e2))
372 this%ef3 = this%a*this%e2*this%e2*this%e2*35.0d0/3072.0d0
373 
374 END SUBROUTINE ellips_compute
375 
376 
378 SUBROUTINE geo_proj_delete(this)
379 TYPE(geo_proj),INTENT(inout) :: this
380 
381 this%proj_type = cmiss
382 this%lov = dmiss
383 this%xoff = dmiss
384 this%yoff = dmiss
385 this%rotated%longitude_south_pole = dmiss
386 this%rotated%latitude_south_pole = dmiss
387 this%rotated%angle_rotation = dmiss
388 this%stretched%longitude_stretch_pole = dmiss
389 this%stretched%latitude_stretch_pole = dmiss
390 this%stretched%stretch_factor = dmiss
391 this%polar%latin1 = dmiss
392 this%polar%latin2 = dmiss
393 this%polar%lad = dmiss
394 this%polar%projection_center_flag = imiss
395 
396 END SUBROUTINE geo_proj_delete
397 
398 
400 SUBROUTINE geo_proj_copy(this, that)
401 TYPE(geo_proj),INTENT(in) :: this
402 TYPE(geo_proj),INTENT(out) :: that
403 
404 that = this
405 
406 END SUBROUTINE geo_proj_copy
407 
408 
409 SUBROUTINE 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)
415 TYPE(geo_proj),INTENT(in) :: this
416 CHARACTER(len=*),OPTIONAL :: proj_type
417 DOUBLE PRECISION,OPTIONAL :: lov
418 INTEGER,OPTIONAL :: zone
419 DOUBLE PRECISION,OPTIONAL :: xoff
420 DOUBLE PRECISION,OPTIONAL :: yoff
421 DOUBLE PRECISION,OPTIONAL :: longitude_south_pole
422 DOUBLE PRECISION,OPTIONAL :: latitude_south_pole
423 DOUBLE PRECISION,OPTIONAL :: angle_rotation
424 DOUBLE PRECISION,OPTIONAL :: longitude_stretch_pole
425 DOUBLE PRECISION,OPTIONAL :: latitude_stretch_pole
426 DOUBLE PRECISION,OPTIONAL :: stretch_factor
427 DOUBLE PRECISION,OPTIONAL :: latin1
428 DOUBLE PRECISION,OPTIONAL :: latin2
429 DOUBLE PRECISION,OPTIONAL :: lad
430 INTEGER,OPTIONAL :: projection_center_flag
431 DOUBLE PRECISION,OPTIONAL :: ellips_smaj_axis
432 DOUBLE PRECISION,OPTIONAL :: ellips_flatt
433 INTEGER,OPTIONAL :: ellips_type
434 INTEGER,OPTIONAL :: unit
435 
436 INTEGER :: i
437 
438 IF (present(proj_type)) proj_type = this%proj_type
439 
440 IF (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
444 ELSE IF (present(lov)) THEN
445  lov = this%lov
446 ELSE IF (present(zone)) THEN
447  zone = nint((this%lov + 183.0d0)/6.0d0)
448 ENDIF
449 
450 IF (present(xoff)) xoff = this%xoff
451 IF (present(yoff)) yoff = this%yoff
452 IF (present(longitude_south_pole)) longitude_south_pole = this%rotated%longitude_south_pole
453 IF (present(latitude_south_pole)) latitude_south_pole = this%rotated%latitude_south_pole
454 IF (present(angle_rotation)) angle_rotation = this%rotated%angle_rotation
455 IF (present(longitude_stretch_pole)) longitude_stretch_pole = this%stretched%longitude_stretch_pole
456 IF (present(latitude_stretch_pole)) latitude_stretch_pole = this%stretched%latitude_stretch_pole
457 IF (present(stretch_factor)) stretch_factor = this%stretched%stretch_factor
458 IF (present(latin1)) latin1 = this%polar%latin1
459 IF (present(latin2)) latin2 = this%polar%latin2
460 IF (present(lad)) lad = this%polar%lad
461 IF (present(projection_center_flag)) projection_center_flag = this%polar%projection_center_flag
462 ! ellipsoid
463 IF (present(ellips_smaj_axis)) ellips_smaj_axis = this%ellips%a
464 IF (present(ellips_flatt)) ellips_flatt = this%ellips%f
465 IF (present(ellips_type)) THEN
466  ellips_type = imiss
467  DO i = 1, nellips
468  IF (this%ellips%f == 1.0d0/rf(i) .AND. this%ellips%a == a(i)) THEN
469  ellips_type = i
470  EXIT
471  ENDIF
472  ENDDO
473 ENDIF
474 
475 IF (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
481  CASE default
482  unit = imiss
483  END SELECT
484 ENDIF
485 
486 
487 END SUBROUTINE geo_proj_get_val
488 
489 
490 SUBROUTINE 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)
496 TYPE(geo_proj),INTENT(inout) :: this
497 CHARACTER(len=*),OPTIONAL :: proj_type
498 DOUBLE PRECISION,OPTIONAL :: lov
499 INTEGER,INTENT(in),OPTIONAL :: zone
500 DOUBLE PRECISION,OPTIONAL :: xoff
501 DOUBLE PRECISION,OPTIONAL :: yoff
502 DOUBLE PRECISION,OPTIONAL :: longitude_south_pole
503 DOUBLE PRECISION,OPTIONAL :: latitude_south_pole
504 DOUBLE PRECISION,OPTIONAL :: angle_rotation
505 DOUBLE PRECISION,OPTIONAL :: longitude_stretch_pole
506 DOUBLE PRECISION,OPTIONAL :: latitude_stretch_pole
507 DOUBLE PRECISION,OPTIONAL :: stretch_factor
508 DOUBLE PRECISION,OPTIONAL :: latin1
509 DOUBLE PRECISION,OPTIONAL :: latin2
510 DOUBLE PRECISION,OPTIONAL :: lad
511 INTEGER,OPTIONAL :: projection_center_flag
512 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
513 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
514 INTEGER,INTENT(in),OPTIONAL :: ellips_type
515 
516 INTEGER :: lzone
517 DOUBLE PRECISION :: llov
518 
519 
520 ! line of view / central meridian
521 llov = optio_d(lov)
522 lzone = optio_i(zone)
523 IF (c_e(llov) .AND. c_e(lzone)) THEN
524  this%lov = llov + zone*6.0d0 - 183.0d0
525 ELSE IF (c_e(llov)) THEN
526  this%lov = llov
527 ELSE IF (c_e(lzone)) THEN
528  this%lov = lzone*6.0d0 - 183.0d0
529 ENDIF
530 
531 ! ellipsoid
532 IF (present(ellips_smaj_axis)) THEN
533 ! explicit ellipsoid parameters provided (sphere if flatt is not present or 0)
534  CALL ellips_compute(this%ellips, ellips_smaj_axis, ellips_flatt)
535 ELSE IF (present(ellips_type)) THEN
536  IF (ellips_type > 0 .AND. ellips_type <= nellips) THEN
537 ! an hard coded ellipsoid has been requested
538  CALL ellips_compute(this%ellips, a(ellips_type), 1.0d0/rf(ellips_type))
539  ELSE ! fallback to default sphere
540  CALL ellips_compute(this%ellips)
541  ENDIF
542 ENDIF
543 
544 ! todo all the rest
545 
546 END SUBROUTINE geo_proj_set_val
547 
548 
549 FUNCTION geo_proj_c_e(this) RESULT(c_e)
550 TYPE(geo_proj),INTENT(in) :: this
551 LOGICAL :: c_e
552 
553 c_e = this%proj_type /= cmiss
554 
555 END FUNCTION geo_proj_c_e
556 
557 
562 SUBROUTINE geo_proj_read_unit(this, unit)
563 TYPE(geo_proj),INTENT(out) :: this
564 INTEGER, INTENT(in) :: unit
565 
566 CHARACTER(len=40) :: form
567 
568 INQUIRE(unit, form=form)
569 IF (form == 'FORMATTED') THEN
570  READ(unit,*)this
571 ELSE
572  READ(unit)this
573 ENDIF
574 
575 END SUBROUTINE geo_proj_read_unit
576 
577 
582 SUBROUTINE geo_proj_write_unit(this, unit)
583 TYPE(geo_proj),INTENT(in) :: this
584 INTEGER, INTENT(in) :: unit
586 CHARACTER(len=40) :: form
587 
588 INQUIRE(unit, form=form)
589 IF (form == 'FORMATTED') THEN
590  WRITE(unit,*)this
591 ELSE
592  WRITE(unit)this
593 ENDIF
595 END SUBROUTINE geo_proj_write_unit
596 
597 
598 SUBROUTINE geo_proj_display(this)
599 TYPE(geo_proj),INTENT(in) :: this
600 
601 print*,"<<<<<<<<<<<<<<< ",t2c(this%proj_type,"undefined projection"), &
602  " >>>>>>>>>>>>>>>>"
603 
604 IF (c_e(this%xoff) .AND. this%xoff /= 0.0d0) THEN
605  print*,"False easting",this%xoff
606 ENDIF
607 IF (c_e(this%yoff) .AND. this%yoff /= 0.0d0) THEN
608  print*,"False northing",this%yoff
609 ENDIF
610 
611 IF (c_e(this%lov)) THEN
612  print*,"Central Meridian",this%lov
613 ENDIF
614 
615 IF (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
620 ENDIF
621 
622 IF (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
627 ENDIF
628 
629 IF (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
633  ENDIF
634  IF (c_e(this%polar%lad)) THEN
635  print*,"isometric latitude",this%polar%lad
636  ENDIF
637  IF (iand(this%polar%projection_center_flag, 128) == 0) THEN
638  print*,"North Pole"
639  ELSE
640  print*,"South Pole"
641  ENDIF
642 ENDIF
643 
644 IF (this%proj_type == 'mercator') THEN
645  IF (c_e(this%polar%lad)) THEN
646  print*,"isometric latitude",this%polar%lad
647  ENDIF
648 ENDIF
649 
650 IF (this%ellips%f == 0.0d0) THEN
651  print*,"Spherical Earth:"
652  print*,"Radius (m)",this%ellips%a
653 ELSE
654  print*,"Ellipsoid:"
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
658 ENDIF
659 
660 
661 END SUBROUTINE geo_proj_display
662 
663 
666 ELEMENTAL SUBROUTINE geo_proj_proj(this, lon, lat, x, y)
667 TYPE(geo_proj),INTENT(in) :: this
669 DOUBLE PRECISION, INTENT(in) :: lon, lat
670 
671 DOUBLE PRECISION, INTENT(out) :: x, y
672 
673 SELECT CASE(this%proj_type)
674 
675 CASE("regular_ll")
676  CALL proj_regular_ll(lon, lat, x, y)
677 
678 CASE("rotated_ll")
679  CALL proj_rotated_ll(lon, lat, x, y, this%rotated%longitude_south_pole, &
680  this%rotated%latitude_south_pole, this%rotated%angle_rotation)
681 
682 CASE("lambert")
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)
686 
687 CASE("polar_stereographic")
688  CALL proj_polar_stereographic(lon, lat, x, y, this%lov, &
689  this%polar%lad, this%polar%projection_center_flag)
690 
691 CASE("mercator")
692  CALL proj_mercator(lon, lat, x, y, this%lov, this%polar%lad)
693 
694 CASE("UTM")
695  CALL proj_utm(lon, lat, x, y, this%lov, this%xoff, this%yoff, this%ellips)
696 
697 CASE default
698  x = dmiss
699  y = dmiss
700 
701 END SELECT
702 
703 END SUBROUTINE geo_proj_proj
704 
705 
708 ELEMENTAL SUBROUTINE geo_proj_unproj(this, x, y, lon, lat)
709 TYPE(geo_proj),INTENT(in) :: this
711 DOUBLE PRECISION, INTENT(in) :: x, y
712 
713 DOUBLE PRECISION, INTENT(out) :: lon, lat
714 
715 SELECT CASE(this%proj_type)
716 
717 CASE("regular_ll")
718  CALL unproj_regular_ll(x, y, lon, lat)
719 
720 CASE("rotated_ll")
721  CALL unproj_rotated_ll(x, y, lon, lat, this%rotated%longitude_south_pole, &
722  this%rotated%latitude_south_pole, this%rotated%angle_rotation)
723 
724 CASE("lambert")
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)
728 
729 CASE("polar_stereographic")
730  CALL unproj_polar_stereographic(x, y, lon, lat, this%lov, &
731  this%polar%lad, this%polar%projection_center_flag)
732 
733 CASE("mercator")
734  CALL unproj_mercator(x, y, lon, lat, this%lov, this%polar%lad)
735 
736 CASE("UTM")
737  CALL unproj_utm(x, y, lon, lat, this%lov, this%xoff, this%yoff, this%ellips)
738 
739 CASE default
740  lon = dmiss
741  lat = dmiss
742 
743 END SELECT
744 
745 END SUBROUTINE geo_proj_unproj
746 
748 ELEMENTAL FUNCTION geo_proj_eq(this, that) RESULT(eq)
749 TYPE(geo_proj),INTENT(in) :: this, that
750 LOGICAL :: eq
751 
752 eq = 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. & ! polar%lon1, polar%lat1
761  this%polar%latin2 == that%polar%latin2 .AND. & ! intentionally not checked
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
765 
766 END FUNCTION geo_proj_eq
768 
769 ELEMENTAL FUNCTION geo_proj_ne(this, that) RESULT(ne)
770 TYPE(geo_proj),INTENT(in) :: this, that
771 LOGICAL :: ne
772 
773 ne = .NOT. (this == that)
774 
775 END FUNCTION geo_proj_ne
776 
777 
780 ELEMENTAL FUNCTION geo_lon_eq(l1, l2) RESULT(eq)
781 DOUBLE PRECISION,INTENT(in) :: l1
782 DOUBLE PRECISION,INTENT(in) :: l2
783 
784 LOGICAL :: eq
785 
786 !eq = (l1 == l2) .OR. (ABS(l2-l1) == 360.0D0)
787 eq = modulo(l2-l1, 360.0d0) == 0.0d0
788 
789 END FUNCTION geo_lon_eq
790 
791 ! =====================
792 ! == transformations ==
793 ! =====================
794 
795 ELEMENTAL SUBROUTINE proj_regular_ll(lon,lat,x,y)
796 DOUBLE PRECISION, INTENT(in) :: lon,lat
797 DOUBLE PRECISION, INTENT(out) :: x,y
798 
799 x = lon
800 y = lat
801 
802 END SUBROUTINE proj_regular_ll
803 
804 ELEMENTAL SUBROUTINE unproj_regular_ll(x,y,lon,lat)
805 DOUBLE PRECISION, INTENT(in) :: x,y
806 DOUBLE PRECISION, INTENT(out) :: lon,lat
807 
808 lon = x
809 lat = y
810 
811 END SUBROUTINE unproj_regular_ll
812 
813 
814 ELEMENTAL SUBROUTINE proj_rotated_ll(lon,lat,x,y, &
815  longitude_south_pole, latitude_south_pole, angle_rotation)
816 DOUBLE PRECISION, INTENT(in) :: lon,lat
817 DOUBLE PRECISION, INTENT(out) :: x,y
818 DOUBLE PRECISION, INTENT(in) :: longitude_south_pole, latitude_south_pole, &
819  angle_rotation
820 
821 DOUBLE PRECISION :: cy0,sy0,rx,srx,crx,sy,cy,l_south_pole
822 
823 ! old hibu formula
824 !l_south_pole = ACOS(-SIN(degrad*latitude_south_pole))
825 l_south_pole = (latitude_south_pole+90.)*degrad
826 
827 rx = degrad*(lon - longitude_south_pole)
828 srx = sin(rx)
829 crx = cos(rx)
830 
831 sy0 = sin(l_south_pole)
832 cy0 = cos(l_south_pole)
833 
834 sy = sin(degrad*lat)
835 cy = cos(degrad*lat)
836 
837 x = raddeg*atan2(cy*srx, cy0*cy*crx+sy0*sy)
838 y = raddeg*asin(cy0*sy - sy0*cy*crx)
839 x = x + angle_rotation ! check
840 
841 END SUBROUTINE proj_rotated_ll
842 
843 ELEMENTAL SUBROUTINE unproj_rotated_ll(x,y,lon,lat,&
844  longitude_south_pole, latitude_south_pole, angle_rotation)
845 DOUBLE PRECISION, INTENT(in) :: x,y
846 DOUBLE PRECISION, INTENT(out) :: lon,lat
847 DOUBLE PRECISION, INTENT(in) :: longitude_south_pole, latitude_south_pole, &
848  angle_rotation
849 
850 DOUBLE PRECISION :: cy0, sy0, l_south_pole, xr
852 xr = (x - angle_rotation)*degrad ! check
853 ! old hibu formula
854 !l_south_pole = ACOS(-SIN(degrad*latitude_south_pole))
855 l_south_pole = (latitude_south_pole+90.)*degrad
856 
857 cy0 = cos(l_south_pole)
858 sy0 = sin(l_south_pole)
859 
860 lat = raddeg*asin(sy0*cos(degrad*y)*cos(xr)+cy0*sin(degrad*y))
861 lon = longitude_south_pole + &
862  raddeg*asin(sin(xr)*cos(degrad*y)/cos(degrad*lat))
863 
864 END SUBROUTINE unproj_rotated_ll
865 
866 ! come usare il polo? ruotare e antiruotare?
867 ELEMENTAL SUBROUTINE proj_stretched_ll(lon,lat,x,y, &
868  longitude_stretch_pole, latitude_stretch_pole, stretch_factor)
869 DOUBLE PRECISION, INTENT(in) :: lon,lat
870 DOUBLE PRECISION, INTENT(out) :: x,y
871 DOUBLE PRECISION, INTENT(in) :: longitude_stretch_pole, latitude_stretch_pole, &
872  stretch_factor
873 
874 DOUBLE PRECISION :: csq
875 
876 csq = stretch_factor**2
877 x = lon
878 y = raddeg*asin((1.0d0 - csq + (1.0d0 + csq)*sin(degrad*lat)) / &
879  (1.0d0 + csq + (1.0d0 - csq)*sin(degrad*lat)))
880 
881 END SUBROUTINE proj_stretched_ll
882 
883 ELEMENTAL SUBROUTINE unproj_stretched_ll(x,y,lon,lat,&
884  longitude_stretch_pole, latitude_stretch_pole, stretch_factor)
885 DOUBLE PRECISION, INTENT(in) :: x,y
886 DOUBLE PRECISION, INTENT(out) :: lon,lat
887 DOUBLE PRECISION, INTENT(in) :: longitude_stretch_pole, latitude_stretch_pole, &
888  stretch_factor
889 
890 DOUBLE PRECISION :: csq
891 
892 csq = stretch_factor**2
893 lon = x
894 ! TODO verificare la formula inversa
895 lat = raddeg*asin((csq - 1.0d0 + (csq + 1.0d0)*sin(degrad*y)) / &
896  (csq + 1.0d0 + (csq - 1.0d0)*sin(degrad*y)))
897 
898 END SUBROUTINE unproj_stretched_ll
899 
900 
901 ! Formulas and notation from:
902 ! http://mathworld.wolfram.com/LambertConformalConicProjection.html
903 ! http://en.wikipedia.org/wiki/Lambert_conformal_conic_projection
904 ! http://fr.wikipedia.org/wiki/Projection_conique_conforme_de_Lambert
905 ! with the following guess:
906 ! projection is always polar, so reference latitude=+-90 according to
907 ! projectionCenterFlag; reference longitude is LoV.
908 ! how coordinates of south pole should be treated? Metview ignores them.
909 ELEMENTAL SUBROUTINE proj_lambert(lon,lat,x,y, &
910  latin1, latin2, lov, lad, projection_center_flag)
911 DOUBLE PRECISION, INTENT(in) :: lon,lat
912 DOUBLE PRECISION, INTENT(out) :: x,y
913 DOUBLE PRECISION, INTENT(in) :: latin1, latin2, lov, lad
914 INTEGER, INTENT(in) :: projection_center_flag
915 
916 DOUBLE PRECISION :: n, f, ro0, ro, cs1, cs2, cs3, pollat, angle, cot
917 DOUBLE PRECISION, PARAMETER :: epsy = 1.0d-100
918 
919 IF (iand(projection_center_flag, 128) == 0) THEN
920  pollat = 90.d0*degrad
921 ELSE
922  pollat = -90.d0*degrad
923 ENDIF
924 cs1 = cos(degrad*latin1)
925 cs2 = tan(pi*.25d0 + degrad*latin1*.5d0)
926 
927 IF (latin1 == latin2) THEN
928  n = sin(degrad*latin1) ! verify that n->sin(latin1) when latin2->latin1
929 ELSE
930  n = log(cs1/cos(degrad*latin2)) / &
931  log(tan(pi*.25d0 + degrad*latin2*.5d0) / cs2)
932 ENDIF
933 f = cs1*cs2**n/n*rearth ! check that rearth is correct here (only if lad==latin1?)
934 angle = pi*.25d0 + pollat*.5d0
935 cot = cos(angle)/sin(angle)
936 IF (cot > epsy) THEN
937  ro0 = f*cot**n
938 ELSE
939  ro0 = 0.0d0
940 ENDIF
941 
942 angle = pi*.25d0 + degrad*lat*.5d0
943 cot = cos(angle)/sin(angle)
944 IF (cot > epsy) THEN
945  ro = f*cot**n
946 ELSE
947  ro = 0.0d0
948 ENDIF
949 
950 cs3 = degrad*n*(lon - lov)
951 
952 x = ro*sin(cs3)
953 y = ro0 - ro*cos(cs3)
954 
955 END SUBROUTINE proj_lambert
956 
957 ELEMENTAL SUBROUTINE unproj_lambert(x,y,lon,lat, &
958  latin1, latin2, lov, lad, projection_center_flag)
959 DOUBLE PRECISION, INTENT(in) :: x,y
960 DOUBLE PRECISION, INTENT(out) :: lon,lat
961 DOUBLE PRECISION, INTENT(in) :: latin1, latin2, lov, lad
962 INTEGER, INTENT(in) :: projection_center_flag
963 
964 DOUBLE PRECISION :: n, f, ro0, ro, theta, cs1, cs2, pollat, angle, cot
965 DOUBLE PRECISION, PARAMETER :: epsy = 1.0d-100
966 
967 ! check, pollat is actually used as the latitude at which
968 ! y=0, may be not correct and is not enough for Southern Hemisphere
969 IF (iand(projection_center_flag, 128) == 0) THEN
970  pollat = 90.d0*degrad
971 ELSE
972  pollat = -90.d0*degrad
973 ENDIF
974 cs1 = cos(degrad*latin1)
975 cs2 = tan(pi*.25d0 + degrad*latin1*.5d0)
976 
977 IF (latin1 == latin2) THEN
978  n = sin(degrad*latin1) ! verify limit
979 ELSE
980  n = log(cs1/cos(degrad*latin2)) / &
981  log(tan(pi*.25d0 + degrad*latin2*.5d0) / cs2)
982 ENDIF
983 f = cs1*cs2**n/n*rearth ! check that rearth is correct here (only if lad==latin1?)
984 angle = pi*.25d0 + pollat*.5d0
985 cot = cos(angle)/sin(angle)
986 IF (cot > epsy) THEN
987  ro0 = f*cot**n
988 ELSE
989  ro0 = 0.0d0
990 ENDIF
991 
992 ro = sign(sqrt(x*x + (ro0-y)*(ro0-y)), n) ! check SIGN
993 theta = raddeg*atan2(x, ro0-y)
994 
995 lon = lov + theta/n
996 lat = raddeg*(2.d0*atan((f/ro)**(1.d0/n)) - pi*.5d0)
997 
998 END SUBROUTINE unproj_lambert
999 
1000 
1001 !http://mathworld.wolfram.com/StereographicProjection.html
1002 ELEMENTAL SUBROUTINE proj_polar_stereographic(lon,lat,x,y, &
1003  lov, lad, projection_center_flag)
1004 DOUBLE PRECISION, INTENT(in) :: lon,lat
1005 DOUBLE PRECISION, INTENT(out) :: x,y
1006 DOUBLE PRECISION, INTENT(in) :: lov, lad
1007 INTEGER, INTENT(in) :: projection_center_flag
1008 
1009 DOUBLE PRECISION :: k, pollat
1010 
1011 IF (iand(projection_center_flag, 128) == 0) THEN
1012  pollat = 90.d0*degrad
1013 ELSE
1014  pollat = -90.d0*degrad
1015 ENDIF
1016 
1017 k = 2.0d0*rearth/(1.0d0 + sin(pollat)*sin(degrad*lat) + &
1018  cos(pollat)*cos(degrad*lat)*cos(degrad*(lon - lov)))
1019 x = k*cos(degrad*lat)*sin(degrad*(lon - lov))
1020 y = k*(cos(pollat)*sin(degrad*lat) - &
1021  sin(pollat)*cos(degrad*lat)*cos(degrad*(lon - lov)))
1022 
1023 END SUBROUTINE proj_polar_stereographic
1024 
1025 ELEMENTAL SUBROUTINE unproj_polar_stereographic(x,y,lon,lat, &
1026  lov, lad, projection_center_flag)
1027 DOUBLE PRECISION, INTENT(in) :: x,y
1028 DOUBLE PRECISION, INTENT(out) :: lon,lat
1029 DOUBLE PRECISION, INTENT(in) :: lov, lad
1030 INTEGER, INTENT(in) :: projection_center_flag
1031 
1032 DOUBLE PRECISION :: ro, c, pollat
1033 
1034 IF (iand(projection_center_flag, 128) == 0) THEN
1035  pollat = 90.d0*degrad
1036 ELSE
1037  pollat = -90.d0*degrad
1038 ENDIF
1039 
1040 ro = sqrt(x**2 + y**2)
1041 c = 2.0d0*atan(ro/(2.0d0*rearth))
1042 lat = raddeg*asin(cos(c)*sin(pollat)+y*sin(c)*cos(pollat)/ro)
1043 lon = lov + raddeg*atan2(x*sin(c), &
1044  (ro*cos(pollat)*cos(c)-y*sin(pollat)*sin(c)))
1045 
1046 END SUBROUTINE unproj_polar_stereographic
1047 
1048 
1049 !http://mathworld.wolfram.com/StereographicProjection.html
1050 ELEMENTAL SUBROUTINE proj_mercator(lon, lat, x, y, lov, lad)
1051 DOUBLE PRECISION,INTENT(in) :: lon,lat
1052 DOUBLE PRECISION,INTENT(out) :: x,y
1053 DOUBLE PRECISION,INTENT(in) :: lov
1054 DOUBLE PRECISION,INTENT(in) :: lad
1055 
1056 DOUBLE PRECISION :: scx, scy
1057 
1058 scy = cos(degrad*lad)*rearth
1059 scx = scy*degrad
1060 x = (lon-lov)*scx
1061 y = log(tan(0.25d0*pi + 0.5d0*lat*degrad))*scy
1062 
1063 END SUBROUTINE proj_mercator
1064 
1065 ELEMENTAL SUBROUTINE unproj_mercator(x, y, lon, lat, lov, lad)
1066 DOUBLE PRECISION,INTENT(in) :: x,y
1067 DOUBLE PRECISION,INTENT(out) :: lon,lat
1068 DOUBLE PRECISION,INTENT(in) :: lov
1069 DOUBLE PRECISION,INTENT(in) :: lad
1070 
1071 DOUBLE PRECISION :: scx, scy
1072 
1073 scy = cos(degrad*lad)*rearth
1074 scx = scy*degrad
1075 lon = x/scx + lov
1076 lat = 2.0d0*atan(exp(y/scy))*raddeg-90.0d0
1077 
1078 END SUBROUTINE unproj_mercator
1079 
1080 
1081 ELEMENTAL SUBROUTINE proj_utm(lon, lat, x, y, lov, false_e, false_n, ellips)
1082 DOUBLE PRECISION,INTENT(in) :: lon,lat
1083 DOUBLE PRECISION,INTENT(out) :: x,y
1084 DOUBLE PRECISION,INTENT(in) :: lov
1085 DOUBLE PRECISION,INTENT(in) :: false_e, false_n
1086 TYPE(geo_proj_ellips),INTENT(in) :: ellips
1087 
1088 DOUBLE PRECISION :: deltalon, p
1089 DOUBLE PRECISION :: n, t, t2, c, m, a1, a2, a3, a4, a5, a6, sinp, cosp, tanp
1090 
1091 ! --- Compute delta longitude in radians
1092 deltalon = degrad*(lon - lov)
1093 
1094 ! --- Convert phi (latitude) to radians
1095 p = degrad*lat
1096 sinp = sin(p)
1097 cosp = cos(p)
1098 tanp = tan(p)
1099 
1100 n = ellips%a/sqrt(1.0d0 - ellips%e2*sinp*sinp)
1101 t = tanp*tanp
1102 c = ellips%ep2*cosp*cosp
1103 a1 = deltalon*cosp
1104 !!$m = 111132.0894_fp_utm*lat - 16216.94_fp_utm*SIN(2.0*p) + 17.21_fp_utm*SIN(4.0*p) &
1105 !!$ - 0.02_fp_utm*SIN(6.0*p)
1106 ! Modificato rispetto alla routine originale, dipende dall'ellissoide
1107 m = ellips%ef0*p - ellips%ef1*sin(2.0d0*p) + ellips%ef2*sin(4.0d0*p) - &
1108  ellips%ef3*sin(6.0d0*p)
1109 
1110 a2 = a1**2
1111 a3 = a2*a1
1112 a4 = a2**2
1113 a5 = a4*a1
1114 a6 = a4*a2
1115 t2 = t**2
1116 
1117 ! --- Compute UTM x and y (m)
1118 x = 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
1121 y = k0*(m + n*tanp &
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))
1125 y = y + false_n
1126 
1127 END SUBROUTINE proj_utm
1128 
1129 ELEMENTAL SUBROUTINE unproj_utm(x, y, lon, lat, lov, false_e, false_n, ellips)
1130 DOUBLE PRECISION,INTENT(IN) :: x, y
1131 DOUBLE PRECISION,INTENT(OUT) :: lon, lat
1132 DOUBLE PRECISION,INTENT(in) :: lov
1133 DOUBLE PRECISION,INTENT(in) :: false_e, false_n
1134 TYPE(geo_proj_ellips),INTENT(in) :: ellips
1135 
1136 DOUBLE 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
1138 
1139 ! --- Correct for false easting, southern hemisphere
1140 xm = x - false_e
1141 !IF (zone < 0) THEN
1142 ym = y - false_n
1143 !ELSE
1144 !ym = utmn
1145 !ENDIF
1146 
1147 m = ym/k0
1148 u = m/(ellips%a*(1.0d0-ellips%e2/4.0d0 - &
1149  3.0d0*ellips%e4/64.0d0 - 5.0d0*ellips%e6/256.0d0))
1150 p1 = 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)
1152 sinp1 = sin(p1)
1153 cosp1 = cos(p1)
1154 tanp1 = tan(p1)
1155 c1 = ellips%ep2*cosp1**2
1156 c2 = c1**2
1157 t1 = tanp1**2
1158 t2 = t1**2
1159 sin2p1 = sinp1**2
1160 r0 = 1.0d0-ellips%e2*sin2p1
1161 n1 = ellips%a/sqrt(r0)
1162 !r1 = ellips%a*(1.0D0-ellips%e2)/SQRT(r0**3)
1163 r1 = r0/(1.0d0-ellips%e2)
1164 
1165 d = xm/(n1*k0)
1166 d2=d**2
1167 d3=d*d2
1168 d4=d*d3
1169 d5=d*d4
1170 d6=d*d5
1171 
1172 ! other variant from mstortini:
1173 !p = p1 - (1.0D0*tanp1/r0) * (d2/2.0D0 &
1174 ! original version with r1 = ellips%a*(1.0D0-ellips%e2)/SQRT(r0**3)
1175 !p = p1 - (n1*tanp1/r1) * (d2/2.0D0 &
1176 ! optimized version with different definition of r1
1177 p = 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)
1182 lat = raddeg*p
1183 l = (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
1186 lon = raddeg*l + lov
1187 
1188 END SUBROUTINE unproj_utm
1189 
1190 END MODULE geo_proj_class
Write the object on a formatted or unformatted file.
Read the object from a formatted or unformatted file.
Definitions of constants and functions for working with missing values.
Functions that return a trimmed CHARACTER representation of the input variable.
Compute forward coordinate transformation from geographical system to projected system.
Operatore di valore assoluto di un intervallo.
Compute backward coordinate transformation from projected system to geographical system.
Copy an object, creating a fully new instance.
Print a brief description on stdout.
Method for returning the contents of the object.
Destructors of the corresponding objects.
Utilities for CHARACTER variables.
Method for testing the existence of the object.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Method for setting the contents of the object.
Costanti fisiche (DOUBLEPRECISION).
Definition: phys_const.f90:72

Generated with Doxygen.