libsim  Versione7.2.1
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, &
139 ellips_grs80 = 3, &
140 ellips_iau76 = 4, &
141 ellips_airy = 5, &
142 ellips_apl4_9 = 6, &
143 ellips_nwl9d = 7, &
144 ellips_mod_airy = 8, &
145 ellips_andrae = 9, &
146 ellips_aust_sa = 10, &
147 ellips_grs67 = 11, &
148 ellips_bessel = 12, &
149 ellips_bess_nam = 13, &
150 ellips_clrk66 = 14, &
151 ellips_clrk80 = 15, &
152 ellips_cpm = 16, &
153 ellips_delmbr = 17, &
154 ellips_engelis = 18, &
155 ellips_evrst30 = 19, &
156 ellips_evrst48 = 20, &
157 ellips_evrst56 = 21, &
158 ellips_evrst69 = 22, &
159 ellips_evrstss = 23, &
160 ellips_fschr60 = 24, &
161 ellips_fschr60m = 25, &
162 ellips_fschr68 = 26, &
163 ellips_helmert = 27, &
164 ellips_hough = 28, &
165 ellips_intl = 29, &
166 ellips_krass = 30, &
167 ellips_kaula = 31, &
168 ellips_lerch = 32, &
169 ellips_mprts = 33, &
170 ellips_new_intl = 34, &
171 ellips_plessis = 35, &
172 ellips_seasia = 36, &
173 ellips_walbeck = 37, &
174 ellips_wgs60 = 38, &
175 ellips_wgs66 = 39, &
176 ellips_wgs72 = 40, &
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(/=)
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
306 
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 
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
594 
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
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
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 
747 
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
767 
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
851 
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
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.
Compute forward coordinate transformation from geographical system to projected system.
Set of functions that return a trimmed CHARACTER representation of the input variable.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
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.
Definitions of constants and functions for working with missing values.
Costanti fisiche (DOUBLEPRECISION).
Definition: phys_const.f90:72
Destructors of the corresponding objects.
Method for testing the existence of the object.
Utilities for CHARACTER variables.
Method for setting the contents of the object.

Generated with Doxygen.