libsim  Versione6.3.0
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 /)
265 
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 
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
306 
307 TYPE(geo_proj) :: this
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
585 
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  " >>>>>>>>>>>>>>>>"
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%ellips%f == 0.0d0) THEN
645  print*,"Spherical Earth:"
646  print*,"Radius (m)",this%ellips%a
647 ELSE
648  print*,"Ellipsoid:"
649  print*,"Flattening",this%ellips%f
650  print*,"Reverse of flattening",1.0d0/this%ellips%f
651  print*,"Semi-major axis (m)",this%ellips%a
652 ENDIF
653 
654 
655 END SUBROUTINE geo_proj_display
656 
657 
660 ELEMENTAL SUBROUTINE geo_proj_proj(this, lon, lat, x, y)
661 TYPE(geo_proj),INTENT(in) :: this
663 DOUBLE PRECISION, INTENT(in) :: lon, lat
664 
665 DOUBLE PRECISION, INTENT(out) :: x, y
666 
667 SELECT CASE(this%proj_type)
668 
669 CASE("regular_ll")
670  CALL proj_regular_ll(lon, lat, x, y)
671 
672 CASE("rotated_ll")
673  CALL proj_rotated_ll(lon, lat, x, y, this%rotated%longitude_south_pole, &
674  this%rotated%latitude_south_pole, this%rotated%angle_rotation)
675 
676 CASE("lambert")
677  CALL proj_lambert(lon, lat, x, y, this%polar%latin1, &
678  this%polar%latin2, this%lov, this%polar%lad, &
679  this%polar%projection_center_flag)
680 
681 CASE("polar_stereographic")
682  CALL proj_polar_stereographic(lon, lat, x, y, this%lov, &
683  this%polar%lad, this%polar%projection_center_flag)
684 
685 CASE("UTM")
686  CALL proj_utm(lon, lat, x, y, this%lov, this%xoff, this%yoff, this%ellips)
687 
688 CASE default
689  x = dmiss
690  y = dmiss
691 
692 END SELECT
693 
694 END SUBROUTINE geo_proj_proj
695 
696 
699 ELEMENTAL SUBROUTINE geo_proj_unproj(this, x, y, lon, lat)
700 TYPE(geo_proj),INTENT(in) :: this
702 DOUBLE PRECISION, INTENT(in) :: x, y
703 
704 DOUBLE PRECISION, INTENT(out) :: lon, lat
705 
706 SELECT CASE(this%proj_type)
707 
708 CASE("regular_ll")
709  CALL unproj_regular_ll(x, y, lon, lat)
710 
711 CASE("rotated_ll")
712  CALL unproj_rotated_ll(x, y, lon, lat, this%rotated%longitude_south_pole, &
713  this%rotated%latitude_south_pole, this%rotated%angle_rotation)
714 
715 CASE("lambert")
716  CALL unproj_lambert(x, y, lon, lat, this%polar%latin1, &
717  this%polar%latin2, this%lov, this%polar%lad, &
718  this%polar%projection_center_flag)
719 
720 CASE("polar_stereographic")
721  CALL unproj_polar_stereographic(x, y, lon, lat, this%lov, &
722  this%polar%lad, this%polar%projection_center_flag)
723 
724 CASE("UTM")
725  CALL unproj_utm(x, y, lon, lat, this%lov, this%xoff, this%yoff, this%ellips)
726 
727 CASE default
728  lon = dmiss
729  lat = dmiss
730 
731 END SELECT
732 
733 END SUBROUTINE geo_proj_unproj
734 
735 
736 ELEMENTAL FUNCTION geo_proj_eq(this, that) RESULT(eq)
737 TYPE(geo_proj),INTENT(in) :: this, that
738 LOGICAL :: eq
739 
740 eq = this%proj_type == that%proj_type .AND. this%xoff == that%xoff .AND. &
741  this%yoff == that%yoff .AND. geo_lon_eq(this%lov, that%lov) .AND. &
742  geo_lon_eq(this%rotated%longitude_south_pole, that%rotated%longitude_south_pole) .AND. &
743  this%rotated%latitude_south_pole == that%rotated%latitude_south_pole .AND. &
744  this%rotated%angle_rotation == that%rotated%angle_rotation .AND. &
745  this%stretched%latitude_stretch_pole == that%stretched%latitude_stretch_pole .AND. &
746  geo_lon_eq(this%stretched%longitude_stretch_pole, that%stretched%longitude_stretch_pole) .AND. &
747  this%stretched%stretch_factor == that%stretched%stretch_factor .AND. &
748  this%polar%latin1 == that%polar%latin1 .AND. & ! polar%lon1, polar%lat1
749  this%polar%latin2 == that%polar%latin2 .AND. & ! intentionally not checked
750  this%polar%lad == that%polar%lad .AND. &
751  this%polar%projection_center_flag == that%polar%projection_center_flag .AND. &
752  this%ellips%f == that%ellips%f .AND. this%ellips%a == that%ellips%a
753 
754 END FUNCTION geo_proj_eq
755 
757 ELEMENTAL FUNCTION geo_proj_ne(this, that) RESULT(ne)
758 TYPE(geo_proj),INTENT(in) :: this, that
759 LOGICAL :: ne
760 
761 ne = .NOT. (this == that)
762 
763 END FUNCTION geo_proj_ne
764 
765 
768 ELEMENTAL FUNCTION geo_lon_eq(l1, l2) RESULT(eq)
769 DOUBLE PRECISION,INTENT(in) :: l1
770 DOUBLE PRECISION,INTENT(in) :: l2
771 
772 LOGICAL :: eq
773 
774 !eq = (l1 == l2) .OR. (ABS(l2-l1) == 360.0D0)
775 eq = modulo(l2-l1, 360.0d0) == 0.0d0
777 END FUNCTION geo_lon_eq
778 
779 ! =====================
780 ! == transformations ==
781 ! =====================
782 
783 ELEMENTAL SUBROUTINE proj_regular_ll(lon,lat,x,y)
784 DOUBLE PRECISION, INTENT(in) :: lon,lat
785 DOUBLE PRECISION, INTENT(out) :: x,y
786 
787 x = lon
788 y = lat
789 
790 END SUBROUTINE proj_regular_ll
791 
792 ELEMENTAL SUBROUTINE unproj_regular_ll(x,y,lon,lat)
793 DOUBLE PRECISION, INTENT(in) :: x,y
794 DOUBLE PRECISION, INTENT(out) :: lon,lat
795 
796 lon = x
797 lat = y
798 
799 END SUBROUTINE unproj_regular_ll
800 
801 
802 ELEMENTAL SUBROUTINE proj_rotated_ll(lon,lat,x,y, &
803  longitude_south_pole, latitude_south_pole, angle_rotation)
804 DOUBLE PRECISION, INTENT(in) :: lon,lat
805 DOUBLE PRECISION, INTENT(out) :: x,y
806 DOUBLE PRECISION, INTENT(in) :: longitude_south_pole, latitude_south_pole, &
807  angle_rotation
808 
809 DOUBLE PRECISION :: cy0,sy0,rx,srx,crx,sy,cy,l_south_pole
810 
811 ! old hibu formula
812 !l_south_pole = ACOS(-SIN(degrad*latitude_south_pole))
813 l_south_pole = (latitude_south_pole+90.)*degrad
814 
815 rx = degrad*(lon - longitude_south_pole)
816 srx = sin(rx)
817 crx = cos(rx)
818 
819 sy0 = sin(l_south_pole)
820 cy0 = cos(l_south_pole)
821 
822 sy = sin(degrad*lat)
823 cy = cos(degrad*lat)
824 
825 x = raddeg*atan2(cy*srx, cy0*cy*crx+sy0*sy)
826 y = raddeg*asin(cy0*sy - sy0*cy*crx)
827 x = x + angle_rotation ! check
828 
829 END SUBROUTINE proj_rotated_ll
830 
831 ELEMENTAL SUBROUTINE unproj_rotated_ll(x,y,lon,lat,&
832  longitude_south_pole, latitude_south_pole, angle_rotation)
833 DOUBLE PRECISION, INTENT(in) :: x,y
834 DOUBLE PRECISION, INTENT(out) :: lon,lat
835 DOUBLE PRECISION, INTENT(in) :: longitude_south_pole, latitude_south_pole, &
836  angle_rotation
837 
838 DOUBLE PRECISION :: cy0, sy0, l_south_pole, xr
839 
840 xr = (x - angle_rotation)*degrad ! check
841 ! old hibu formula
842 !l_south_pole = ACOS(-SIN(degrad*latitude_south_pole))
843 l_south_pole = (latitude_south_pole+90.)*degrad
844 
845 cy0 = cos(l_south_pole)
846 sy0 = sin(l_south_pole)
847 
848 lat = raddeg*asin(sy0*cos(degrad*y)*cos(xr)+cy0*sin(degrad*y))
849 lon = longitude_south_pole + &
850  raddeg*asin(sin(xr)*cos(degrad*y)/cos(degrad*lat))
851 
852 END SUBROUTINE unproj_rotated_ll
853 
854 ! come usare il polo? ruotare e antiruotare?
855 ELEMENTAL SUBROUTINE proj_stretched_ll(lon,lat,x,y, &
856  longitude_stretch_pole, latitude_stretch_pole, stretch_factor)
857 DOUBLE PRECISION, INTENT(in) :: lon,lat
858 DOUBLE PRECISION, INTENT(out) :: x,y
859 DOUBLE PRECISION, INTENT(in) :: longitude_stretch_pole, latitude_stretch_pole, &
860  stretch_factor
861 
862 DOUBLE PRECISION :: csq
863 
864 csq = stretch_factor**2
865 x = lon
866 y = raddeg*asin((1.0d0 - csq + (1.0d0 + csq)*sin(degrad*lat)) / &
867  (1.0d0 + csq + (1.0d0 - csq)*sin(degrad*lat)))
868 
869 END SUBROUTINE proj_stretched_ll
870 
871 ELEMENTAL SUBROUTINE unproj_stretched_ll(x,y,lon,lat,&
872  longitude_stretch_pole, latitude_stretch_pole, stretch_factor)
873 DOUBLE PRECISION, INTENT(in) :: x,y
874 DOUBLE PRECISION, INTENT(out) :: lon,lat
875 DOUBLE PRECISION, INTENT(in) :: longitude_stretch_pole, latitude_stretch_pole, &
876  stretch_factor
877 
878 DOUBLE PRECISION :: csq
879 
880 csq = stretch_factor**2
881 lon = x
882 ! TODO verificare la formula inversa
883 lat = raddeg*asin((csq - 1.0d0 + (csq + 1.0d0)*sin(degrad*y)) / &
884  (csq + 1.0d0 + (csq - 1.0d0)*sin(degrad*y)))
885 
886 END SUBROUTINE unproj_stretched_ll
887 
888 
889 ! Formulas and notation from:
890 ! http://mathworld.wolfram.com/LambertConformalConicProjection.html
891 ! http://en.wikipedia.org/wiki/Lambert_conformal_conic_projection
892 ! http://fr.wikipedia.org/wiki/Projection_conique_conforme_de_Lambert
893 ! with the following guess:
894 ! projection is always polar, so reference latitude=+-90 according to
895 ! projectionCenterFlag; reference longitude is LoV.
896 ! how coordinates of south pole should be treated? Metview ignores them.
897 ELEMENTAL SUBROUTINE proj_lambert(lon,lat,x,y, &
898  latin1, latin2, lov, lad, projection_center_flag)
899 DOUBLE PRECISION, INTENT(in) :: lon,lat
900 DOUBLE PRECISION, INTENT(out) :: x,y
901 DOUBLE PRECISION, INTENT(in) :: latin1, latin2, lov, lad
902 INTEGER, INTENT(in) :: projection_center_flag
903 
904 DOUBLE PRECISION :: n, f, ro0, ro, cs1, cs2, cs3, pollat, angle, cot
905 DOUBLE PRECISION, PARAMETER :: epsy = 1.0d-100
906 
907 IF (iand(projection_center_flag, 128) == 0) THEN
908  pollat = 90.d0*degrad
909 ELSE
910  pollat = -90.d0*degrad
911 ENDIF
912 cs1 = cos(degrad*latin1)
913 cs2 = tan(pi*.25d0 + degrad*latin1*.5d0)
914 
915 IF (latin1 == latin2) THEN
916  n = sin(degrad*latin1) ! verify that n->sin(latin1) when latin2->latin1
917 ELSE
918  n = log(cs1/cos(degrad*latin2)) / &
919  log(tan(pi*.25d0 + degrad*latin2*.5d0) / cs2)
920 ENDIF
921 f = cs1*cs2**n/n*rearth ! check that rearth is correct here (only if lad==latin1?)
922 angle = pi*.25d0 + pollat*.5d0
923 cot = cos(angle)/sin(angle)
924 IF (cot > epsy) THEN
925  ro0 = f*cot**n
926 ELSE
927  ro0 = 0.0d0
928 ENDIF
929 
930 angle = pi*.25d0 + degrad*lat*.5d0
931 cot = cos(angle)/sin(angle)
932 IF (cot > epsy) THEN
933  ro = f*cot**n
934 ELSE
935  ro = 0.0d0
936 ENDIF
937 
938 cs3 = degrad*n*(lon - lov)
939 
940 x = ro*sin(cs3)
941 y = ro0 - ro*cos(cs3)
942 
943 END SUBROUTINE proj_lambert
944 
945 ELEMENTAL SUBROUTINE unproj_lambert(x,y,lon,lat, &
946  latin1, latin2, lov, lad, projection_center_flag)
947 DOUBLE PRECISION, INTENT(in) :: x,y
948 DOUBLE PRECISION, INTENT(out) :: lon,lat
949 DOUBLE PRECISION, INTENT(in) :: latin1, latin2, lov, lad
950 INTEGER, INTENT(in) :: projection_center_flag
951 
952 DOUBLE PRECISION :: n, f, ro0, ro, theta, cs1, cs2, pollat, angle, cot
953 DOUBLE PRECISION, PARAMETER :: epsy = 1.0d-100
954 
955 ! check, pollat is actually used as the latitude at which
956 ! y=0, may be not correct and is not enough for Southern Hemisphere
957 IF (iand(projection_center_flag, 128) == 0) THEN
958  pollat = 90.d0*degrad
959 ELSE
960  pollat = -90.d0*degrad
961 ENDIF
962 cs1 = cos(degrad*latin1)
963 cs2 = tan(pi*.25d0 + degrad*latin1*.5d0)
964 
965 IF (latin1 == latin2) THEN
966  n = sin(degrad*latin1) ! verify limit
967 ELSE
968  n = log(cs1/cos(degrad*latin2)) / &
969  log(tan(pi*.25d0 + degrad*latin2*.5d0) / cs2)
970 ENDIF
971 f = cs1*cs2**n/n*rearth ! check that rearth is correct here (only if lad==latin1?)
972 angle = pi*.25d0 + pollat*.5d0
973 cot = cos(angle)/sin(angle)
974 IF (cot > epsy) THEN
975  ro0 = f*cot**n
976 ELSE
977  ro0 = 0.0d0
978 ENDIF
979 
980 ro = sign(sqrt(x*x + (ro0-y)*(ro0-y)), n) ! check SIGN
981 theta = raddeg*atan2(x, ro0-y)
982 
983 lon = lov + theta/n
984 lat = raddeg*(2.d0*atan((f/ro)**(1.d0/n)) - pi*.5d0)
985 
986 END SUBROUTINE unproj_lambert
987 
988 
989 !http://mathworld.wolfram.com/StereographicProjection.html
990 ELEMENTAL SUBROUTINE proj_polar_stereographic(lon,lat,x,y, &
991  lov, lad, projection_center_flag)
992 DOUBLE PRECISION, INTENT(in) :: lon,lat
993 DOUBLE PRECISION, INTENT(out) :: x,y
994 DOUBLE PRECISION, INTENT(in) :: lov, lad
995 INTEGER, INTENT(in) :: projection_center_flag
996 
997 DOUBLE PRECISION :: k, pollat
998 
999 IF (iand(projection_center_flag, 128) == 0) THEN
1000  pollat = 90.d0*degrad
1001 ELSE
1002  pollat = -90.d0*degrad
1003 ENDIF
1004 
1005 k = 2.0d0*rearth/(1.0d0 + sin(pollat)*sin(degrad*lat) + &
1006  cos(pollat)*cos(degrad*lat)*cos(degrad*(lon - lov)))
1007 x = k*cos(degrad*lat)*sin(degrad*(lon - lov))
1008 y = k*(cos(pollat)*sin(degrad*lat) - &
1009  sin(pollat)*cos(degrad*lat)*cos(degrad*(lon - lov)))
1010 
1011 END SUBROUTINE proj_polar_stereographic
1012 
1013 ELEMENTAL SUBROUTINE unproj_polar_stereographic(x,y,lon,lat, &
1014  lov, lad, projection_center_flag)
1015 DOUBLE PRECISION, INTENT(in) :: x,y
1016 DOUBLE PRECISION, INTENT(out) :: lon,lat
1017 DOUBLE PRECISION, INTENT(in) :: lov, lad
1018 INTEGER, INTENT(in) :: projection_center_flag
1019 
1020 DOUBLE PRECISION :: ro, c, pollat
1021 
1022 IF (iand(projection_center_flag, 128) == 0) THEN
1023  pollat = 90.d0*degrad
1024 ELSE
1025  pollat = -90.d0*degrad
1026 ENDIF
1027 
1028 ro = sqrt(x**2 + y**2)
1029 c = 2.0d0*atan(ro/(2.0d0*rearth))
1030 lat = raddeg*asin(cos(c)*sin(pollat)+y*sin(c)*cos(pollat)/ro)
1031 lon = lov + raddeg*atan2(x*sin(c), &
1032  (ro*cos(pollat)*cos(c)-y*sin(pollat)*sin(c)))
1033 
1034 END SUBROUTINE unproj_polar_stereographic
1035 
1036 
1037 ELEMENTAL SUBROUTINE proj_utm(lon, lat, x, y, lov, false_e, false_n, ellips)
1038 DOUBLE PRECISION,INTENT(in) :: lon,lat
1039 DOUBLE PRECISION,INTENT(out) :: x,y
1040 DOUBLE PRECISION,INTENT(in) :: lov
1041 DOUBLE PRECISION,INTENT(in) :: false_e, false_n
1042 TYPE(geo_proj_ellips),INTENT(in) :: ellips
1043 
1044 DOUBLE PRECISION :: deltalon, p
1045 DOUBLE PRECISION :: n, t, t2, c, m, a1, a2, a3, a4, a5, a6, sinp, cosp, tanp
1046 
1047 ! --- Compute delta longitude in radians
1048 deltalon = degrad*(lon - lov)
1049 
1050 ! --- Convert phi (latitude) to radians
1051 p = degrad*lat
1052 sinp = sin(p)
1053 cosp = cos(p)
1054 tanp = tan(p)
1055 
1056 n = ellips%a/sqrt(1.0d0 - ellips%e2*sinp*sinp)
1057 t = tanp*tanp
1058 c = ellips%ep2*cosp*cosp
1059 a1 = deltalon*cosp
1060 !!$m = 111132.0894_fp_utm*lat - 16216.94_fp_utm*SIN(2.0*p) + 17.21_fp_utm*SIN(4.0*p) &
1061 !!$ - 0.02_fp_utm*SIN(6.0*p)
1062 ! Modificato rispetto alla routine originale, dipende dall'ellissoide
1063 m = ellips%ef0*p - ellips%ef1*sin(2.0d0*p) + ellips%ef2*sin(4.0d0*p) - &
1064  ellips%ef3*sin(6.0d0*p)
1065 
1066 a2 = a1**2
1067 a3 = a2*a1
1068 a4 = a2**2
1069 a5 = a4*a1
1070 a6 = a4*a2
1071 t2 = t**2
1072 
1073 ! --- Compute UTM x and y (m)
1074 x = k0*n*(a1 + (1.0d0 - t + c)*a3/6.0d0 &
1075  + (5.0d0 - 18.0d0*t + t2 + 72.0d0*c - 58.0d0*ellips%ep2) &
1076  *a5/120.0d0) + false_e
1077 y = k0*(m + n*tanp &
1078  *(a2/2.0d0 + (5.0d0 - t + 9.0d0*c + &
1079  4.0d0*c*c)*a4/24.0d0 + (61.0d0 - 58.0d0*t + t2 + &
1080  600.0d0*c - 330.0d0*ellips%ep2)*a6/720.0d0))
1081 y = y + false_n
1082 
1083 END SUBROUTINE proj_utm
1084 
1085 ELEMENTAL SUBROUTINE unproj_utm(x, y, lon, lat, lov, false_e, false_n, ellips)
1086 DOUBLE PRECISION,INTENT(IN) :: x, y
1087 DOUBLE PRECISION,INTENT(OUT) :: lon, lat
1088 DOUBLE PRECISION,INTENT(in) :: lov
1089 DOUBLE PRECISION,INTENT(in) :: false_e, false_n
1090 TYPE(geo_proj_ellips),INTENT(in) :: ellips
1091 
1092 DOUBLE PRECISION :: xm, ym, m, u, p1, c1, c2, t1, t2, n1, &
1093  sinp1, cosp1, tanp1, sin2p1, r0, r1, d, d2, d3, d4, d5, d6, p, l
1094 
1095 ! --- Correct for false easting, southern hemisphere
1096 xm = x - false_e
1097 !IF (zone < 0) THEN
1098 ym = y - false_n
1099 !ELSE
1100 !ym = utmn
1101 !ENDIF
1102 
1103 m = ym/k0
1104 u = m/(ellips%a*(1.0d0-ellips%e2/4.0d0 - &
1105  3.0d0*ellips%e4/64.0d0 - 5.0d0*ellips%e6/256.0d0))
1106 p1 = u + ellips%e11*sin(2.0d0*u) + ellips%e12*sin(4.0d0*u) + &
1107  ellips%e13*sin(6.0d0*u) + ellips%e14*sin(8.0d0*u)
1108 sinp1 = sin(p1)
1109 cosp1 = cos(p1)
1110 tanp1 = tan(p1)
1111 c1 = ellips%ep2*cosp1**2
1112 c2 = c1**2
1113 t1 = tanp1**2
1114 t2 = t1**2
1115 sin2p1 = sinp1**2
1116 r0 = 1.0d0-ellips%e2*sin2p1
1117 n1 = ellips%a/sqrt(r0)
1118 !r1 = ellips%a*(1.0D0-ellips%e2)/SQRT(r0**3)
1119 r1 = r0/(1.0d0-ellips%e2)
1120 
1121 d = xm/(n1*k0)
1122 d2=d**2
1123 d3=d*d2
1124 d4=d*d3
1125 d5=d*d4
1126 d6=d*d5
1127 
1128 ! other variant from mstortini:
1129 !p = p1 - (1.0D0*tanp1/r0) * (d2/2.0D0 &
1130 ! original version with r1 = ellips%a*(1.0D0-ellips%e2)/SQRT(r0**3)
1131 !p = p1 - (n1*tanp1/r1) * (d2/2.0D0 &
1132 ! optimized version with different definition of r1
1133 p = p1 - (r1*tanp1) * (d2/2.0d0 &
1134  - (5.0d0 + 3.0d0*t1 + 10.0d0*c1 - 4.0d0*c2 &
1135  - 9.0d0*ellips%ep2)*d4/24.0d0 &
1136  + (61.0d0 + 90.0d0*t1 + 298.0d0*c1 + 45.0d0*t2 &
1137  - 252d0*ellips%ep2 - 3.0d0*c2)*d6/720.0d0)
1138 lat = raddeg*p
1139 l = (d - (1.0d0 + 2.0d0*t1 + c1)*d3/6.0d0 &
1140  + (5.0d0 - 2.0d0*c1 + 28.0d0*t1 - 3.0d0*c2 &
1141  + 8.0d0*ellips%ep2 + 24.0d0*t2)*d5/120.0d0)/cosp1
1142 lon = raddeg*l + lov
1143 
1144 END SUBROUTINE unproj_utm
1145 
1146 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:67

Generated with Doxygen.