[cig-commits] [commit] devel, master: fixed important topography bug detected by Bernhard Schubert: longitude blocks were inverted in topography/bathymetry file (6a567c2)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:01:08 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit 6a567c257818fc47b7548acded1fe9cd727aa49e
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Mon Mar 20 12:50:30 2006 +0000
fixed important topography bug detected by Bernhard Schubert: longitude blocks were inverted in topography/bathymetry file
>---------------------------------------------------------------
6a567c257818fc47b7548acded1fe9cd727aa49e
topo_bathy/convert_etopo2_bin2ascii_Sun.f90 | 2 ++
topo_bathy/image_topo_bathy_no_smooth.jpg | Bin 587385 -> 0 bytes
topo_bathy/image_topo_bathy_smooth_7.jpg | Bin 320697 -> 0 bytes
topo_bathy/smooth_topo_bathy_PPM_image.f90 | 28 +++++++++----------
topo_bathy/swap_topo_bathy_Sun.f90 | 40 ++++++++++++++++++++++++++++
5 files changed, 54 insertions(+), 16 deletions(-)
diff --git a/topo_bathy/convert_etopo2_bin2ascii_Sun.f90 b/topo_bathy/convert_etopo2_bin2ascii_Sun.f90
index 5971520..25f2ee1 100644
--- a/topo_bathy/convert_etopo2_bin2ascii_Sun.f90
+++ b/topo_bathy/convert_etopo2_bin2ascii_Sun.f90
@@ -37,6 +37,8 @@
print *,'min and max of original topography = ',minval(ibathy_topo),maxval(ibathy_topo)
! save the ASCII file
+!! DK DK beware, longitude blocks MUST be swapped
+ stop 'beware, longitude blocks MUST be swapped, see file swap_topo_bathy_Sun.f90 for details'
open(unit=13,file='topo_bathy_etopo2.dat',status='unknown')
do itopo_y=1,NY_BATHY
do itopo_x=1,NX_BATHY
diff --git a/topo_bathy/image_topo_bathy_no_smooth.jpg b/topo_bathy/image_topo_bathy_no_smooth.jpg
deleted file mode 100644
index e54c273..0000000
Binary files a/topo_bathy/image_topo_bathy_no_smooth.jpg and /dev/null differ
diff --git a/topo_bathy/image_topo_bathy_smooth_7.jpg b/topo_bathy/image_topo_bathy_smooth_7.jpg
deleted file mode 100644
index e1e61c7..0000000
Binary files a/topo_bathy/image_topo_bathy_smooth_7.jpg and /dev/null differ
diff --git a/topo_bathy/smooth_topo_bathy_PPM_image.f90 b/topo_bathy/smooth_topo_bathy_PPM_image.f90
index 7279f60..3572775 100644
--- a/topo_bathy/smooth_topo_bathy_PPM_image.f90
+++ b/topo_bathy/smooth_topo_bathy_PPM_image.f90
@@ -20,11 +20,11 @@
!--- read topography and bathymetry ASCII file, smooth it and create PNM image to check it
!
-!! DK DK this program uses a very simple window filter
-!! DK DK it works fine because the model has a very high resolution
-!! DK DK in principle it would be better to use a circular cap to take
-!! DK DK the notion of angular distance between points into account in the filter
-!! DK DK in practice though this simple filter works perfectly fine
+! this program uses a very simple window filter
+! it works fine because the model has a very high resolution
+! in principle it would be better to use a circular cap to take
+! the notion of angular distance between points into account in the filter
+! in practice though this simple filter works perfectly fine
implicit none
@@ -38,8 +38,7 @@
integer ibathy_topo(NX_BATHY,NY_BATHY)
integer ibathy_topo_ori(NX_BATHY,NY_BATHY)
- integer ivalue,icurrent_rec,ix,iy
- integer minvalue,maxvalue
+ integer ix,iy,minvalue,maxvalue
integer ix_current,iy_current,ix_min,ix_max,iy_min,iy_max,ix_value,iy_value
double precision value_sum,area_window
@@ -53,7 +52,7 @@
print *
print *,'reading topo file'
- open(unit=13,file='topo_bathy_etopo4_ori_from_etopo2_subsampled.dat',status='old')
+ open(unit=13,file='topo_bathy_etopo4_from_etopo2_subsampled.dat',status='old')
do iy=1,NY_BATHY
do ix=1,NX_BATHY
read(13,*) ibathy_topo_ori(ix,iy)
@@ -61,15 +60,14 @@
enddo
close(13)
-!! DK DK compute min and max before smoothing
+! compute min and max before smoothing
minvalue = minval(ibathy_topo_ori)
maxvalue = maxval(ibathy_topo_ori)
print *,'min and max of topography before smoothing = ',minvalue,maxvalue
!----
-!! DK DK smooth topography/bathymetry model
-
+! smooth topography/bathymetry model
if(SMOOTH_THE_MODEL) then
print *
@@ -131,16 +129,16 @@
!----
-!! DK DK compute min and max after smoothing
+! compute min and max after smoothing
minvalue = minval(ibathy_topo)
maxvalue = maxval(ibathy_topo)
print *,'min and max of topography after smoothing = ',minvalue,maxvalue
-!! DK DK save the smoothed model
+! save the smoothed model
if(SMOOTH_THE_MODEL) then
print *
print *,'saving the smoothed model'
- open(unit=13,file='topo_bathy_etopo4_smoothed_window7.dat',status='unknown')
+ open(unit=13,file='topo_bathy_etopo4_smoothed_window_7.dat',status='unknown')
do iy=1,NY_BATHY
do ix=1,NX_BATHY
write(13,*) ibathy_topo(ix,iy)
@@ -168,12 +166,10 @@
! creating the header
open(unit=27,file='image_topo_bathy.pnm',status='unknown')
write(27,100)
- write(27,101)
write(27,102) NX_BATHY,NY_BATHY
write(27,103)
100 format('P3')
- 101 format('# creator DK')
102 format(i6,' ',i6)
103 format('255')
104 format(i3)
diff --git a/topo_bathy/swap_topo_bathy_Sun.f90 b/topo_bathy/swap_topo_bathy_Sun.f90
new file mode 100644
index 0000000..d04d386
--- /dev/null
+++ b/topo_bathy/swap_topo_bathy_Sun.f90
@@ -0,0 +1,40 @@
+
+!---- swap longitude blocks in topography and bathymetry file once and for all
+!---- to switch from the convention used in the original file on the Sun
+!---- to the convention used in SPECFEM3D
+
+ program swap_topo_bathy_Sun
+
+ implicit none
+
+ include "../../constants.h"
+
+! use integer array to store values
+ integer ibathy_topo(NX_BATHY,NY_BATHY)
+
+ integer itopo_x,itopo_y
+
+ open(unit=13,file='topo_bathy_etopo4_from_etopo2_subsampled.dat',status='old')
+ do itopo_y=1,NY_BATHY
+ do itopo_x=1,NX_BATHY
+ read(13,*) ibathy_topo(itopo_x,itopo_y)
+ enddo
+ enddo
+ close(13)
+
+! blocks of longitude in [0,180[ and [180,360[ must be swapped
+! in the final file, itopo_x = 1 should correspond to longitude = 0
+! therefore one should see Africa on the left of the JPEG image of topography
+ open(unit=13,file='topo_bathy_etopo4_from_etopo2_subsampled_2.dat',status='unknown')
+ do itopo_y=1,NY_BATHY
+ do itopo_x=NX_BATHY/2+1,NX_BATHY
+ write(13,*) ibathy_topo(itopo_x,itopo_y)
+ enddo
+ do itopo_x=1,NX_BATHY/2
+ write(13,*) ibathy_topo(itopo_x,itopo_y)
+ enddo
+ enddo
+ close(13)
+
+ end program swap_topo_bathy_Sun
+
More information about the CIG-COMMITS
mailing list