libsim Versione 7.2.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"
19MODULE geo_proj_class
24IMPLICIT NONE
25
26TYPE 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
31END TYPE geo_proj_polar
32
33TYPE geo_proj_rotated
34 DOUBLE PRECISION :: longitude_south_pole, latitude_south_pole, angle_rotation
35END TYPE geo_proj_rotated
36
37TYPE geo_proj_stretched
38 DOUBLE PRECISION :: latitude_stretch_pole, longitude_stretch_pole, stretch_factor
39END TYPE geo_proj_stretched
40
41TYPE 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
45END TYPE geo_proj_ellips
46
47TYPE 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
56END TYPE geo_proj
57
58
60INTERFACE delete
61 MODULE PROCEDURE geo_proj_delete
62END INTERFACE
63
65INTERFACE copy
66 MODULE PROCEDURE geo_proj_copy
67END INTERFACE
68
70INTERFACE get_val
71 MODULE PROCEDURE geo_proj_get_val
72END INTERFACE
73
75INTERFACE set_val
76 MODULE PROCEDURE geo_proj_set_val
77END INTERFACE
78
80INTERFACE c_e
81 MODULE PROCEDURE geo_proj_c_e
82END INTERFACE
83
85INTERFACE write_unit
86 MODULE PROCEDURE geo_proj_write_unit
87END INTERFACE
88
90INTERFACE read_unit
91 MODULE PROCEDURE geo_proj_read_unit
92END INTERFACE
93
95INTERFACE display
96 MODULE PROCEDURE geo_proj_display
97END INTERFACE
98
101INTERFACE proj
102 MODULE PROCEDURE geo_proj_proj
103END INTERFACE
104
107INTERFACE unproj
108 MODULE PROCEDURE geo_proj_unproj
109END INTERFACE
110
114INTERFACE OPERATOR (==)
115 MODULE PROCEDURE geo_proj_eq
116END INTERFACE
117
121INTERFACE OPERATOR (/=)
122 MODULE PROCEDURE geo_proj_ne
123END INTERFACE
124
125
126INTEGER,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
131INTEGER,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 ).
136INTEGER,PARAMETER,PUBLIC :: &
137ellips_merit = 1, & !< constants for predefined ellipsoids MERIT 1983
138ellips_sgs85 = 2, &
139ellips_grs80 = 3, &
140ellips_iau76 = 4, &
141ellips_airy = 5, &
142ellips_apl4_9 = 6, &
143ellips_nwl9d = 7, &
144ellips_mod_airy = 8, &
145ellips_andrae = 9, &
146ellips_aust_sa = 10, &
147ellips_grs67 = 11, &
148ellips_bessel = 12, &
149ellips_bess_nam = 13, &
150ellips_clrk66 = 14, &
151ellips_clrk80 = 15, &
152ellips_cpm = 16, &
153ellips_delmbr = 17, &
154ellips_engelis = 18, &
155ellips_evrst30 = 19, &
156ellips_evrst48 = 20, &
157ellips_evrst56 = 21, &
158ellips_evrst69 = 22, &
159ellips_evrstss = 23, &
160ellips_fschr60 = 24, &
161ellips_fschr60m = 25, &
162ellips_fschr68 = 26, &
163ellips_helmert = 27, &
164ellips_hough = 28, &
165ellips_intl = 29, &
166ellips_krass = 30, &
167ellips_kaula = 31, &
168ellips_lerch = 32, &
169ellips_mprts = 33, &
170ellips_new_intl = 34, &
171ellips_plessis = 35, &
172ellips_seasia = 36, &
173ellips_walbeck = 37, &
174ellips_wgs60 = 38, &
175ellips_wgs66 = 39, &
176ellips_wgs72 = 40, &
177ellips_wgs84 = 41
178
179DOUBLE 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 /)
222DOUBLE 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
266DOUBLE PRECISION,PARAMETER,PRIVATE :: k0=0.9996d0 ! scale factor at central meridian (check whether this is correct and constant)
267
268PRIVATE
269PUBLIC 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
275CONTAINS
276
283FUNCTION geo_proj_new(proj_type, lov, zone, xoff, yoff, &
284 longitude_south_pole, latitude_south_pole, angle_rotation, &
285 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
286 latin1, latin2, lad, projection_center_flag, &
287 ellips_smaj_axis, ellips_flatt, ellips_type) RESULT(this)
288CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type
289DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov
290INTEGER,INTENT(in),OPTIONAL :: zone
291DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff
292DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff
293DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole
294DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole
295DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation
296DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole
297DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole
298DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor
299DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1
300DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2
301DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad
302INTEGER,INTENT(in),OPTIONAL :: projection_center_flag
303DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
304DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
305INTEGER,INTENT(in),OPTIONAL :: ellips_type
306
307TYPE(geo_proj) :: this
308
309this%proj_type = optio_c(proj_type, len(this%proj_type))
310
311! line of view / central meridian, use set_val
312this%lov = optio_d(lov)
313CALL set_val(this, zone=zone, lov=lov)
314
315! offset / false *ing
316this%xoff = optio_d(xoff)
317IF (.NOT.c_e(this%xoff)) this%xoff = 0.0d0
318this%yoff = optio_d(yoff)
319IF (.NOT.c_e(this%yoff)) this%yoff = 0.0d0
320
321this%rotated%longitude_south_pole = optio_d(longitude_south_pole)
322this%rotated%latitude_south_pole = optio_d(latitude_south_pole)
323this%rotated%angle_rotation = optio_d(angle_rotation)
324this%stretched%longitude_stretch_pole = optio_d(longitude_stretch_pole)
325this%stretched%latitude_stretch_pole = optio_d(latitude_stretch_pole)
326this%stretched%stretch_factor = optio_d(stretch_factor)
327this%polar%latin1 = optio_d(latin1)
328this%polar%latin2 = optio_d(latin2)
329this%polar%lad = optio_d(lad)
330this%polar%projection_center_flag = optio_l(projection_center_flag)
331
332! ellipsoid, start from sphere, then use set_val
333CALL ellips_compute(this%ellips)
334CALL set_val(this, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt, &
335 ellips_type=ellips_type)
336
337END 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
342SUBROUTINE ellips_compute(this, a, f)
343TYPE(geo_proj_ellips),INTENT(inout) :: this
344DOUBLE PRECISION,INTENT(in),OPTIONAL :: a, f
345
346IF (PRESENT(a) .AND. PRESENT(f)) THEN
347 this%f = f
348 this%a = a
349ELSE IF (PRESENT(a)) THEN ! parameters for a spherical Earth with given radius
350 this%f = 0.0d0
351 this%a = a
352ELSE ! parameters for a standard spherical Earth
353 this%f = 0.0d0
354 this%a = rearth
355ENDIF
356
357this%e2 = 2.0d0*this%f - this%f*this%f ! Eccentricity
358this%e1 = this%f/(2.0d0 - this%f)
359this%ep2 = this%e2/(1.0d0 - this%e2)
360this%e11 = 3.0d0*this%e1/2.0d0 - 27.0d0*this%e1*this%e1*this%e1/32.0d0
361this%e12 = 21.0d0*this%e1*this%e1/16.0d0 &
362 - 55.0d0*this%e1*this%e1*this%e1*this%e1/32.0d0
363this%e13 = 151.0d0*this%e1*this%e1*this%e1/96.0d0
364this%e14 = 1097.0d0*this%e1*this%e1*this%e1*this%e1/512.0d0
365this%e4 = this%e2*this%e2
366this%e6 = this%e2*this%e4
367this%ef0 = this%a*(1.0d0 - 0.25d0*this%e2&
368 *(1.0d0+this%e2/16.0d0*(3.0d0 + 1.25d0*this%e2)))
369this%ef1 = this%a*(0.375d0*this%e2 &
370 *(1.0d0 + 0.25d0*this%e2*(1.0d0 + 0.46875d0*this%e2)))
371this%ef2 = this%a*(0.05859375d0*this%e2*this%e2*(1.0d0 + 0.75d0*this%e2))
372this%ef3 = this%a*this%e2*this%e2*this%e2*35.0d0/3072.0d0
373
374END SUBROUTINE ellips_compute
375
376
378SUBROUTINE geo_proj_delete(this)
379TYPE(geo_proj),INTENT(inout) :: this
380
381this%proj_type = cmiss
382this%lov = dmiss
383this%xoff = dmiss
384this%yoff = dmiss
385this%rotated%longitude_south_pole = dmiss
386this%rotated%latitude_south_pole = dmiss
387this%rotated%angle_rotation = dmiss
388this%stretched%longitude_stretch_pole = dmiss
389this%stretched%latitude_stretch_pole = dmiss
390this%stretched%stretch_factor = dmiss
391this%polar%latin1 = dmiss
392this%polar%latin2 = dmiss
393this%polar%lad = dmiss
394this%polar%projection_center_flag = imiss
395
396END SUBROUTINE geo_proj_delete
397
398
400SUBROUTINE geo_proj_copy(this, that)
401TYPE(geo_proj),INTENT(in) :: this
402TYPE(geo_proj),INTENT(out) :: that
403
404that = this
405
406END SUBROUTINE geo_proj_copy
407
408
409SUBROUTINE geo_proj_get_val(this, &
410 proj_type, lov, zone, xoff, yoff, &
411 longitude_south_pole, latitude_south_pole, angle_rotation, &
412 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
413 latin1, latin2, lad, projection_center_flag, &
414 ellips_smaj_axis, ellips_flatt, ellips_type, unit)
415TYPE(geo_proj),INTENT(in) :: this
416CHARACTER(len=*),OPTIONAL :: proj_type
417DOUBLE PRECISION,OPTIONAL :: lov
418INTEGER,OPTIONAL :: zone
419DOUBLE PRECISION,OPTIONAL :: xoff
420DOUBLE PRECISION,OPTIONAL :: yoff
421DOUBLE PRECISION,OPTIONAL :: longitude_south_pole
422DOUBLE PRECISION,OPTIONAL :: latitude_south_pole
423DOUBLE PRECISION,OPTIONAL :: angle_rotation
424DOUBLE PRECISION,OPTIONAL :: longitude_stretch_pole
425DOUBLE PRECISION,OPTIONAL :: latitude_stretch_pole
426DOUBLE PRECISION,OPTIONAL :: stretch_factor
427DOUBLE PRECISION,OPTIONAL :: latin1
428DOUBLE PRECISION,OPTIONAL :: latin2
429DOUBLE PRECISION,OPTIONAL :: lad
430INTEGER,OPTIONAL :: projection_center_flag
431DOUBLE PRECISION,OPTIONAL :: ellips_smaj_axis
432DOUBLE PRECISION,OPTIONAL :: ellips_flatt
433INTEGER,OPTIONAL :: ellips_type
434INTEGER,OPTIONAL :: unit
435
436INTEGER :: i
437
438IF (PRESENT(proj_type)) proj_type = this%proj_type
439
440IF (PRESENT(lov) .AND. PRESENT(zone)) THEN
441 zone = nint((this%lov + 183.0d0)/6.0d0)
442 lov = this%lov - zone*6.0d0 - 183.0d0
443 IF (abs(lov) < 1.0d-6) lov = 0.0d-6
444ELSE IF (PRESENT(lov)) THEN
445 lov = this%lov
446ELSE IF (PRESENT(zone)) THEN
447 zone = nint((this%lov + 183.0d0)/6.0d0)
448ENDIF
449
450IF (PRESENT(xoff)) xoff = this%xoff
451IF (PRESENT(yoff)) yoff = this%yoff
452IF (PRESENT(longitude_south_pole)) longitude_south_pole = this%rotated%longitude_south_pole
453IF (PRESENT(latitude_south_pole)) latitude_south_pole = this%rotated%latitude_south_pole
454IF (PRESENT(angle_rotation)) angle_rotation = this%rotated%angle_rotation
455IF (PRESENT(longitude_stretch_pole)) longitude_stretch_pole = this%stretched%longitude_stretch_pole
456IF (PRESENT(latitude_stretch_pole)) latitude_stretch_pole = this%stretched%latitude_stretch_pole
457IF (PRESENT(stretch_factor)) stretch_factor = this%stretched%stretch_factor
458IF (PRESENT(latin1)) latin1 = this%polar%latin1
459IF (PRESENT(latin2)) latin2 = this%polar%latin2
460IF (PRESENT(lad)) lad = this%polar%lad
461IF (PRESENT(projection_center_flag)) projection_center_flag = this%polar%projection_center_flag
462! ellipsoid
463IF (PRESENT(ellips_smaj_axis)) ellips_smaj_axis = this%ellips%a
464IF (PRESENT(ellips_flatt)) ellips_flatt = this%ellips%f
465IF (PRESENT(ellips_type)) THEN
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
473ENDIF
474
475IF (PRESENT(unit)) THEN
476 SELECT CASE(this%proj_type)
477 CASE("regular_ll", "rotated_ll")
478 unit = geo_proj_unit_degree
479 CASE("lambert", "polar_stereographic", "UTM")
480 unit = geo_proj_unit_meter
481 CASE default
482 unit = imiss
483 END SELECT
484ENDIF
485
486
487END SUBROUTINE geo_proj_get_val
488
489
490SUBROUTINE geo_proj_set_val(this, &
491 proj_type, lov, zone, xoff, yoff, &
492 longitude_south_pole, latitude_south_pole, angle_rotation, &
493 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
494 latin1, latin2, lad, projection_center_flag, &
495 ellips_smaj_axis, ellips_flatt, ellips_type)
496TYPE(geo_proj),INTENT(inout) :: this
497CHARACTER(len=*),OPTIONAL :: proj_type
498DOUBLE PRECISION,OPTIONAL :: lov
499INTEGER,INTENT(in),OPTIONAL :: zone
500DOUBLE PRECISION,OPTIONAL :: xoff
501DOUBLE PRECISION,OPTIONAL :: yoff
502DOUBLE PRECISION,OPTIONAL :: longitude_south_pole
503DOUBLE PRECISION,OPTIONAL :: latitude_south_pole
504DOUBLE PRECISION,OPTIONAL :: angle_rotation
505DOUBLE PRECISION,OPTIONAL :: longitude_stretch_pole
506DOUBLE PRECISION,OPTIONAL :: latitude_stretch_pole
507DOUBLE PRECISION,OPTIONAL :: stretch_factor
508DOUBLE PRECISION,OPTIONAL :: latin1
509DOUBLE PRECISION,OPTIONAL :: latin2
510DOUBLE PRECISION,OPTIONAL :: lad
511INTEGER,OPTIONAL :: projection_center_flag
512DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis
513DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt
514INTEGER,INTENT(in),OPTIONAL :: ellips_type
515
516INTEGER :: lzone
517DOUBLE PRECISION :: llov
518
519
520! line of view / central meridian
521llov = optio_d(lov)
522lzone = optio_i(zone)
523IF (c_e(llov) .AND. c_e(lzone)) THEN
524 this%lov = llov + zone*6.0d0 - 183.0d0
525ELSE IF (c_e(llov)) THEN
526 this%lov = llov
527ELSE IF (c_e(lzone)) THEN
528 this%lov = lzone*6.0d0 - 183.0d0
529ENDIF
530
531! ellipsoid
532IF (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)
535ELSE 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
542ENDIF
543
544! todo all the rest
545
546END SUBROUTINE geo_proj_set_val
547
548
549FUNCTION geo_proj_c_e(this) RESULT(c_e)
550TYPE(geo_proj),INTENT(in) :: this
551LOGICAL :: c_e
552
553c_e = this%proj_type /= cmiss
554
555END FUNCTION geo_proj_c_e
556
557
562SUBROUTINE geo_proj_read_unit(this, unit)
563TYPE(geo_proj),INTENT(out) :: this
564INTEGER, INTENT(in) :: unit
565
566CHARACTER(len=40) :: form
567
568INQUIRE(unit, form=form)
569IF (form == 'FORMATTED') THEN
570 READ(unit,*)this
571ELSE
572 READ(unit)this
573ENDIF
574
575END SUBROUTINE geo_proj_read_unit
576
577
582SUBROUTINE geo_proj_write_unit(this, unit)
583TYPE(geo_proj),INTENT(in) :: this
584INTEGER, INTENT(in) :: unit
585
586CHARACTER(len=40) :: form
587
588INQUIRE(unit, form=form)
589IF (form == 'FORMATTED') THEN
590 WRITE(unit,*)this
591ELSE
592 WRITE(unit)this
593ENDIF
594
595END SUBROUTINE geo_proj_write_unit
596
598SUBROUTINE geo_proj_display(this)
599TYPE(geo_proj),INTENT(in) :: this
600
601print*,"<<<<<<<<<<<<<<< ",t2c(this%proj_type,"undefined projection"), &
602 " >>>>>>>>>>>>>>>>"
603
604IF (c_e(this%xoff) .AND. this%xoff /= 0.0d0) THEN
605 print*,"False easting",this%xoff
606ENDIF
607IF (c_e(this%yoff) .AND. this%yoff /= 0.0d0) THEN
608 print*,"False northing",this%yoff
609ENDIF
610
611IF (c_e(this%lov)) THEN
612 print*,"Central Meridian",this%lov
613ENDIF
614
615IF (this%proj_type == 'rotated_ll' .OR. this%proj_type == 'stretched_rotated_ll') THEN
616 print*,"Rotated projection:"
617 print*,"lon of south pole",this%rotated%longitude_south_pole
618 print*,"lat of south pole",this%rotated%latitude_south_pole
619 print*,"angle of rotation",this%rotated%angle_rotation
620ENDIF
621
622IF (this%proj_type == 'stretched_ll' .OR. this%proj_type == 'stretched_rotated_ll') THEN
623 print*,"Stretched projection:"
624 print*,"lon of stretch pole",this%stretched%longitude_stretch_pole
625 print*,"lat of stretch pole",this%stretched%latitude_stretch_pole
626 print*,"stretching factor",this%stretched%stretch_factor
627ENDIF
628
629IF (this%proj_type == 'lambert' .OR. this%proj_type == 'polar_stereographic') THEN
630 print*,"Polar projection:"
631 IF (c_e(this%polar%latin1) .OR. c_e(this%polar%latin2)) THEN
632 print*,"lat of intersections",this%polar%latin1,this%polar%latin2
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
642ENDIF
643
644IF (this%proj_type == 'mercator') THEN
645 IF (c_e(this%polar%lad)) THEN
646 print*,"isometric latitude",this%polar%lad
647 ENDIF
648ENDIF
649
650IF (this%ellips%f == 0.0d0) THEN
651 print*,"Spherical Earth:"
652 print*,"Radius (m)",this%ellips%a
653ELSE
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
658ENDIF
659
660
661END SUBROUTINE geo_proj_display
662
663
666ELEMENTAL SUBROUTINE geo_proj_proj(this, lon, lat, x, y)
667TYPE(geo_proj),INTENT(in) :: this
669DOUBLE PRECISION, INTENT(in) :: lon, lat
671DOUBLE PRECISION, INTENT(out) :: x, y
672
673SELECT CASE(this%proj_type)
674
675CASE("regular_ll")
676 CALL proj_regular_ll(lon, lat, x, y)
677
678CASE("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
682CASE("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
687CASE("polar_stereographic")
688 CALL proj_polar_stereographic(lon, lat, x, y, this%lov, &
689 this%polar%lad, this%polar%projection_center_flag)
690
691CASE("mercator")
692 CALL proj_mercator(lon, lat, x, y, this%lov, this%polar%lad)
693
694CASE("UTM")
695 CALL proj_utm(lon, lat, x, y, this%lov, this%xoff, this%yoff, this%ellips)
696
697CASE default
698 x = dmiss
699 y = dmiss
700
701END SELECT
702
703END SUBROUTINE geo_proj_proj
704
705
708ELEMENTAL SUBROUTINE geo_proj_unproj(this, x, y, lon, lat)
709TYPE(geo_proj),INTENT(in) :: this
711DOUBLE PRECISION, INTENT(in) :: x, y
713DOUBLE PRECISION, INTENT(out) :: lon, lat
714
715SELECT CASE(this%proj_type)
716
717CASE("regular_ll")
718 CALL unproj_regular_ll(x, y, lon, lat)
719
720CASE("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
724CASE("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
729CASE("polar_stereographic")
730 CALL unproj_polar_stereographic(x, y, lon, lat, this%lov, &
731 this%polar%lad, this%polar%projection_center_flag)
732
733CASE("mercator")
734 CALL unproj_mercator(x, y, lon, lat, this%lov, this%polar%lad)
735
736CASE("UTM")
737 CALL unproj_utm(x, y, lon, lat, this%lov, this%xoff, this%yoff, this%ellips)
738
739CASE default
740 lon = dmiss
741 lat = dmiss
742
743END SELECT
744
745END SUBROUTINE geo_proj_unproj
746
747
748ELEMENTAL FUNCTION geo_proj_eq(this, that) RESULT(eq)
749TYPE(geo_proj),INTENT(in) :: this, that
750LOGICAL :: eq
751
752eq = this%proj_type == that%proj_type .AND. this%xoff == that%xoff .AND. &
753 this%yoff == that%yoff .AND. geo_lon_eq(this%lov, that%lov) .AND. &
754 geo_lon_eq(this%rotated%longitude_south_pole, that%rotated%longitude_south_pole) .AND. &
755 this%rotated%latitude_south_pole == that%rotated%latitude_south_pole .AND. &
756 this%rotated%angle_rotation == that%rotated%angle_rotation .AND. &
757 this%stretched%latitude_stretch_pole == that%stretched%latitude_stretch_pole .AND. &
758 geo_lon_eq(this%stretched%longitude_stretch_pole, that%stretched%longitude_stretch_pole) .AND. &
759 this%stretched%stretch_factor == that%stretched%stretch_factor .AND. &
760 this%polar%latin1 == that%polar%latin1 .AND. & ! 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
766END FUNCTION geo_proj_eq
767
768
769ELEMENTAL FUNCTION geo_proj_ne(this, that) RESULT(ne)
770TYPE(geo_proj),INTENT(in) :: this, that
771LOGICAL :: ne
772
773ne = .NOT. (this == that)
774
775END FUNCTION geo_proj_ne
776
777
780ELEMENTAL FUNCTION geo_lon_eq(l1, l2) RESULT(eq)
781DOUBLE PRECISION,INTENT(in) :: l1
782DOUBLE PRECISION,INTENT(in) :: l2
783
784LOGICAL :: eq
785
786!eq = (l1 == l2) .OR. (ABS(l2-l1) == 360.0D0)
787eq = modulo(l2-l1, 360.0d0) == 0.0d0
788
789END FUNCTION geo_lon_eq
790
791! =====================
792! == transformations ==
793! =====================
794
795ELEMENTAL SUBROUTINE proj_regular_ll(lon,lat,x,y)
796DOUBLE PRECISION, INTENT(in) :: lon,lat
797DOUBLE PRECISION, INTENT(out) :: x,y
798
799x = lon
800y = lat
801
802END SUBROUTINE proj_regular_ll
803
804ELEMENTAL SUBROUTINE unproj_regular_ll(x,y,lon,lat)
805DOUBLE PRECISION, INTENT(in) :: x,y
806DOUBLE PRECISION, INTENT(out) :: lon,lat
807
808lon = x
809lat = y
810
811END SUBROUTINE unproj_regular_ll
812
813
814ELEMENTAL SUBROUTINE proj_rotated_ll(lon,lat,x,y, &
815 longitude_south_pole, latitude_south_pole, angle_rotation)
816DOUBLE PRECISION, INTENT(in) :: lon,lat
817DOUBLE PRECISION, INTENT(out) :: x,y
818DOUBLE PRECISION, INTENT(in) :: longitude_south_pole, latitude_south_pole, &
819 angle_rotation
820
821DOUBLE 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))
825l_south_pole = (latitude_south_pole+90.)*degrad
826
827rx = degrad*(lon - longitude_south_pole)
828srx = sin(rx)
829crx = cos(rx)
830
831sy0 = sin(l_south_pole)
832cy0 = cos(l_south_pole)
833
834sy = sin(degrad*lat)
835cy = cos(degrad*lat)
836
837x = raddeg*atan2(cy*srx, cy0*cy*crx+sy0*sy)
838y = raddeg*asin(cy0*sy - sy0*cy*crx)
839x = x + angle_rotation ! check
840
841END SUBROUTINE proj_rotated_ll
842
843ELEMENTAL SUBROUTINE unproj_rotated_ll(x,y,lon,lat,&
844 longitude_south_pole, latitude_south_pole, angle_rotation)
845DOUBLE PRECISION, INTENT(in) :: x,y
846DOUBLE PRECISION, INTENT(out) :: lon,lat
847DOUBLE PRECISION, INTENT(in) :: longitude_south_pole, latitude_south_pole, &
848 angle_rotation
849
850DOUBLE PRECISION :: cy0, sy0, l_south_pole, xr
851
852xr = (x - angle_rotation)*degrad ! check
853! old hibu formula
854!l_south_pole = ACOS(-SIN(degrad*latitude_south_pole))
855l_south_pole = (latitude_south_pole+90.)*degrad
856
857cy0 = cos(l_south_pole)
858sy0 = sin(l_south_pole)
859
860lat = raddeg*asin(sy0*cos(degrad*y)*cos(xr)+cy0*sin(degrad*y))
861lon = longitude_south_pole + &
862 raddeg*asin(sin(xr)*cos(degrad*y)/cos(degrad*lat))
863
864END SUBROUTINE unproj_rotated_ll
865
866! come usare il polo? ruotare e antiruotare?
867ELEMENTAL SUBROUTINE proj_stretched_ll(lon,lat,x,y, &
868 longitude_stretch_pole, latitude_stretch_pole, stretch_factor)
869DOUBLE PRECISION, INTENT(in) :: lon,lat
870DOUBLE PRECISION, INTENT(out) :: x,y
871DOUBLE PRECISION, INTENT(in) :: longitude_stretch_pole, latitude_stretch_pole, &
872 stretch_factor
873
874DOUBLE PRECISION :: csq
875
876csq = stretch_factor**2
877x = lon
878y = raddeg*asin((1.0d0 - csq + (1.0d0 + csq)*sin(degrad*lat)) / &
879 (1.0d0 + csq + (1.0d0 - csq)*sin(degrad*lat)))
880
881END SUBROUTINE proj_stretched_ll
882
883ELEMENTAL SUBROUTINE unproj_stretched_ll(x,y,lon,lat,&
884 longitude_stretch_pole, latitude_stretch_pole, stretch_factor)
885DOUBLE PRECISION, INTENT(in) :: x,y
886DOUBLE PRECISION, INTENT(out) :: lon,lat
887DOUBLE PRECISION, INTENT(in) :: longitude_stretch_pole, latitude_stretch_pole, &
888 stretch_factor
889
890DOUBLE PRECISION :: csq
891
892csq = stretch_factor**2
893lon = x
894! TODO verificare la formula inversa
895lat = raddeg*asin((csq - 1.0d0 + (csq + 1.0d0)*sin(degrad*y)) / &
896 (csq + 1.0d0 + (csq - 1.0d0)*sin(degrad*y)))
897
898END 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.
909ELEMENTAL SUBROUTINE proj_lambert(lon,lat,x,y, &
910 latin1, latin2, lov, lad, projection_center_flag)
911DOUBLE PRECISION, INTENT(in) :: lon,lat
912DOUBLE PRECISION, INTENT(out) :: x,y
913DOUBLE PRECISION, INTENT(in) :: latin1, latin2, lov, lad
914INTEGER, INTENT(in) :: projection_center_flag
915
916DOUBLE PRECISION :: n, f, ro0, ro, cs1, cs2, cs3, pollat, angle, cot
917DOUBLE PRECISION, PARAMETER :: epsy = 1.0d-100
918
919IF (iand(projection_center_flag, 128) == 0) THEN
920 pollat = 90.d0*degrad
921ELSE
922 pollat = -90.d0*degrad
923ENDIF
924cs1 = cos(degrad*latin1)
925cs2 = tan(pi*.25d0 + degrad*latin1*.5d0)
926
927IF (latin1 == latin2) THEN
928 n = sin(degrad*latin1) ! verify that n->sin(latin1) when latin2->latin1
929ELSE
930 n = log(cs1/cos(degrad*latin2)) / &
931 log(tan(pi*.25d0 + degrad*latin2*.5d0) / cs2)
932ENDIF
933f = cs1*cs2**n/n*rearth ! check that rearth is correct here (only if lad==latin1?)
934angle = pi*.25d0 + pollat*.5d0
935cot = cos(angle)/sin(angle)
936IF (cot > epsy) THEN
937 ro0 = f*cot**n
938ELSE
939 ro0 = 0.0d0
940ENDIF
941
942angle = pi*.25d0 + degrad*lat*.5d0
943cot = cos(angle)/sin(angle)
944IF (cot > epsy) THEN
945 ro = f*cot**n
946ELSE
947 ro = 0.0d0
948ENDIF
949
950cs3 = degrad*n*(lon - lov)
951
952x = ro*sin(cs3)
953y = ro0 - ro*cos(cs3)
954
955END SUBROUTINE proj_lambert
956
957ELEMENTAL SUBROUTINE unproj_lambert(x,y,lon,lat, &
958 latin1, latin2, lov, lad, projection_center_flag)
959DOUBLE PRECISION, INTENT(in) :: x,y
960DOUBLE PRECISION, INTENT(out) :: lon,lat
961DOUBLE PRECISION, INTENT(in) :: latin1, latin2, lov, lad
962INTEGER, INTENT(in) :: projection_center_flag
963
964DOUBLE PRECISION :: n, f, ro0, ro, theta, cs1, cs2, pollat, angle, cot
965DOUBLE 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
969IF (iand(projection_center_flag, 128) == 0) THEN
970 pollat = 90.d0*degrad
971ELSE
972 pollat = -90.d0*degrad
973ENDIF
974cs1 = cos(degrad*latin1)
975cs2 = tan(pi*.25d0 + degrad*latin1*.5d0)
976
977IF (latin1 == latin2) THEN
978 n = sin(degrad*latin1) ! verify limit
979ELSE
980 n = log(cs1/cos(degrad*latin2)) / &
981 log(tan(pi*.25d0 + degrad*latin2*.5d0) / cs2)
982ENDIF
983f = cs1*cs2**n/n*rearth ! check that rearth is correct here (only if lad==latin1?)
984angle = pi*.25d0 + pollat*.5d0
985cot = cos(angle)/sin(angle)
986IF (cot > epsy) THEN
987 ro0 = f*cot**n
988ELSE
989 ro0 = 0.0d0
990ENDIF
991
992ro = sign(sqrt(x*x + (ro0-y)*(ro0-y)), n) ! check SIGN
993theta = raddeg*atan2(x, ro0-y)
994
995lon = lov + theta/n
996lat = raddeg*(2.d0*atan((f/ro)**(1.d0/n)) - pi*.5d0)
997
998END SUBROUTINE unproj_lambert
999
1000
1001!http://mathworld.wolfram.com/StereographicProjection.html
1002ELEMENTAL SUBROUTINE proj_polar_stereographic(lon,lat,x,y, &
1003 lov, lad, projection_center_flag)
1004DOUBLE PRECISION, INTENT(in) :: lon,lat
1005DOUBLE PRECISION, INTENT(out) :: x,y
1006DOUBLE PRECISION, INTENT(in) :: lov, lad
1007INTEGER, INTENT(in) :: projection_center_flag
1008
1009DOUBLE PRECISION :: k, pollat
1010
1011IF (iand(projection_center_flag, 128) == 0) THEN
1012 pollat = 90.d0*degrad
1013ELSE
1014 pollat = -90.d0*degrad
1015ENDIF
1016
1017k = 2.0d0*rearth/(1.0d0 + sin(pollat)*sin(degrad*lat) + &
1018 cos(pollat)*cos(degrad*lat)*cos(degrad*(lon - lov)))
1019x = k*cos(degrad*lat)*sin(degrad*(lon - lov))
1020y = k*(cos(pollat)*sin(degrad*lat) - &
1021 sin(pollat)*cos(degrad*lat)*cos(degrad*(lon - lov)))
1022
1023END SUBROUTINE proj_polar_stereographic
1024
1025ELEMENTAL SUBROUTINE unproj_polar_stereographic(x,y,lon,lat, &
1026 lov, lad, projection_center_flag)
1027DOUBLE PRECISION, INTENT(in) :: x,y
1028DOUBLE PRECISION, INTENT(out) :: lon,lat
1029DOUBLE PRECISION, INTENT(in) :: lov, lad
1030INTEGER, INTENT(in) :: projection_center_flag
1031
1032DOUBLE PRECISION :: ro, c, pollat
1033
1034IF (iand(projection_center_flag, 128) == 0) THEN
1035 pollat = 90.d0*degrad
1036ELSE
1037 pollat = -90.d0*degrad
1038ENDIF
1039
1040ro = sqrt(x**2 + y**2)
1041c = 2.0d0*atan(ro/(2.0d0*rearth))
1042lat = raddeg*asin(cos(c)*sin(pollat)+y*sin(c)*cos(pollat)/ro)
1043lon = lov + raddeg*atan2(x*sin(c), &
1044 (ro*cos(pollat)*cos(c)-y*sin(pollat)*sin(c)))
1045
1046END SUBROUTINE unproj_polar_stereographic
1047
1048
1049!http://mathworld.wolfram.com/StereographicProjection.html
1050ELEMENTAL SUBROUTINE proj_mercator(lon, lat, x, y, lov, lad)
1051DOUBLE PRECISION,INTENT(in) :: lon,lat
1052DOUBLE PRECISION,INTENT(out) :: x,y
1053DOUBLE PRECISION,INTENT(in) :: lov
1054DOUBLE PRECISION,INTENT(in) :: lad
1055
1056DOUBLE PRECISION :: scx, scy
1057
1058scy = cos(degrad*lad)*rearth
1059scx = scy*degrad
1060x = (lon-lov)*scx
1061y = log(tan(0.25d0*pi + 0.5d0*lat*degrad))*scy
1062
1063END SUBROUTINE proj_mercator
1064
1065ELEMENTAL SUBROUTINE unproj_mercator(x, y, lon, lat, lov, lad)
1066DOUBLE PRECISION,INTENT(in) :: x,y
1067DOUBLE PRECISION,INTENT(out) :: lon,lat
1068DOUBLE PRECISION,INTENT(in) :: lov
1069DOUBLE PRECISION,INTENT(in) :: lad
1070
1071DOUBLE PRECISION :: scx, scy
1072
1073scy = cos(degrad*lad)*rearth
1074scx = scy*degrad
1075lon = x/scx + lov
1076lat = 2.0d0*atan(exp(y/scy))*raddeg-90.0d0
1077
1078END SUBROUTINE unproj_mercator
1079
1080
1081ELEMENTAL SUBROUTINE proj_utm(lon, lat, x, y, lov, false_e, false_n, ellips)
1082DOUBLE PRECISION,INTENT(in) :: lon,lat
1083DOUBLE PRECISION,INTENT(out) :: x,y
1084DOUBLE PRECISION,INTENT(in) :: lov
1085DOUBLE PRECISION,INTENT(in) :: false_e, false_n
1086TYPE(geo_proj_ellips),INTENT(in) :: ellips
1087
1088DOUBLE PRECISION :: deltalon, p
1089DOUBLE PRECISION :: n, t, t2, c, m, a1, a2, a3, a4, a5, a6, sinp, cosp, tanp
1090
1091! --- Compute delta longitude in radians
1092deltalon = degrad*(lon - lov)
1093
1094! --- Convert phi (latitude) to radians
1095p = degrad*lat
1096sinp = sin(p)
1097cosp = cos(p)
1098tanp = tan(p)
1099
1100n = ellips%a/sqrt(1.0d0 - ellips%e2*sinp*sinp)
1101t = tanp*tanp
1102c = ellips%ep2*cosp*cosp
1103a1 = 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
1107m = ellips%ef0*p - ellips%ef1*sin(2.0d0*p) + ellips%ef2*sin(4.0d0*p) - &
1108 ellips%ef3*sin(6.0d0*p)
1109
1110a2 = a1**2
1111a3 = a2*a1
1112a4 = a2**2
1113a5 = a4*a1
1114a6 = a4*a2
1115t2 = t**2
1116
1117! --- Compute UTM x and y (m)
1118x = k0*n*(a1 + (1.0d0 - t + c)*a3/6.0d0 &
1119 + (5.0d0 - 18.0d0*t + t2 + 72.0d0*c - 58.0d0*ellips%ep2) &
1120 *a5/120.0d0) + false_e
1121y = 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))
1125y = y + false_n
1126
1127END SUBROUTINE proj_utm
1128
1129ELEMENTAL SUBROUTINE unproj_utm(x, y, lon, lat, lov, false_e, false_n, ellips)
1130DOUBLE PRECISION,INTENT(IN) :: x, y
1131DOUBLE PRECISION,INTENT(OUT) :: lon, lat
1132DOUBLE PRECISION,INTENT(in) :: lov
1133DOUBLE PRECISION,INTENT(in) :: false_e, false_n
1134TYPE(geo_proj_ellips),INTENT(in) :: ellips
1135
1136DOUBLE PRECISION :: xm, ym, m, u, p1, c1, c2, t1, t2, n1, &
1137 sinp1, cosp1, tanp1, sin2p1, r0, r1, d, d2, d3, d4, d5, d6, p, l
1138
1139! --- Correct for false easting, southern hemisphere
1140xm = x - false_e
1141!IF (zone < 0) THEN
1142ym = y - false_n
1143!ELSE
1144!ym = utmn
1145!ENDIF
1146
1147m = ym/k0
1148u = m/(ellips%a*(1.0d0-ellips%e2/4.0d0 - &
1149 3.0d0*ellips%e4/64.0d0 - 5.0d0*ellips%e6/256.0d0))
1150p1 = u + ellips%e11*sin(2.0d0*u) + ellips%e12*sin(4.0d0*u) + &
1151 ellips%e13*sin(6.0d0*u) + ellips%e14*sin(8.0d0*u)
1152sinp1 = sin(p1)
1153cosp1 = cos(p1)
1154tanp1 = tan(p1)
1155c1 = ellips%ep2*cosp1**2
1156c2 = c1**2
1157t1 = tanp1**2
1158t2 = t1**2
1159sin2p1 = sinp1**2
1160r0 = 1.0d0-ellips%e2*sin2p1
1161n1 = ellips%a/sqrt(r0)
1162!r1 = ellips%a*(1.0D0-ellips%e2)/SQRT(r0**3)
1163r1 = r0/(1.0d0-ellips%e2)
1164
1165d = xm/(n1*k0)
1166d2=d**2
1167d3=d*d2
1168d4=d*d3
1169d5=d*d4
1170d6=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
1177p = p1 - (r1*tanp1) * (d2/2.0d0 &
1178 - (5.0d0 + 3.0d0*t1 + 10.0d0*c1 - 4.0d0*c2 &
1179 - 9.0d0*ellips%ep2)*d4/24.0d0 &
1180 + (61.0d0 + 90.0d0*t1 + 298.0d0*c1 + 45.0d0*t2 &
1181 - 252d0*ellips%ep2 - 3.0d0*c2)*d6/720.0d0)
1182lat = raddeg*p
1183l = (d - (1.0d0 + 2.0d0*t1 + c1)*d3/6.0d0 &
1184 + (5.0d0 - 2.0d0*c1 + 28.0d0*t1 - 3.0d0*c2 &
1185 + 8.0d0*ellips%ep2 + 24.0d0*t2)*d5/120.0d0)/cosp1
1186lon = raddeg*l + lov
1187
1188END SUBROUTINE unproj_utm
1189
1190END MODULE geo_proj_class
Set of functions that return a trimmed CHARACTER representation of the input variable.
Method for testing the existence of the object.
Copy an object, creating a fully new instance.
Destructors of the corresponding objects.
Print a brief description on stdout.
Method for returning the contents of the object.
Compute forward coordinate transformation from geographical system to projected system.
Read the object from a formatted or unformatted file.
Method for setting the contents of the object.
Compute backward coordinate transformation from projected system to geographical system.
Write the object on a formatted or unformatted file.
Utilities for CHARACTER variables.
Costanti fisiche (DOUBLEPRECISION).
Definition: phys_const.f90:72
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.

Generated with Doxygen.