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, &
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
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)
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
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
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
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
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
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
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 
756 
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.
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:67
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.