From b994d558086f0ef2e70f8749654f85bf27a512b3 Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Fri, 30 Apr 2021 22:09:02 +0300 Subject: [PATCH] fixes utm_geo errors in formula for utm/geo conversions (a 20-year old mistake, finally fixed...) --- EXAMPLES/Mount_StHelens/convert_lonlat2utm.pl | 4 +- .../CMT3D/cmt3d/utm_geo.f90 | 171 +++++++++--------- src/shared/utm_geo.f90 | 129 +++++++++---- .../compute_coupling_acoustic_el.f90 | 14 +- .../compute_coupling_viscoelastic_ac.f90 | 14 +- .../codes_18_06Aug2006/wave2d_sub2.f90 | 8 +- .../SEM2D_iterate/src/wave2d_sub2.f90 | 8 +- .../cluster/model_slice/utm_geo.f90 | 8 +- utils/lib/utm_geo.f90 | 7 +- .../create_movie_GMT/utm_geo.f90 | 8 +- 10 files changed, 222 insertions(+), 149 deletions(-) diff --git a/EXAMPLES/Mount_StHelens/convert_lonlat2utm.pl b/EXAMPLES/Mount_StHelens/convert_lonlat2utm.pl index 45ae5790d..655070b48 100755 --- a/EXAMPLES/Mount_StHelens/convert_lonlat2utm.pl +++ b/EXAMPLES/Mount_StHelens/convert_lonlat2utm.pl @@ -179,7 +179,7 @@ sub geo2utm { $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 = 15.*$e4/256. + 45.*$e6/1024.; $f3 = $f3*sin(4.*$rlat); $f4 = 35.*$e6/3072.; $f4 = $f4*sin(6.*$rlat); @@ -244,7 +244,7 @@ sub geo2utm { $f1 = $rn1*tan($rlat1)/$r1; $f2 = $d**2/2.; - $f3 = 5.*3.*$t1 + 10.*$c1 - 4.*$c1**2 - 9.*$ep2; + $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.; diff --git a/src/inverse_problem_for_source/CMT3D/cmt3d/utm_geo.f90 b/src/inverse_problem_for_source/CMT3D/cmt3d/utm_geo.f90 index f805fed09..a97c3c22f 100644 --- a/src/inverse_problem_for_source/CMT3D/cmt3d/utm_geo.f90 +++ b/src/inverse_problem_for_source/CMT3D/cmt3d/utm_geo.f90 @@ -41,7 +41,8 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway,SUPPRESS_UTM_PROJECT logical SUPPRESS_UTM_PROJECTION double precision, parameter :: degrad=PI/180.d0, raddeg=180.d0/PI - double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0 + !double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0 ! Clarke 1866 ellipsoid + double precision, parameter :: semimaj=6378137.0d0, semimin=6356752.314245d0 ! WGS-84 ellipsoid double precision, parameter :: scfa=0.9996d0 double precision, parameter :: north=0.d0, east=500000.d0 @@ -90,96 +91,96 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway,SUPPRESS_UTM_PROJECT ! 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 + 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 + 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 diff --git a/src/shared/utm_geo.f90 b/src/shared/utm_geo.f90 index 0e663aea4..b59a7867a 100644 --- a/src/shared/utm_geo.f90 +++ b/src/shared/utm_geo.f90 @@ -41,6 +41,33 @@ ! 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 + +! original publication: +! Snyder, J.P., Map projections - a working manual, USGS professional paper 1395, 1987 +! https://doi.org/10.3133/pp1395 +! +! (Universal) Transverse Mercator projection - formulas for the ellipsoid (pages 61-63) are implemented +! +! note: the original version of this utm_geo(..) routine contained 2 errors in the formula +! in factor M (rm) when converting from long/lat to UTM, and +! in factor f3 for computing latitude (lat) when converting from UTM to long/lat. +! +! mistakes size in original routine: +! from UTM x/y to lon/lat and converted back again: +! 477415.5 5712313.5 -> 2.6741959317615298 51.561449486527074 (zone 31) +! 477415.50000327773 5712320.9415952200 <- 2.6741959317615298 51.561449486527074 +! has maximum UTM error ~ 7.44 m +! +! corrected version: +! 477415.5 5712313.5 -> 2.6741959317615298 51.561449479910003 (zone 31) +! 477415.49999999994 5712313.5002520066 <- 2.6741959317615298 51.561449479910003 +! has maximum UTM error ~ 0.00025 m +! +! error introduced in converting lon/lat to UTM: +! 2.6741959317615298 51.561449486527074 -> 477415.50000327773 5712320.9415952200 ! wrong formula +! 477415.50579994149 5712313.6646418981 ! corrected + + subroutine utm_geo(rlon4,rlat4,rx4,ry4,iway) ! @@ -139,16 +166,17 @@ subroutine utm_geo(rlon4,rlat4,rx4,ry4,iway) return endif -! save original parameters + ! save original parameters rlon_save = rlon4 rlat_save = rlat4 rx_save = rx4 ry_save = ry4 - e2 = 1.d0-(SEMI_MINOR_AXIS/SEMI_MAJOR_AXIS)**2 + ! defines parameters of reference ellipsoid + e2 = 1.d0 - (SEMI_MINOR_AXIS/SEMI_MAJOR_AXIS)**2 e4 = e2*e2 e6 = e2*e4 - ep2 = e2/(1.d0-e2) + ep2 = e2 / (1.d0-e2) ! !---- Set Zone parameters @@ -157,8 +185,10 @@ subroutine utm_geo(rlon4,rlat4,rx4,ry4,iway) lsouth = .false. if (UTM_PROJECTION_ZONE < 0) lsouth = .true. zone = abs(UTM_PROJECTION_ZONE) - cm = zone*6.0d0 - 183.d0 - cmr = cm*DEGREES_TO_RADIANS + + ! set central meridian for this zone + cm = zone * 6.0d0 - 183.d0 + cmr = cm * DEGREES_TO_RADIANS if (iway == IUTM2LONGLAT) then xx = rx4 @@ -180,36 +210,53 @@ subroutine utm_geo(rlon4,rlat4,rx4,ry4,iway) delam = dlon - cm if (delam < -180.d0) delam = delam + 360.d0 if (delam > 180.d0) delam = delam - 360.d0 + delam = delam*DEGREES_TO_RADIANS + ! Snyder, J.P., Map projections - a working manual, USGS professional paper 1395, 1987 + ! https://doi.org/10.3133/pp1395 + ! + ! (Universal) Transverse Mercator projection + ! + ! page 61, eq. 3-21 for M f1 = (1.d0 - e2/4.d0 - 3.d0*e4/64.d0 - 5.d0*e6/256d0)*rlat f2 = 3.d0*e2/8.d0 + 3.d0*e4/32.d0 + 45.d0*e6/1024.d0 - f2 = f2*sin(2.d0*rlat) - f3 = 15.d0*e4/256.d0*45.d0*e6/1024.d0 - f3 = f3*sin(4.d0*rlat) + f2 = f2 * sin(2.d0*rlat) + ! corrected: using .. + 45 e6 / 1024 instead of .. * 45 e6 / 1024 + !!wrong: f3 = 15.0 * e4 / 256.0 * 45.0 * e6 /1024.0 + f3 = 15.d0*e4/256.d0 + 45.d0*e6/1024.d0 + f3 = f3 * sin(4.d0*rlat) f4 = 35.d0*e6/3072.d0 - f4 = f4*sin(6.d0*rlat) + f4 = f4 * sin(6.d0*rlat) rm = SEMI_MAJOR_AXIS*(f1 - f2 + f3 - f4) + if (dlat == 90.d0 .or. dlat == -90.d0) then xx = 0.d0 yy = scfa*rm else + ! page 61, eq. 4-20 rn = SEMI_MAJOR_AXIS/sqrt(1.d0 - e2*sin(rlat)**2) + ! page 61, eq. 8-13 t = tan(rlat)**2 + ! page 61, eq. 8-14 c = ep2*cos(rlat)**2 a = cos(rlat)*delam - f1 = (1.d0 - t + c)*a**3/6.d0 + ! page 61, eq. 8-9 for x + f1 = (1.d0 - t + c) * a**3 / 6.d0 f2 = 5.d0 - 18.d0*t + t**2 + 72.d0*c - 58.d0*ep2 - f2 = f2*a**5/120.d0 - xx = scfa*rn*(a + f1 + f2) - f1 = a**2/2.d0 - f2 = 5.d0 - t + 9.d0*c + 4.d0*c**2 - f2 = f2*a**4/24.d0 + f2 = f2 * a**5 / 120.d0 + xx = scfa * rn * (a + f1 + f2) + + ! page 61, eq. 8-10 for y + f1 = a**2 / 2.d0 + f2 = 5.d0 - t + 9.d0*c + 4.d0 * c**2 + f2 = f2 * a**4 / 24.d0 f3 = 61.d0 - 58.d0*t + t**2 + 600.d0*c - 330.d0*ep2 - f3 = f3*a**6/720.d0 - yy = scfa*(rm + rn*tan(rlat)*(f1 + f2 + f3)) + f3 = f3 * a**6 / 720.d0 + yy = scfa * (rm + rn * tan(rlat) * (f1 + f2 + f3)) endif + xx = xx + east yy = yy + north @@ -220,20 +267,31 @@ subroutine utm_geo(rlon4,rlat4,rx4,ry4,iway) xx = xx - east yy = yy - north + + ! Snyder, J.P., Map projections - a working manual, USGS professional paper 1395, 1987 + ! https://doi.org/10.3133/pp1395 + ! + ! (Universal) Transverse Mercator projection + ! inverse formulas + ! + ! page 63, eq. 3-24 for e_1 e1 = sqrt(1.d0 - e2) e1 = (1.d0 - e1)/(1.d0 + e1) rm = yy/scfa + ! page 63, eq. 7-19 for mu u = 1.d0 - e2/4.d0 - 3.d0*e4/64.d0 - 5.d0*e6/256.d0 u = rm/(SEMI_MAJOR_AXIS*u) - f1 = 3.d0*e1/2.d0 - 27.d0*e1**3.d0/32.d0 - f1 = f1*sin(2.d0*u) - f2 = 21.d0*e1**2/16.d0 - 55.d0*e1**4/32.d0 - f2 = f2*sin(4.d0*u) - f3 = 151.d0*e1**3.d0/96.d0 - f3 = f3*sin(6.d0*u) + ! page 63, eq. 3-26 for phi_1 + f1 = 3.d0 * e1/2.d0 - 27.d0 * e1**3 / 32.d0 + f1 = f1 * sin(2.d0*u) + f2 = 21.d0 * e1**2 / 16.d0 - 55.d0 * e1**4 / 32.d0 + f2 = f2 * sin(4.d0*u) + f3 = 151.d0 * e1**3 / 96.d0 + f3 = f3 * sin(6.d0*u) rlat1 = u + f1 + f2 + f3 - dlat1 = rlat1*RADIANS_TO_DEGREES + dlat1 = rlat1 * RADIANS_TO_DEGREES + if (dlat1 >= 90.d0 .or. dlat1 <= -90.d0) then dlat1 = dmin1(dlat1,90.d0) dlat1 = dmax1(dlat1,-90.d0) @@ -246,21 +304,26 @@ subroutine utm_geo(rlon4,rlat4,rx4,ry4,iway) r1 = SEMI_MAJOR_AXIS*(1.d0 - e2)/sqrt(f1**3) d = xx/(rn1*scfa) - f1 = rn1*tan(rlat1)/r1 - f2 = d**2/2.d0 - f3 = 5.d0*3.d0*t1 + 10.d0*c1 - 4.d0*c1**2 - 9.d0*ep2 - f3 = f3*d**2*d**2/24.d0 + ! page 63, eq. 8-17 for phi + f1 = rn1 * tan(rlat1)/r1 + f2 = d**2 / 2.d0 + ! corrected: factor 5 + 3 T1 + .. instead of 5 * 3 T1 .. + !!wrong: f3 = 5.d0*3.d0*t1 + 10.d0*c1 - 4.d0*c1**2 - 9.d0*ep2 + f3 = 5.d0 + 3.d0*t1 + 10.d0*c1 - 4.d0*c1**2 - 9.d0*ep2 + f3 = f3 * d**4 / 24.d0 f4 = 61.d0 + 90.d0*t1 + 298.d0*c1 + 45.d0*t1**2 - 252.d0*ep2 - 3.d0*c1**2 - f4 = f4*(d**2)**3.d0/720.d0 + f4 = f4 * d**6 / 720.d0 rlat = rlat1 - f1*(f2 - f3 + f4) - dlat = rlat*RADIANS_TO_DEGREES + dlat = rlat * RADIANS_TO_DEGREES + ! page 63, eq. 8-18 for lambda f1 = 1.d0 + 2.d0*t1 + c1 - f1 = f1*d**2*d/6.d0 + f1 = f1 * d**3 / 6.d0 f2 = 5.d0 - 2.d0*c1 + 28.d0*t1 - 3.d0*c1**2 + 8.d0*ep2 + 24.d0*t1**2 - f2 = f2*(d**2)**2*d/120.d0 + f2 = f2 * d**5 / 120.d0 rlon = cmr + (d - f1 + f2)/cos(rlat1) - dlon = rlon*RADIANS_TO_DEGREES + dlon = rlon * RADIANS_TO_DEGREES + if (dlon < -180.d0) dlon = dlon + 360.d0 if (dlon > 180.d0) dlon = dlon - 360.d0 endif diff --git a/src/specfem3D/compute_coupling_acoustic_el.f90 b/src/specfem3D/compute_coupling_acoustic_el.f90 index 4937ed9e3..6361cdf21 100644 --- a/src/specfem3D/compute_coupling_acoustic_el.f90 +++ b/src/specfem3D/compute_coupling_acoustic_el.f90 @@ -46,29 +46,29 @@ subroutine compute_coupling_acoustic_el(NSPEC_AB,NGLOB_AB, & integer,intent(in) :: NSPEC_AB,NGLOB_AB -! displacement and pressure + ! displacement and pressure real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB),intent(in) :: displ real(kind=CUSTOM_REAL), dimension(NGLOB_AB),intent(inout) :: potential_dot_dot_acoustic integer,intent(in) :: SIMULATION_TYPE logical,intent(in) :: backward_simulation -! global indexing + ! global indexing integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: ibool -! acoustic-elastic coupling surface + ! acoustic-elastic coupling surface integer,intent(in) :: num_coupling_ac_el_faces real(kind=CUSTOM_REAL),intent(in) :: coupling_ac_el_normal(NDIM,NGLLSQUARE,num_coupling_ac_el_faces) real(kind=CUSTOM_REAL),intent(in) :: coupling_ac_el_jacobian2Dw(NGLLSQUARE,num_coupling_ac_el_faces) integer,intent(in) :: coupling_ac_el_ijk(3,NGLLSQUARE,num_coupling_ac_el_faces) integer,intent(in) :: coupling_ac_el_ispec(num_coupling_ac_el_faces) -! communication overlap + ! communication overlap integer,intent(in) :: iphase -! CPML + ! CPML logical,intent(in) :: PML_CONDITIONS -! local parameters + ! local parameters real(kind=CUSTOM_REAL) :: displ_x,displ_y,displ_z,displ_n real(kind=CUSTOM_REAL) :: nx,ny,nz,jacobianw integer :: iface,igll,ispec,iglob,ispec_CPML,i,j,k @@ -78,7 +78,7 @@ subroutine compute_coupling_acoustic_el(NSPEC_AB,NGLOB_AB, & ! checks if anything to do if (num_coupling_ac_el_faces == 0) return -! loops on all coupling faces + ! loops on all coupling faces do iface = 1,num_coupling_ac_el_faces ! gets corresponding elements diff --git a/src/specfem3D/compute_coupling_viscoelastic_ac.f90 b/src/specfem3D/compute_coupling_viscoelastic_ac.f90 index efe50a894..5d665f4cd 100644 --- a/src/specfem3D/compute_coupling_viscoelastic_ac.f90 +++ b/src/specfem3D/compute_coupling_viscoelastic_ac.f90 @@ -48,38 +48,38 @@ subroutine compute_coupling_viscoelastic_ac(NSPEC_AB,NGLOB_AB, & integer,intent(in) :: NSPEC_AB,NGLOB_AB,SIMULATION_TYPE logical,intent(in) :: backward_simulation -! displacement and pressure + ! displacement and pressure real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB),intent(inout) :: accel real(kind=CUSTOM_REAL), dimension(NGLOB_AB),intent(in) :: potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic -! global indexing + ! global indexing integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB),intent(in) :: ibool -! acoustic-elastic coupling surface + ! acoustic-elastic coupling surface integer,intent(in) :: num_coupling_ac_el_faces real(kind=CUSTOM_REAL),intent(in) :: coupling_ac_el_normal(NDIM,NGLLSQUARE,num_coupling_ac_el_faces) real(kind=CUSTOM_REAL),intent(in) :: coupling_ac_el_jacobian2Dw(NGLLSQUARE,num_coupling_ac_el_faces) integer,intent(in) :: coupling_ac_el_ijk(3,NGLLSQUARE,num_coupling_ac_el_faces) integer,intent(in) :: coupling_ac_el_ispec(num_coupling_ac_el_faces) -! communication overlap + ! communication overlap integer,intent(in) :: iphase -! local parameters + ! local parameters real(kind=CUSTOM_REAL) :: pressure_x,pressure_y,pressure_z real(kind=CUSTOM_REAL) :: nx,ny,nz,jacobianw integer :: iface,igll,ispec,iglob integer :: i,j,k -! CPML + ! CPML integer :: ispec_CPML logical :: PML_CONDITIONS ! only add these contributions in first pass if (iphase /= 1) return -! loops on all coupling faces + ! loops on all coupling faces do iface = 1,num_coupling_ac_el_faces ! gets corresponding spectral element diff --git a/utils/ADJOINT_TOMOGRAPHY_TOOLS/iterate_adj/SEM2D_iterate/gji_paper/codes_18_06Aug2006/wave2d_sub2.f90 b/utils/ADJOINT_TOMOGRAPHY_TOOLS/iterate_adj/SEM2D_iterate/gji_paper/codes_18_06Aug2006/wave2d_sub2.f90 index ba209e422..9be5319e4 100644 --- a/utils/ADJOINT_TOMOGRAPHY_TOOLS/iterate_adj/SEM2D_iterate/gji_paper/codes_18_06Aug2006/wave2d_sub2.f90 +++ b/utils/ADJOINT_TOMOGRAPHY_TOOLS/iterate_adj/SEM2D_iterate/gji_paper/codes_18_06Aug2006/wave2d_sub2.f90 @@ -635,7 +635,9 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway) double precision rx,ry,rlon,rlat double precision, parameter :: degrad=PI/180., raddeg=180./PI - double precision, parameter :: semimaj=6378206.4, semimin=6356583.8 + double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0 ! Clarke 1866 original uses this one + !double precision, parameter :: semimaj=6378137.0d0, semimin=6356752.314245d0 ! WGS-84 + double precision, parameter :: scfa=.9996 double precision, parameter :: north=0., east=500000. @@ -695,7 +697,7 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway) 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 = 15.*e4/256.*45.*e6/1024. ! note: original used here contains an error, formula should be corrected to .. + 45.0*e6/1024. f3 = f3*sin(4.*rlat) f4 = 35.*e6/3072. f4 = f4*sin(6.*rlat) @@ -770,7 +772,7 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway) ! latitude f1 = rn1*tan(rlat1) / r1 f2 = d**2 / 2. - f3 = 5.*3.*t1 + 10.*c1 - 4.*c1**2 - 9.*ep2 + f3 = 5.*3.*t1 + 10.*c1 - 4.*c1**2 - 9.*ep2 ! note: original contains an error, should be 5. + 3.t1 .. instead of 5. * .. 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. diff --git a/utils/ADJOINT_TOMOGRAPHY_TOOLS/iterate_adj/SEM2D_iterate/src/wave2d_sub2.f90 b/utils/ADJOINT_TOMOGRAPHY_TOOLS/iterate_adj/SEM2D_iterate/src/wave2d_sub2.f90 index 33d565146..239a029a2 100644 --- a/utils/ADJOINT_TOMOGRAPHY_TOOLS/iterate_adj/SEM2D_iterate/src/wave2d_sub2.f90 +++ b/utils/ADJOINT_TOMOGRAPHY_TOOLS/iterate_adj/SEM2D_iterate/src/wave2d_sub2.f90 @@ -647,7 +647,9 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway) double precision rx,ry,rlon,rlat double precision, parameter :: degrad=PI/180., raddeg=180./PI - double precision, parameter :: semimaj=6378206.4, semimin=6356583.8 + !double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0 ! Clarke 1866 + double precision, parameter :: semimaj=6378137.0d0, semimin=6356752.314245d0 ! WGS-84 + double precision, parameter :: scfa=.9996 double precision, parameter :: north=0., east=500000. @@ -707,7 +709,7 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway) 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 = 15.*e4/256. + 45.*e6/1024. f3 = f3*sin(4.*rlat) f4 = 35.*e6/3072. f4 = f4*sin(6.*rlat) @@ -782,7 +784,7 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway) ! latitude f1 = rn1*tan(rlat1) / r1 f2 = d**2 / 2. - f3 = 5.*3.*t1 + 10.*c1 - 4.*c1**2 - 9.*ep2 + 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. diff --git a/utils/ADJOINT_TOMOGRAPHY_TOOLS/iterate_adj/cluster/model_slice/utm_geo.f90 b/utils/ADJOINT_TOMOGRAPHY_TOOLS/iterate_adj/cluster/model_slice/utm_geo.f90 index 3bbc01a55..ff9179adb 100644 --- a/utils/ADJOINT_TOMOGRAPHY_TOOLS/iterate_adj/cluster/model_slice/utm_geo.f90 +++ b/utils/ADJOINT_TOMOGRAPHY_TOOLS/iterate_adj/cluster/model_slice/utm_geo.f90 @@ -41,7 +41,9 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway,SUPPRESS_UTM_PROJECT logical SUPPRESS_UTM_PROJECTION double precision, parameter :: degrad=PI/180.d0, raddeg=180.d0/PI - double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0 + !double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0 ! Clarke 1866 + double precision, parameter :: semimaj=6378137.0d0, semimin=6356752.314245d0 ! WGS-84 + double precision, parameter :: scfa=0.9996d0 double precision, parameter :: north=0.d0, east=500000.d0 @@ -101,7 +103,7 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway,SUPPRESS_UTM_PROJECT 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 = 15.*e4/256. + 45.*e6/1024. f3 = f3*sin(4.*rlat) f4 = 35.*e6/3072. f4 = f4*sin(6.*rlat) @@ -164,7 +166,7 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway,SUPPRESS_UTM_PROJECT f1 = rn1*tan(rlat1)/r1 f2 = d**2/2. - f3 = 5.*3.*t1 + 10.*c1 - 4.*c1**2 - 9.*ep2 + 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. diff --git a/utils/lib/utm_geo.f90 b/utils/lib/utm_geo.f90 index 603cd1d56..72e736322 100644 --- a/utils/lib/utm_geo.f90 +++ b/utils/lib/utm_geo.f90 @@ -40,7 +40,8 @@ subroutine utm_geo(rlon,rlat,rx,ry,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 :: semimaj=6378206.4d0, semimin=6356583.8d0 ! Clarke 1866 + double precision, parameter :: semimaj=6378137.0d0, semimin=6356752.314245d0 ! WGS-84 double precision, parameter :: scfa=.9996d0 double precision, parameter :: north=0.d0, east=500000.d0 @@ -100,7 +101,7 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway) 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 = 15.*e4/256. + 45.*e6/1024. f3 = f3*sin(4.*rlat) f4 = 35.*e6/3072. f4 = f4*sin(6.*rlat) @@ -163,7 +164,7 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway) f1 = rn1*tan(rlat1)/r1 f2 = d**2/2. - f3 = 5.*3.*t1 + 10.*c1 - 4.*c1**2 - 9.*ep2 + 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. diff --git a/utils/unused_routines/create_movie_GMT/utm_geo.f90 b/utils/unused_routines/create_movie_GMT/utm_geo.f90 index be0b23213..46a8e9ce4 100644 --- a/utils/unused_routines/create_movie_GMT/utm_geo.f90 +++ b/utils/unused_routines/create_movie_GMT/utm_geo.f90 @@ -41,7 +41,9 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway,SUPPRESS_UTM_PROJECT logical SUPPRESS_UTM_PROJECTION double precision, parameter :: degrad=PI/180., raddeg=180./PI - double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0 + !double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0 ! Clarke 1866 + double precision, parameter :: semimaj=6378137.0d0, semimin=6356752.314245d0 ! WGS-84 + double precision, parameter :: scfa=.9996d0 double precision, parameter :: north=0.d0, east=500000.d0 @@ -101,7 +103,7 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway,SUPPRESS_UTM_PROJECT 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 = 15.*e4/256. + 45.*e6/1024. f3 = f3*sin(4.*rlat) f4 = 35.*e6/3072. f4 = f4*sin(6.*rlat) @@ -164,7 +166,7 @@ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway,SUPPRESS_UTM_PROJECT f1 = rn1*tan(rlat1)/r1 f2 = d**2/2. - f3 = 5.*3.*t1 + 10.*c1 - 4.*c1**2 - 9.*ep2 + 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.