[cig-commits] r19159 - seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Tue Nov 8 11:19:55 PST 2011
Author: danielpeter
Date: 2011-11-08 11:19:55 -0800 (Tue, 08 Nov 2011)
New Revision: 19159
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90
Log:
adds choice for horizontal outputs in create_movie_AVS_DX.f90
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90 2011-11-08 18:42:40 UTC (rev 19158)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90 2011-11-08 19:19:55 UTC (rev 19159)
@@ -42,7 +42,8 @@
integer NSTEP,NTSTEP_BETWEEN_FRAMES,NCHUNKS
integer NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA
logical MOVIE_SURFACE
-
+ integer USE_COMPONENT
+
! ************** PROGRAM STARTS HERE **************
call read_AVS_DX_parameters(NEX_XI,NEX_ETA, &
@@ -55,7 +56,7 @@
print *,'1 = create files in OpenDX format'
print *,'2 = create files in AVS UCD format with individual files'
print *,'3 = create files in AVS UCD format with one time-dependent file'
- print *,'4 = create files in GMT xyz Ascii long/lat/Uz format'
+ print *,'4 = create files in GMT xyz Ascii long/lat/U format'
print *,'any other value = exit'
print *
print *,'enter value:'
@@ -71,12 +72,16 @@
print *,'enter last time step of movie (e.g. ',NSTEP,')'
read(5,*) it2
+ print *,'enter component (e.g. 1=Z, 2=N, 3=E)'
+ read(5,*) USE_COMPONENT
+ if( USE_COMPONENT < 1 .or. USE_COMPONENT > 3 ) stop 'component must be 1, 2 or 3'
+
! run the main program
call create_movie_AVS_DX(iformat,it1,it2, &
- NEX_XI,NEX_ETA, &
- NSTEP,NTSTEP_BETWEEN_FRAMES, &
- NCHUNKS, &
- NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA)
+ NEX_XI,NEX_ETA, &
+ NSTEP,NTSTEP_BETWEEN_FRAMES, &
+ NCHUNKS,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ USE_COMPONENT)
end program xcreate_movie_AVS_DX
@@ -84,8 +89,11 @@
!=====================================================================
!
- subroutine create_movie_AVS_DX(iformat,it1,it2,NEX_XI,NEX_ETA,NSTEP,NTSTEP_BETWEEN_FRAMES, &
- NCHUNKS,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA)
+ subroutine create_movie_AVS_DX(iformat,it1,it2, &
+ NEX_XI,NEX_ETA, &
+ NSTEP,NTSTEP_BETWEEN_FRAMES, &
+ NCHUNKS,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ USE_COMPONENT)
implicit none
@@ -123,8 +131,15 @@
real(kind=CUSTOM_REAL) xcoord,ycoord,zcoord,rval,thetaval,phival,lat,long
real(kind=CUSTOM_REAL) displx,disply,displz
real(kind=CUSTOM_REAL) normal_x,normal_y,normal_z
+ real(kind=CUSTOM_REAL) RRval,rhoval
+ real(kind=CUSTOM_REAL) thetahat_x,thetahat_y,thetahat_z
+ real(kind=CUSTOM_REAL) phihat_x,phihat_y
+
double precision min_field_current,max_field_current,max_absol
+
logical USE_OPENDX,UNIQUE_FILE,USE_GMT,USE_AVS
+ integer USE_COMPONENT
+
integer iformat,nframes,iframe
character(len=150) outputname
@@ -319,12 +334,18 @@
! read all the elements from the same file
write(outputname,"('/moviedata',i6.6)") it
open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='old',action='read',form='unformatted')
+
+ ! reads in point locations
+ ! (given as r theta phi for geocentric coordinate system)
read(IOUT) store_val_x
read(IOUT) store_val_y
read(IOUT) store_val_z
+
+ ! reads in associated values (displacement or velocity..)
read(IOUT) store_val_ux
read(IOUT) store_val_uy
read(IOUT) store_val_uz
+
close(IOUT)
! clear number of elements kept
@@ -343,6 +364,7 @@
ipoin = ipoin + 1
+ ! coordinates actually contain r theta phi
xcoord = store_val_x(ipoin,iproc)
ycoord = store_val_y(ipoin,iproc)
zcoord = store_val_z(ipoin,iproc)
@@ -351,23 +373,63 @@
disply = store_val_uy(ipoin,iproc)
displz = store_val_uz(ipoin,iproc)
-! coordinates actually contain r theta phi, therefore convert back to x y z
+ ! coordinates actually contain r theta phi, therefore convert back to x y z
rval = xcoord
thetaval = ycoord
phival = zcoord
call rthetaphi_2_xyz(xcoord,ycoord,zcoord,rval,thetaval,phival)
-! compute unit normal vector to the surface
- normal_x = xcoord / sqrt(xcoord**2 + ycoord**2 + zcoord**2)
- normal_y = ycoord / sqrt(xcoord**2 + ycoord**2 + zcoord**2)
- normal_z = zcoord / sqrt(xcoord**2 + ycoord**2 + zcoord**2)
-
-! save the results for this element
+ ! save the results for this element
x(i,j) = xcoord
y(i,j) = ycoord
z(i,j) = zcoord
- displn(i,j) = displx*normal_x + disply*normal_y + displz*normal_z
+ ! saves the desired component
+ if(USE_COMPONENT == 1) then
+ ! compute unit normal vector to the surface
+ RRval = sqrt(xcoord**2 + ycoord**2 + zcoord**2)
+ if( RRval < 1.e-10 ) stop 'error unit normal vector'
+ normal_x = xcoord / RRval
+ normal_y = ycoord / RRval
+ normal_z = zcoord / RRval
+
+ displn(i,j) = displx*normal_x + disply*normal_y + displz*normal_z
+
+ elseif(USE_COMPONENT == 2) then
+
+ ! compute unit tangent vector to the surface (N-S)
+ RRval = sqrt(xcoord**2 + ycoord**2 + zcoord**2)
+ if( RRval < 1.e-10 ) stop 'error unit normal vector'
+ rhoval = sqrt(xcoord**2 + ycoord**2)
+ if( rhoval < 1.e-10 ) then
+ ! location at pole
+ thetahat_x = 0.0
+ thetahat_y = 0.0
+ else
+ thetahat_x = (zcoord*xcoord) / (rhoval*RRval)
+ thetahat_y = (zcoord*ycoord) / (rhoval*RRval)
+ endif
+ thetahat_z = - rhoval/RRval
+
+ displn(i,j) = - (displx*thetahat_x + disply*thetahat_y + displz*thetahat_z)
+
+ elseif(USE_COMPONENT == 3) then
+
+ ! compute unit tangent to the surface (E-W)
+ rhoval = sqrt(xcoord**2 + ycoord**2)
+ if( rhoval < 1.e-10 ) then
+ ! location at pole
+ phihat_x = 0.0
+ phihat_y = 0.0
+ else
+ phihat_x = -ycoord / rhoval
+ phihat_y = xcoord / rhoval
+ endif
+
+ displn(i,j) = displx*phihat_x + disply*phihat_y
+ endif
+
+
enddo
enddo
@@ -490,24 +552,49 @@
! create file name and open file
if(USE_OPENDX) then
- write(outputname,"('/DX_movie_',i6.6,'.dx')") it
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
+ if( USE_COMPONENT == 1) then
+ write(outputname,"('/DX_movie_',i6.6,'.Z.dx')") it
+ elseif( USE_COMPONENT == 2) then
+ write(outputname,"('/DX_movie_',i6.6,'.N.dx')") it
+ elseif( USE_COMPONENT == 3) then
+ write(outputname,"('/DX_movie_',i6.6,'.E.dx')") it
+ endif
+ open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows'
else if(USE_AVS) then
if(UNIQUE_FILE .and. iframe == 1) then
- open(unit=11,file=trim(OUTPUT_FILES)//'/AVS_movie_all.inp',status='unknown')
+ if( USE_COMPONENT == 1) then
+ outputname = '/AVS_movie_all.Z.inp'
+ elseif( USE_COMPONENT == 2) then
+ outputname = '/AVS_movie_all.N.inp'
+ elseif( USE_COMPONENT == 3) then
+ outputname = '/AVS_movie_all.E.inp'
+ endif
+ open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
write(11,*) nframes
write(11,*) 'data'
write(11,"('step',i1,' image',i1)") 1,1
write(11,*) nglob,' ',nspectot_AVS_max
else if(.not. UNIQUE_FILE) then
- write(outputname,"('/AVS_movie_',i6.6,'.inp')") it
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
+ if( USE_COMPONENT == 1) then
+ write(outputname,"('/AVS_movie_',i6.6,'.Z.inp')") it
+ elseif( USE_COMPONENT == 2) then
+ write(outputname,"('/AVS_movie_',i6.6,'.N.inp')") it
+ elseif( USE_COMPONENT == 3) then
+ write(outputname,"('/AVS_movie_',i6.6,'.E.inp')") it
+ endif
+ open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
endif
else if(USE_GMT) then
- write(outputname,"('/gmt_movie_',i6.6,'.xyz')") it
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
+ if( USE_COMPONENT == 1) then
+ write(outputname,"('/gmt_movie_',i6.6,'.Z.xyz')") it
+ elseif( USE_COMPONENT == 2) then
+ write(outputname,"('/gmt_movie_',i6.6,'.N.xyz')") it
+ elseif( USE_COMPONENT == 3) then
+ write(outputname,"('/gmt_movie_',i6.6,'.E.xyz')") it
+ endif
+ open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
else
stop 'wrong output format selected'
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90 2011-11-08 18:42:40 UTC (rev 19158)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90 2011-11-08 19:19:55 UTC (rev 19159)
@@ -267,6 +267,7 @@
print *,'enter component (e.g. 1=Z, 2=N, 3=E)'
read(5,*) USE_COMPONENT
+ if( USE_COMPONENT < 1 .or. USE_COMPONENT > 3 ) stop 'component must be 1, 2 or 3'
print *,'enter output ascii (F) or binary (T)'
read(5,*) OUTPUT_BINARY
@@ -524,7 +525,7 @@
y(i,j) = ycoord
z(i,j) = zcoord
-
+ ! saves the desired component
if(USE_COMPONENT == 1) then
! compute unit normal vector to the surface
RRval = sqrt(xcoord**2 + ycoord**2 + zcoord**2)
@@ -552,6 +553,7 @@
thetahat_z = - rhoval/RRval
displn(i,j) = - (displx*thetahat_x + disply*thetahat_y + displz*thetahat_z)
+
elseif(USE_COMPONENT == 3) then
! compute unit tangent to the surface (E-W)
More information about the CIG-COMMITS
mailing list