[cig-commits] r17967 - seismo/3D/SPECFEM3D/trunk/utils/lib
liuqy at geodynamics.org
liuqy at geodynamics.org
Wed Feb 23 19:35:46 PST 2011
Author: liuqy
Date: 2011-02-23 19:35:46 -0800 (Wed, 23 Feb 2011)
New Revision: 17967
Added:
seismo/3D/SPECFEM3D/trunk/utils/lib/sph2utm.f90
seismo/3D/SPECFEM3D/trunk/utils/lib/utm2sph.f90
seismo/3D/SPECFEM3D/trunk/utils/lib/utm_geo.f90
Log:
Add sph2utm utilities to utils/lib directory.
Added: seismo/3D/SPECFEM3D/trunk/utils/lib/sph2utm.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/lib/sph2utm.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/lib/sph2utm.f90 2011-02-24 03:35:46 UTC (rev 17967)
@@ -0,0 +1,24 @@
+
+program testf
+
+implicit none
+
+include 'constants.h'
+
+double precision lon,lat
+double precision x,y
+
+integer iway, iproject
+
+iway = ILONGLAT2UTM
+iproject = 11
+
+print *, 'input lon, lat'
+read(*,*) lon, lat
+
+call utm_geo(lon,lat,x,y,iproject,iway)
+
+print *, 'x = ', x, ' y = ', y
+
+end program testf
+
Added: seismo/3D/SPECFEM3D/trunk/utils/lib/utm2sph.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/lib/utm2sph.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/lib/utm2sph.f90 2011-02-24 03:35:46 UTC (rev 17967)
@@ -0,0 +1,24 @@
+
+program testf
+
+implicit none
+
+include 'constants.h'
+
+double precision lon,lat
+double precision x,y
+
+integer iway, iproject
+
+iway = IUTM2LONGLAT
+iproject = 11
+
+print *, 'input x, y'
+read(*,*) x, y
+
+call utm_geo(lon,lat,x,y,iproject,iway)
+
+print *, 'lon = ', lon, ' lat = ', lat
+
+end program testf
+
Added: seismo/3D/SPECFEM3D/trunk/utils/lib/utm_geo.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/lib/utm_geo.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/lib/utm_geo.f90 2011-02-24 03:35:46 UTC (rev 17967)
@@ -0,0 +1,197 @@
+!=====================================================================
+!
+! UTM (Universal Transverse Mercator) projection from the USGS
+!
+!=====================================================================
+
+ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway)
+
+! convert geodetic longitude and latitude to UTM, and back
+! use iway = ILONGLAT2UTM for long/lat to UTM, IUTM2LONGLAT for UTM to lat/long
+! a list of UTM zones of the world is available at www.dmap.co.uk/utmworld.htm
+
+ implicit none
+
+ include "constants.h"
+
+!
+!-----CAMx v2.03
+!
+! UTM_GEO performs UTM to geodetic (long/lat) translation, and back.
+!
+! This is a Fortran version of the BASIC program "Transverse Mercator
+! Conversion", Copyright 1986, Norman J. Berls (Stefan Musarra, 2/94)
+! Based on algorithm taken from "Map Projections Used by the USGS"
+! by John P. Snyder, Geological Survey Bulletin 1532, USDI.
+!
+! Input/Output arguments:
+!
+! rlon Longitude (deg, negative for West)
+! rlat Latitude (deg)
+! rx UTM easting (m)
+! ry UTM northing (m)
+! UTM_PROJECTION_ZONE UTM zone
+! iway Conversion type
+! ILONGLAT2UTM = geodetic to UTM
+! IUTM2LONGLAT = UTM to geodetic
+!
+
+ integer UTM_PROJECTION_ZONE,iway
+ double precision rx,ry,rlon,rlat
+
+ double precision, parameter :: degrad=PI/180., raddeg=180./PI
+ double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0
+ double precision, parameter :: scfa=.9996d0
+ double precision, parameter :: north=0.d0, east=500000.d0
+
+ double precision e2,e4,e6,ep2,xx,yy,dlat,dlon,zone,cm,cmr,delam
+ double precision f1,f2,f3,f4,rm,rn,t,c,a,e1,u,rlat1,dlat1,c1,t1,rn1,r1,d
+ double precision rx_save,ry_save,rlon_save,rlat_save
+
+ if(SUPPRESS_UTM_PROJECTION) then
+ if (iway == ILONGLAT2UTM) then
+ rx = rlon
+ ry = rlat
+ else
+ rlon = rx
+ rlat = ry
+ endif
+ return
+ endif
+
+! save original parameters
+ rlon_save = rlon
+ rlat_save = rlat
+ rx_save = rx
+ ry_save = ry
+
+! define parameters of reference ellipsoid
+ e2=1.0-(semimin/semimaj)**2.0
+ e4=e2*e2
+ e6=e2*e4
+ ep2=e2/(1.-e2)
+
+ if (iway == IUTM2LONGLAT) then
+ xx = rx
+ yy = ry
+ else
+ dlon = rlon
+ dlat = rlat
+ endif
+!
+!----- Set Zone parameters
+!
+ zone = dble(UTM_PROJECTION_ZONE)
+ cm = zone*6.0 - 183.0
+ cmr = cm*degrad
+!
+!---- Lat/Lon to UTM conversion
+!
+ if (iway == ILONGLAT2UTM) then
+
+ rlon = degrad*dlon
+ rlat = degrad*dlat
+
+ delam = dlon - cm
+ if (delam < -180.) delam = delam + 360.
+ if (delam > 180.) delam = delam - 360.
+ delam = delam*degrad
+
+ f1 = (1. - e2/4. - 3.*e4/64. - 5.*e6/256)*rlat
+ f2 = 3.*e2/8. + 3.*e4/32. + 45.*e6/1024.
+ f2 = f2*sin(2.*rlat)
+ f3 = 15.*e4/256.*45.*e6/1024.
+ f3 = f3*sin(4.*rlat)
+ f4 = 35.*e6/3072.
+ f4 = f4*sin(6.*rlat)
+ rm = semimaj*(f1 - f2 + f3 - f4)
+ if (dlat == 90. .or. dlat == -90.) then
+ xx = 0.
+ yy = scfa*rm
+ else
+ rn = semimaj/sqrt(1. - e2*sin(rlat)**2)
+ t = tan(rlat)**2
+ c = ep2*cos(rlat)**2
+ a = cos(rlat)*delam
+
+ f1 = (1. - t + c)*a**3/6.
+ f2 = 5. - 18.*t + t**2 + 72.*c - 58.*ep2
+ f2 = f2*a**5/120.
+ xx = scfa*rn*(a + f1 + f2)
+ f1 = a**2/2.
+ f2 = 5. - t + 9.*c + 4.*c**2
+ f2 = f2*a**4/24.
+ f3 = 61. - 58.*t + t**2 + 600.*c - 330.*ep2
+ f3 = f3*a**6/720.
+ yy = scfa*(rm + rn*tan(rlat)*(f1 + f2 + f3))
+ endif
+ xx = xx + east
+ yy = yy + north
+
+!
+!---- UTM to Lat/Lon conversion
+!
+ else
+
+ xx = xx - east
+ yy = yy - north
+ e1 = sqrt(1. - e2)
+ e1 = (1. - e1)/(1. + e1)
+ rm = yy/scfa
+ u = 1. - e2/4. - 3.*e4/64. - 5.*e6/256.
+ u = rm/(semimaj*u)
+
+ f1 = 3.*e1/2. - 27.*e1**3./32.
+ f1 = f1*sin(2.*u)
+ f2 = 21.*e1**2/16. - 55.*e1**4/32.
+ f2 = f2*sin(4.*u)
+ f3 = 151.*e1**3./96.
+ f3 = f3*sin(6.*u)
+ rlat1 = u + f1 + f2 + f3
+ dlat1 = rlat1*raddeg
+ if (dlat1 >= 90. .or. dlat1 <= -90.) then
+ dlat1 = dmin1(dlat1,dble(90.) )
+ dlat1 = dmax1(dlat1,dble(-90.) )
+ dlon = cm
+ else
+ c1 = ep2*cos(rlat1)**2.
+ t1 = tan(rlat1)**2.
+ f1 = 1. - e2*sin(rlat1)**2.
+ rn1 = semimaj/sqrt(f1)
+ r1 = semimaj*(1. - e2)/sqrt(f1**3)
+ d = xx/(rn1*scfa)
+
+ f1 = rn1*tan(rlat1)/r1
+ f2 = d**2/2.
+ f3 = 5.*3.*t1 + 10.*c1 - 4.*c1**2 - 9.*ep2
+ f3 = f3*d**2*d**2/24.
+ f4 = 61. + 90.*t1 + 298.*c1 + 45.*t1**2. - 252.*ep2 - 3.*c1**2
+ f4 = f4*(d**2)**3./720.
+ rlat = rlat1 - f1*(f2 - f3 + f4)
+ dlat = rlat*raddeg
+
+ f1 = 1. + 2.*t1 + c1
+ f1 = f1*d**2*d/6.
+ f2 = 5. - 2.*c1 + 28.*t1 - 3.*c1**2 + 8.*ep2 + 24.*t1**2.
+ f2 = f2*(d**2)**2*d/120.
+ rlon = cmr + (d - f1 + f2)/cos(rlat1)
+ dlon = rlon*raddeg
+ if (dlon < -180.) dlon = dlon + 360.
+ if (dlon > 180.) dlon = dlon - 360.
+ endif
+ endif
+
+ if (iway == IUTM2LONGLAT) then
+ rlon = dlon
+ rlat = dlat
+ rx = rx_save
+ ry = ry_save
+ else
+ rx = xx
+ ry = yy
+ rlon = rlon_save
+ rlat = rlat_save
+ endif
+
+ end subroutine utm_geo
+
More information about the CIG-COMMITS
mailing list