[cig-commits] r22751 - in seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src: auxiliaries specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Sep 2 03:52:05 PDT 2013
Author: danielpeter
Date: 2013-09-02 03:52:04 -0700 (Mon, 02 Sep 2013)
New Revision: 22751
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/create_movie_AVS_DX.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/create_movie_GMT_global.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_stability.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic_calling_routine.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_output.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_surface.f90
Log:
updates movie surface routines in solver, and auxiliary programs create_movie_AVS_DX and create_movie_GMT_global to read in new moviedata format
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/create_movie_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/create_movie_AVS_DX.f90 2013-08-31 00:33:04 UTC (rev 22750)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/create_movie_AVS_DX.f90 2013-09-02 10:52:04 UTC (rev 22751)
@@ -32,74 +32,17 @@
program xcreate_movie_AVS_DX
- implicit none
-
- integer it1,it2
- integer iformat
-
-! parameters read from parameter file
- integer :: NEX_XI,NEX_ETA
- 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, &
- NSTEP,NTSTEP_BETWEEN_FRAMES, &
- NCHUNKS,MOVIE_SURFACE, &
- NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA)
-
- if(.not. MOVIE_SURFACE) stop 'movie frames were not saved by the solver'
-
- 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/U format'
- print *,'any other value = exit'
- print *
- print *,'enter value:'
- read(5,*) iformat
- if(iformat<1 .or. iformat>4) stop 'exiting...'
-
- print *,'movie frames have been saved every ',NTSTEP_BETWEEN_FRAMES,' time steps'
- print *
-
- print *,'enter first time step of movie (e.g. 1)'
- read(5,*) it1
-
- 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, &
- USE_COMPONENT)
-
- end program xcreate_movie_AVS_DX
-
-!
-!=====================================================================
-!
-
- 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)
-
use constants
+ use shared_parameters,only: &
+ NEX_XI,NEX_ETA,NSTEP,NTSTEP_BETWEEN_FRAMES, &
+ NCHUNKS,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ OUTPUT_FILES,MOVIE_COARSE
implicit none
!---------------------
+! USER PARAMETER
+
! threshold and normalization parameters
! single parameter to turn off all thresholding
@@ -110,8 +53,9 @@
! flag to apply non linear scaling to normalized norm of displacement
logical, parameter :: NONLINEAR_SCALING = .false.
+
! uses fixed max_value to normalize instead of max of current wavefield
- logical, parameter :: FIX_SCALING = .true.
+ logical, parameter :: FIX_SCALING = .false.
real,parameter:: MAX_VALUE = 1.0
! coefficient of power law used for non linear scaling
@@ -122,53 +66,89 @@
!---------------------
- integer i,j,it
- integer it1,it2
- integer nspectot_AVS_max
- integer ispec
- integer ibool_number,ibool_number1,ibool_number2,ibool_number3,ibool_number4
+ ! user input
+ integer :: it1,it2
+ integer :: iformat
+ integer :: USE_COMPONENT
+
+ integer :: i,j,it
+
+ integer :: nspectot_AVS_max
+ integer :: ispec
+ integer :: ibool_number,ibool_number1,ibool_number2,ibool_number3,ibool_number4
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: x,y,z,displn
- 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
+ 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
+ double precision :: min_field_current,max_field_current,max_absol
- logical USE_OPENDX,UNIQUE_FILE,USE_GMT,USE_AVS
- integer USE_COMPONENT
+ logical :: USE_OPENDX,UNIQUE_FILE,USE_GMT,USE_AVS
- integer iformat,nframes,iframe
+ integer :: nframes,iframe
- character(len=150) outputname
+ character(len=150) :: outputname
- integer iproc,ipoin
+ integer :: iproc,ipoin
! for sorting routine
- integer npointot,ilocnum,nglob,ielm,ieoff,ispecloc
+ integer :: npointot,ilocnum,nglob,ielm,ieoff,ispecloc
+ integer :: NIT
integer, dimension(:), allocatable :: iglob,loc,ireorder
logical, dimension(:), allocatable :: ifseg,mask_point
double precision, dimension(:), allocatable :: xp,yp,zp,xp_save,yp_save,zp_save,field_display
! for dynamic memory allocation
- integer ierror
+ integer :: ierror
! movie files stored by solver
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
store_val_x,store_val_y,store_val_z, &
store_val_ux,store_val_uy,store_val_uz
-! parameters read from file or deduced from parameters read from file
- integer NEX_XI,NEX_ETA
- integer NSTEP,NTSTEP_BETWEEN_FRAMES,NCHUNKS
- integer NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA
- character(len=150) OUTPUT_FILES
+! ************** PROGRAM STARTS HERE **************
+ call read_AVS_DX_parameters()
+
+ ! user input
+ 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/U format'
+ print *,'any other value = exit'
+ print *
+ print *,'enter value:'
+ read(5,*) iformat
+ if( iformat < 1 .or. iformat > 4 ) stop 'exiting...'
+
+ print *,'movie frames have been saved every ',NTSTEP_BETWEEN_FRAMES,' time steps'
+ print *
+
+ print *,'enter first time step of movie (e.g. 1)'
+ read(5,*) it1
+
+ print *,'enter last time step of movie (e.g. ',NSTEP,'or -1 for all)'
+ 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'
+
+ print *
+ print *,'--------'
+ print *
+
+ ! checks options
+ if( it2 == -1 ) it2 = NSTEP
+
! --------------------------------------
+ ! runs the main program
+ ! sets flags
if(iformat == 1) then
USE_OPENDX = .true.
USE_AVS = .false.
@@ -196,20 +176,28 @@
print *
print *,'Recombining all movie frames to create a movie'
print *
-
-! get the base pathname for output files
- call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
-
print *
print *,'There are ',NPROCTOT,' slices numbered from 0 to ',NPROCTOT-1
print *
- ilocnum = NGLLX * NGLLY * NEX_PER_PROC_XI * NEX_PER_PROC_ETA
+ if(MOVIE_COARSE) then
+ ! note:
+ ! nex_per_proc_xi*nex_per_proc_eta = nex_xi*nex_eta/nproc = nspec2d_top(iregion_crust_mantle)
+ ! used in specfem3D.f90
+ ! and ilocnum = nmovie_points = 2 * 2 * NEX_XI * NEX_ETA / NPROC
+ ilocnum = 2 * 2 * NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ NIT = NGLLX-1
+ else
+ ilocnum = NGLLX*NGLLY*NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ NIT = 1
+ endif
print *
- print *,'Allocating arrays of size ',ilocnum*NPROCTOT
+ print *,'Allocating arrays for reading data of size ',ilocnum*NPROCTOT,'=', &
+ 6*ilocnum*NPROCTOT*CUSTOM_REAL/1024./1024.,'MB'
print *
+
allocate(store_val_x(ilocnum,0:NPROCTOT-1),stat=ierror)
if(ierror /= 0) stop 'error while allocating store_val_x'
@@ -257,13 +245,17 @@
! spectral element one should have (NGLL-1)^2 OpenDX "elements"
! define the total number of OpenDX "elements" at the surface
- nspectot_AVS_max = NCHUNKS * NEX_XI * NEX_ETA * (NGLLX-1) * (NGLLY-1)
+ if(MOVIE_COARSE) then
+ nspectot_AVS_max = NCHUNKS * NEX_XI * NEX_ETA
+ else
+ nspectot_AVS_max = NCHUNKS * NEX_XI * NEX_ETA * (NGLLX-1) * (NGLLY-1)
+ endif
print *
print *,'there are a total of ',nspectot_AVS_max,' OpenDX "elements" at the surface'
print *
-! maximum theoretical number of points at the surface
+ ! maximum theoretical number of points at the surface
npointot = NGNOD2D_AVS_DX * nspectot_AVS_max
print *
@@ -317,26 +309,12 @@
! --------------------------------------
- iframe = 0
-
-! loop on all the time steps in the range entered
- do it = it1,it2
-
-! check if time step corresponds to a movie frame
- if(mod(it,NTSTEP_BETWEEN_FRAMES) == 0) then
-
- iframe = iframe + 1
-
- print *
- print *,'reading snapshot time step ',it,' out of ',NSTEP
- print *
-
-! read all the elements from the same file
- write(outputname,"('/moviedata',i6.6)") it
- open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname, &
+ ! movie point locations
+ outputname = "/moviedata_xyz.bin"
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//trim(outputname), &
status='old',action='read',form='unformatted',iostat=ierror)
if(ierror /= 0 ) then
- print*,'error opening file: ',trim(OUTPUT_FILES)//outputname
+ print*,'error opening file: ',trim(OUTPUT_FILES)//trim(outputname)
stop 'error opening moviedata file'
endif
@@ -345,434 +323,507 @@
read(IOUT) store_val_x
read(IOUT) store_val_y
read(IOUT) store_val_z
+ close(IOUT)
- ! reads in associated values (displacement or velocity..)
- read(IOUT) store_val_ux
- read(IOUT) store_val_uy
- read(IOUT) store_val_uz
+! --------------------------------------
- close(IOUT)
+ iframe = 0
- !debug
- print*
- print*,"data ux min/max: ",minval(store_val_ux),maxval(store_val_ux)
- print*,"data uy min/max: ",minval(store_val_uy),maxval(store_val_uy)
- print*,"data uz min/max: ",minval(store_val_uz),maxval(store_val_uz)
- print*
+ ! loop on all the time steps in the range entered
+ do it = it1,it2
-! clear number of elements kept
- ispec = 0
+ ! check if time step corresponds to a movie frame
+ if(mod(it,NTSTEP_BETWEEN_FRAMES) /= 0) cycle
-! read points for all the slices
- do iproc = 0,NPROCTOT-1
+ iframe = iframe + 1
-! reset point number
- ipoin = 0
+ print *
+ print *,'reading snapshot time step ',it,' out of ',NSTEP
+ print *
- do ispecloc = 1,NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ ! 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',iostat=ierror)
+ if(ierror /= 0 ) then
+ print*,'error opening file: ',trim(OUTPUT_FILES)//outputname
+ stop 'error opening moviedata file'
+ endif
- do j = 1,NGLLY
- do i = 1,NGLLX
+ ! reads in associated values (displacement or velocity..)
+ read(IOUT) store_val_ux
+ read(IOUT) store_val_uy
+ read(IOUT) store_val_uz
- ipoin = ipoin + 1
+ close(IOUT)
- ! coordinates actually contain r theta phi
- xcoord = store_val_x(ipoin,iproc)
- ycoord = store_val_y(ipoin,iproc)
- zcoord = store_val_z(ipoin,iproc)
+ print *,' file ',trim(OUTPUT_FILES)//trim(outputname)
+ print *
+ !debug
+ print *," data ux min/max: ",minval(store_val_ux),maxval(store_val_ux)
+ print *," data uy min/max: ",minval(store_val_uy),maxval(store_val_uy)
+ print *," data uz min/max: ",minval(store_val_uz),maxval(store_val_uz)
+ print *
- displx = store_val_ux(ipoin,iproc)
- disply = store_val_uy(ipoin,iproc)
- displz = store_val_uz(ipoin,iproc)
+ ! clear number of elements kept
+ ispec = 0
- ! 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)
+ ! read points for all the slices
+ do iproc = 0,NPROCTOT-1
- ! save the results for this element
- x(i,j) = xcoord
- y(i,j) = ycoord
- z(i,j) = zcoord
+ ! reset point number
+ ipoin = 0
- ! 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
+ do ispecloc = 1,NEX_PER_PROC_XI*NEX_PER_PROC_ETA
- displn(i,j) = displx*normal_x + disply*normal_y + displz*normal_z
+ ! gets element values
+ do j = 1,NGLLY,NIT
+ do i = 1,NGLLX,NIT
+ ipoin = ipoin + 1
- else if(USE_COMPONENT == 2) then
+ ! coordinates actually contain r theta phi
+ xcoord = store_val_x(ipoin,iproc)
+ ycoord = store_val_y(ipoin,iproc)
+ zcoord = store_val_z(ipoin,iproc)
- ! 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
+ displx = store_val_ux(ipoin,iproc)
+ disply = store_val_uy(ipoin,iproc)
+ displz = store_val_uz(ipoin,iproc)
- displn(i,j) = - (displx*thetahat_x + disply*thetahat_y + displz*thetahat_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)
- else if(USE_COMPONENT == 3) then
+ ! save the results for this element
+ x(i,j) = xcoord
+ y(i,j) = ycoord
+ z(i,j) = zcoord
- ! 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
+ ! 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*phihat_x + disply*phihat_y
- endif
+ displn(i,j) = displx*normal_x + disply*normal_y + displz*normal_z
+ else if(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)
+
+ else if(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
- enddo
-! assign the values of the corners of the OpenDX "elements"
- ispec = ispec + 1
- ielm = (NGLLX-1)*(NGLLY-1)*(ispec-1)
- do j = 1,NGLLY-1
- do i = 1,NGLLX-1
- ieoff = NGNOD2D_AVS_DX*(ielm+(i-1)+(j-1)*(NGLLX-1))
- do ilocnum = 1,NGNOD2D_AVS_DX
- if(ilocnum == 1) then
- xp(ieoff+ilocnum) = dble(x(i,j))
- yp(ieoff+ilocnum) = dble(y(i,j))
- zp(ieoff+ilocnum) = dble(z(i,j))
- field_display(ieoff+ilocnum) = dble(displn(i,j))
- else if(ilocnum == 2) then
- xp(ieoff+ilocnum) = dble(x(i+1,j))
- yp(ieoff+ilocnum) = dble(y(i+1,j))
- zp(ieoff+ilocnum) = dble(z(i+1,j))
- field_display(ieoff+ilocnum) = dble(displn(i+1,j))
- else if(ilocnum == 3) then
- xp(ieoff+ilocnum) = dble(x(i+1,j+1))
- yp(ieoff+ilocnum) = dble(y(i+1,j+1))
- zp(ieoff+ilocnum) = dble(z(i+1,j+1))
- field_display(ieoff+ilocnum) = dble(displn(i+1,j+1))
+ ! assign the values of the corners of the OpenDX "elements"
+ ispec = ispec + 1
+ if(MOVIE_COARSE) then
+ ielm = ispec-1
+ else
+ ielm = (NGLLX-1)*(NGLLY-1)*(ispec-1)
+ endif
+
+ do j = 1,NGLLY-NIT
+ do i = 1,NGLLX-NIT
+ ! offset (starts at 0)
+ if(MOVIE_COARSE) then
+ ieoff = NGNOD2D_AVS_DX * ielm
else
- xp(ieoff+ilocnum) = dble(x(i,j+1))
- yp(ieoff+ilocnum) = dble(y(i,j+1))
- zp(ieoff+ilocnum) = dble(z(i,j+1))
- field_display(ieoff+ilocnum) = dble(displn(i,j+1))
+ ieoff = NGNOD2D_AVS_DX * (ielm+(i-1)+(j-1)*(NGLLX-1))
endif
+
+ ! stores 4 points as element corners
+ do ilocnum = 1,NGNOD2D_AVS_DX
+ ! for movie_coarse e.g. x(i,j) is defined at x(1,1), x(1,NGLLY), x(NGLLX,1) and x(NGLLX,NGLLY)
+ ! be aware that for the cubed sphere, the mapping changes for different chunks,
+ ! i.e. e.g. x(1,1) and x(5,5) flip left and right sides of the elements in geographical coordinates
+ if(MOVIE_COARSE) then
+ if(ilocnum == 1) then
+ xp(ieoff+ilocnum) = dble(x(1,1))
+ yp(ieoff+ilocnum) = dble(y(1,1))
+ zp(ieoff+ilocnum) = dble(z(1,1))
+ field_display(ieoff+ilocnum) = dble(displn(1,1))
+ else if(ilocnum == 2) then
+ xp(ieoff+ilocnum) = dble(x(NGLLX,1))
+ yp(ieoff+ilocnum) = dble(y(NGLLX,1))
+ zp(ieoff+ilocnum) = dble(z(NGLLX,1))
+ field_display(ieoff+ilocnum) = dble(displn(NGLLX,1))
+ else if(ilocnum == 3) then
+ xp(ieoff+ilocnum) = dble(x(NGLLX,NGLLY))
+ yp(ieoff+ilocnum) = dble(y(NGLLX,NGLLY))
+ zp(ieoff+ilocnum) = dble(z(NGLLX,NGLLY))
+ field_display(ieoff+ilocnum) = dble(displn(NGLLX,NGLLY))
+ else
+ xp(ieoff+ilocnum) = dble(x(1,NGLLY))
+ yp(ieoff+ilocnum) = dble(y(1,NGLLY))
+ zp(ieoff+ilocnum) = dble(z(1,NGLLY))
+ field_display(ieoff+ilocnum) = dble(displn(1,NGLLY))
+ endif
+ else
+ ! movie saved on fine grid (all gll points)
+ if(ilocnum == 1) then
+ xp(ieoff+ilocnum) = dble(x(i,j))
+ yp(ieoff+ilocnum) = dble(y(i,j))
+ zp(ieoff+ilocnum) = dble(z(i,j))
+ field_display(ieoff+ilocnum) = dble(displn(i,j))
+ else if(ilocnum == 2) then
+ xp(ieoff+ilocnum) = dble(x(i+1,j))
+ yp(ieoff+ilocnum) = dble(y(i+1,j))
+ zp(ieoff+ilocnum) = dble(z(i+1,j))
+ field_display(ieoff+ilocnum) = dble(displn(i+1,j))
+ else if(ilocnum == 3) then
+ xp(ieoff+ilocnum) = dble(x(i+1,j+1))
+ yp(ieoff+ilocnum) = dble(y(i+1,j+1))
+ zp(ieoff+ilocnum) = dble(z(i+1,j+1))
+ field_display(ieoff+ilocnum) = dble(displn(i+1,j+1))
+ else
+ xp(ieoff+ilocnum) = dble(x(i,j+1))
+ yp(ieoff+ilocnum) = dble(y(i,j+1))
+ zp(ieoff+ilocnum) = dble(z(i,j+1))
+ field_display(ieoff+ilocnum) = dble(displn(i,j+1))
+ endif
+ endif ! MOVIE_COARSE
+ enddo ! ilocnum
+
enddo
enddo
- enddo
+ enddo ! ispecloc
+
enddo
- enddo
+ !---------------------------------
+ ! apply normalization and thresholding, if desired
-!---------------------------------
-! apply normalization and thresholding, if desired
+ if (.not. ALL_THRESHOLD_OFF) then
- if (.not. ALL_THRESHOLD_OFF) then
+ ! compute min and max of data value to normalize
+ min_field_current = minval(field_display(:))
+ max_field_current = maxval(field_display(:))
- ! compute min and max of data value to normalize
- min_field_current = minval(field_display(:))
- max_field_current = maxval(field_display(:))
+ ! make sure range is always symmetric and center is in zero
+ ! this assumption works only for fields that can be negative
+ ! would not work for norm of vector for instance
+ ! (we would lose half of the color palette if no negative values)
+ max_absol = max(abs(min_field_current),abs(max_field_current))
+ min_field_current = - max_absol
+ max_field_current = + max_absol
- ! make sure range is always symmetric and center is in zero
- ! this assumption works only for fields that can be negative
- ! would not work for norm of vector for instance
- ! (we would lose half of the color palette if no negative values)
- max_absol = max(abs(min_field_current),abs(max_field_current))
- min_field_current = - max_absol
- max_field_current = + max_absol
+ ! print minimum and maximum amplitude in current snapshot
+ print *
+ print *,'minimum amplitude in current snapshot = ',min_field_current
+ print *,'maximum amplitude in current snapshot = ',max_field_current
+ if( FIX_SCALING ) then
+ print *,' to be normalized by : ',MAX_VALUE
+ if( max_field_current > MAX_VALUE ) stop 'increase MAX_VALUE'
+ endif
+ print *
- ! print minimum and maximum amplitude in current snapshot
- print *
- print *,'minimum amplitude in current snapshot = ',min_field_current
- print *,'maximum amplitude in current snapshot = ',max_field_current
- if( FIX_SCALING ) then
- print *,' to be normalized by : ',MAX_VALUE
- if( max_field_current > MAX_VALUE ) stop 'increase MAX_VALUE'
- endif
- print *
+ ! normalize field to [0:1]
+ print *,'normalizing... '
+ if( FIX_SCALING ) then
+ field_display(:) = (field_display(:) + MAX_VALUE) / (2.0*MAX_VALUE)
+ else
+ field_display(:) = (field_display(:) - min_field_current) / (max_field_current - min_field_current)
+ endif
+ ! rescale to [-1,1]
+ field_display(:) = 2.*field_display(:) - 1.
- ! normalize field to [0:1]
- print *,'normalizing... '
- if( FIX_SCALING ) then
- field_display(:) = (field_display(:) + MAX_VALUE) / (2.0*MAX_VALUE)
- else
- field_display(:) = (field_display(:) - min_field_current) / (max_field_current - min_field_current)
- endif
- ! rescale to [-1,1]
- field_display(:) = 2.*field_display(:) - 1.
+ ! apply threshold to normalized field
+ if(APPLY_THRESHOLD) then
+ print *,'thresholding... '
+ where(abs(field_display(:)) <= THRESHOLD) field_display = 0.
+ endif
- ! apply threshold to normalized field
- if(APPLY_THRESHOLD) then
- print *,'thresholding... '
- where(abs(field_display(:)) <= THRESHOLD) field_display = 0.
- endif
+ ! apply non linear scaling to normalized field if needed
+ if(NONLINEAR_SCALING) then
+ print *,'nonlinear scaling... '
+ where(field_display(:) >= 0.)
+ field_display = field_display ** POWER_SCALING
+ elsewhere
+ field_display = - abs(field_display) ** POWER_SCALING
+ endwhere
+ endif
- ! apply non linear scaling to normalized field if needed
- if(NONLINEAR_SCALING) then
- print *,'nonlinear scaling... '
- where(field_display(:) >= 0.)
- field_display = field_display ** POWER_SCALING
- elsewhere
- field_display = - abs(field_display) ** POWER_SCALING
- endwhere
- endif
+ print *,'color scaling... '
+ ! map back to [0,1]
+ field_display(:) = (field_display(:) + 1.) / 2.
- print *,'color scaling... '
- ! map back to [0,1]
- field_display(:) = (field_display(:) + 1.) / 2.
+ ! map field to [0:255] for AVS color scale
+ field_display(:) = 255. * field_display(:)
- ! map field to [0:255] for AVS color scale
- field_display(:) = 255. * field_display(:)
+ else
- else
+ ! compute min and max of data value
+ min_field_current = minval(field_display(:))
+ max_field_current = maxval(field_display(:))
- ! compute min and max of data value
- min_field_current = minval(field_display(:))
- max_field_current = maxval(field_display(:))
+ ! print minimum and maximum amplitude in current snapshot
+ print *
+ print *,'minimum amplitude in current snapshot = ',min_field_current
+ print *,'maximum amplitude in current snapshot = ',max_field_current
+ print *
- ! print minimum and maximum amplitude in current snapshot
- print *
- print *,'minimum amplitude in current snapshot = ',min_field_current
- print *,'maximum amplitude in current snapshot = ',max_field_current
- print *
+ endif
- endif
+ !---------------------------------
-!---------------------------------
+ ! copy coordinate arrays since the sorting routine does not preserve them
+ print *,'sorting... '
+ xp_save(:) = xp(:)
+ yp_save(:) = yp(:)
+ zp_save(:) = zp(:)
-! copy coordinate arrays since the sorting routine does not preserve them
- print *,'sorting... '
- xp_save(:) = xp(:)
- yp_save(:) = yp(:)
- zp_save(:) = zp(:)
+ !--- sort the list based upon coordinates to get rid of multiples
+ print *,'sorting list of points'
+ call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot)
-!--- sort the list based upon coordinates to get rid of multiples
- print *,'sorting list of points'
- call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot)
+ !--- print total number of points found
+ print *
+ print *,'found a total of ',nglob,' points'
+ print *,'initial number of points (with multiples) was ',npointot
-!--- print total number of points found
- print *
- print *,'found a total of ',nglob,' points'
- print *,'initial number of points (with multiples) was ',npointot
+ !--- ****** create AVS file using sorted list ******
-!--- ****** create AVS file using sorted list ******
-
-! create file name and open file
- if(USE_OPENDX) then
- if( USE_COMPONENT == 1) then
- write(outputname,"('/DX_movie_',i6.6,'.Z.dx')") it
- else if( USE_COMPONENT == 2) then
- write(outputname,"('/DX_movie_',i6.6,'.N.dx')") it
- else if( 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
+ ! create file name and open file
+ if(USE_OPENDX) then
if( USE_COMPONENT == 1) then
- outputname = '/AVS_movie_all.Z.inp'
+ write(outputname,"('/DX_movie_',i6.6,'.Z.dx')") it
else if( USE_COMPONENT == 2) then
- outputname = '/AVS_movie_all.N.inp'
+ write(outputname,"('/DX_movie_',i6.6,'.N.dx')") it
else if( USE_COMPONENT == 3) then
- outputname = '/AVS_movie_all.E.inp'
+ write(outputname,"('/DX_movie_',i6.6,'.E.dx')") it
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(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
+ if( USE_COMPONENT == 1) then
+ outputname = '/AVS_movie_all.Z.inp'
+ else if( USE_COMPONENT == 2) then
+ outputname = '/AVS_movie_all.N.inp'
+ else if( 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
+ if( USE_COMPONENT == 1) then
+ write(outputname,"('/AVS_movie_',i6.6,'.Z.inp')") it
+ else if( USE_COMPONENT == 2) then
+ write(outputname,"('/AVS_movie_',i6.6,'.N.inp')") it
+ else if( 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
if( USE_COMPONENT == 1) then
- write(outputname,"('/AVS_movie_',i6.6,'.Z.inp')") it
+ write(outputname,"('/gmt_movie_',i6.6,'.Z.xyz')") it
else if( USE_COMPONENT == 2) then
- write(outputname,"('/AVS_movie_',i6.6,'.N.inp')") it
+ write(outputname,"('/gmt_movie_',i6.6,'.N.xyz')") it
else if( USE_COMPONENT == 3) then
- write(outputname,"('/AVS_movie_',i6.6,'.E.inp')") it
+ write(outputname,"('/gmt_movie_',i6.6,'.E.xyz')") it
endif
open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
- write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
+ else
+ stop 'wrong output format selected'
endif
- else if(USE_GMT) then
- if( USE_COMPONENT == 1) then
- write(outputname,"('/gmt_movie_',i6.6,'.Z.xyz')") it
- else if( USE_COMPONENT == 2) then
- write(outputname,"('/gmt_movie_',i6.6,'.N.xyz')") it
- else if( 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
- if(USE_GMT) then
+ if(USE_GMT) then
- ! output list of points
- mask_point = .false.
- do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
- ! four points for each element
- do ilocnum = 1,NGNOD2D_AVS_DX
- ibool_number = iglob(ilocnum+ieoff)
- if(.not. mask_point(ibool_number)) then
- xcoord = sngl(xp_save(ilocnum+ieoff))
- ycoord = sngl(yp_save(ilocnum+ieoff))
- zcoord = sngl(zp_save(ilocnum+ieoff))
- call xyz_2_rthetaphi(xcoord,ycoord,zcoord,rval,thetaval,phival)
-
- ! note: converts the geocentric colatitude to a geographic colatitude
- if(.not. ASSUME_PERFECT_SPHERE) then
- thetaval = PI_OVER_TWO - &
- datan(1.006760466d0*dcos(dble(thetaval))/dmax1(TINYVAL,dble(sin(thetaval))))
- endif
-
- lat = (PI_OVER_TWO-thetaval)*RADIANS_TO_DEGREES
- long = phival*RADIANS_TO_DEGREES
- if(long > 180.0) long = long-360.0
- write(11,*) long,lat,sngl(field_display(ilocnum+ieoff))
- endif
- mask_point(ibool_number) = .true.
- enddo
- enddo
-
- else
-! if unique file, output geometry only once
- if(.not. UNIQUE_FILE .or. iframe == 1) then
-
-! output list of points
+ ! output list of points
mask_point = .false.
- ipoin = 0
do ispec=1,nspectot_AVS_max
ieoff = NGNOD2D_AVS_DX*(ispec-1)
-! four points for each element
+ ! four points for each element
do ilocnum = 1,NGNOD2D_AVS_DX
ibool_number = iglob(ilocnum+ieoff)
if(.not. mask_point(ibool_number)) then
- ipoin = ipoin + 1
- ireorder(ibool_number) = ipoin
- if(USE_OPENDX) then
- write(11,"(f10.7,1x,f10.7,1x,f10.7)") &
- xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
- else if(USE_AVS) then
- write(11,"(i10,1x,f10.7,1x,f10.7,1x,f10.7)") ireorder(ibool_number), &
- xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
+ ! gets cartesian coordinates
+ xcoord = sngl(xp_save(ilocnum+ieoff))
+ ycoord = sngl(yp_save(ilocnum+ieoff))
+ zcoord = sngl(zp_save(ilocnum+ieoff))
+
+ ! location latitude/longitude (with geocentric colatitude theta )
+ call xyz_2_rthetaphi(xcoord,ycoord,zcoord,rval,thetaval,phival)
+
+ ! note: converts the geocentric colatitude to a geographic colatitude
+ if(.not. ASSUME_PERFECT_SPHERE) then
+ thetaval = PI_OVER_TWO - &
+ datan(1.006760466d0*dcos(dble(thetaval))/dmax1(TINYVAL,dble(sin(thetaval))))
endif
+
+ ! gets geographic latitude and longitude in degrees
+ lat = (PI_OVER_TWO-thetaval)*RADIANS_TO_DEGREES
+ long = phival*RADIANS_TO_DEGREES
+ if(long > 180.0) long = long-360.0
+
+ write(11,*) long,lat,sngl(field_display(ilocnum+ieoff))
endif
mask_point(ibool_number) = .true.
enddo
enddo
- if(USE_OPENDX) &
- write(11,*) 'object 2 class array type int rank 1 shape 4 items ',nspectot_AVS_max,' data follows'
+ else
+ ! if unique file, output geometry only once
+ if(.not. UNIQUE_FILE .or. iframe == 1) then
-! output list of elements
- do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
-! four points for each element
- ibool_number1 = iglob(ieoff + 1)
- ibool_number2 = iglob(ieoff + 2)
- ibool_number3 = iglob(ieoff + 3)
- ibool_number4 = iglob(ieoff + 4)
- if(USE_OPENDX) then
-! point order in OpenDX is 1,4,2,3 *not* 1,2,3,4 as in AVS
- write(11,"(i10,1x,i10,1x,i10,1x,i10)") ireorder(ibool_number1)-1, &
- ireorder(ibool_number4)-1,ireorder(ibool_number2)-1,ireorder(ibool_number3)-1
- else
- write(11,"(i10,' 1 quad ',i10,1x,i10,1x,i10,1x,i10)") ispec,ireorder(ibool_number1), &
- ireorder(ibool_number2),ireorder(ibool_number3),ireorder(ibool_number4)
- endif
- enddo
+ ! output list of points
+ mask_point = .false.
+ ipoin = 0
+ do ispec=1,nspectot_AVS_max
+ ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ! four points for each element
+ do ilocnum = 1,NGNOD2D_AVS_DX
+ ibool_number = iglob(ilocnum+ieoff)
+ if(.not. mask_point(ibool_number)) then
+ ipoin = ipoin + 1
+ ireorder(ibool_number) = ipoin
+ if(USE_OPENDX) then
+ write(11,"(f10.7,1x,f10.7,1x,f10.7)") &
+ xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
+ else if(USE_AVS) then
+ write(11,"(i10,1x,f10.7,1x,f10.7,1x,f10.7)") ireorder(ibool_number), &
+ xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
+ endif
+ endif
+ mask_point(ibool_number) = .true.
+ enddo
+ enddo
- endif
+ if(USE_OPENDX) &
+ write(11,*) 'object 2 class array type int rank 1 shape 4 items ',nspectot_AVS_max,' data follows'
- if(USE_OPENDX) then
- write(11,*) 'attribute "element type" string "quads"'
- write(11,*) 'attribute "ref" string "positions"'
- write(11,*) 'object 3 class array type float rank 0 items ',nglob,' data follows'
- else
- if(UNIQUE_FILE) then
-! step number for AVS multistep file
- if(iframe > 1) then
- if(iframe < 10) then
- write(11,"('step',i1,' image',i1)") iframe,iframe
- else if(iframe < 100) then
- write(11,"('step',i2,' image',i2)") iframe,iframe
- else if(iframe < 1000) then
- write(11,"('step',i3,' image',i3)") iframe,iframe
+ ! output list of elements
+ do ispec=1,nspectot_AVS_max
+ ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ! four points for each element
+ ibool_number1 = iglob(ieoff + 1)
+ ibool_number2 = iglob(ieoff + 2)
+ ibool_number3 = iglob(ieoff + 3)
+ ibool_number4 = iglob(ieoff + 4)
+ if(USE_OPENDX) then
+ ! point order in OpenDX is 1,4,2,3 *not* 1,2,3,4 as in AVS
+ write(11,"(i10,1x,i10,1x,i10,1x,i10)") ireorder(ibool_number1)-1, &
+ ireorder(ibool_number4)-1,ireorder(ibool_number2)-1,ireorder(ibool_number3)-1
else
- write(11,"('step',i4,' image',i4)") iframe,iframe
+ write(11,"(i10,' 1 quad ',i10,1x,i10,1x,i10,1x,i10)") ispec,ireorder(ibool_number1), &
+ ireorder(ibool_number2),ireorder(ibool_number3),ireorder(ibool_number4)
endif
+ enddo
+
+ endif
+
+ if(USE_OPENDX) then
+ write(11,*) 'attribute "element type" string "quads"'
+ write(11,*) 'attribute "ref" string "positions"'
+ write(11,*) 'object 3 class array type float rank 0 items ',nglob,' data follows'
+ else
+ if(UNIQUE_FILE) then
+ ! step number for AVS multistep file
+ if(iframe > 1) then
+ if(iframe < 10) then
+ write(11,"('step',i1,' image',i1)") iframe,iframe
+ else if(iframe < 100) then
+ write(11,"('step',i2,' image',i2)") iframe,iframe
+ else if(iframe < 1000) then
+ write(11,"('step',i3,' image',i3)") iframe,iframe
+ else
+ write(11,"('step',i4,' image',i4)") iframe,iframe
+ endif
+ endif
+ write(11,*) '1 0'
endif
- write(11,*) '1 0'
+ ! dummy text for labels
+ write(11,*) '1 1'
+ write(11,*) 'a, b'
endif
-! dummy text for labels
- write(11,*) '1 1'
- write(11,*) 'a, b'
- endif
-! output data values
- mask_point = .false.
+ ! output data values
+ mask_point = .false.
-! output point data
- do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
-! four points for each element
- do ilocnum = 1,NGNOD2D_AVS_DX
- ibool_number = iglob(ilocnum+ieoff)
- if(.not. mask_point(ibool_number)) then
- if( .not. ALL_THRESHOLD_OFF ) then
- if(USE_OPENDX) then
- write(11,"(f7.2)") field_display(ilocnum+ieoff)
+ ! output point data
+ do ispec=1,nspectot_AVS_max
+ ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ! four points for each element
+ do ilocnum = 1,NGNOD2D_AVS_DX
+ ibool_number = iglob(ilocnum+ieoff)
+ if(.not. mask_point(ibool_number)) then
+ if( .not. ALL_THRESHOLD_OFF ) then
+ if(USE_OPENDX) then
+ write(11,"(f7.4)") field_display(ilocnum+ieoff)
+ else
+ write(11,"(i10,1x,f7.2)") ireorder(ibool_number),field_display(ilocnum+ieoff)
+ endif
else
- write(11,"(i10,1x,f7.2)") ireorder(ibool_number),field_display(ilocnum+ieoff)
+ if(USE_OPENDX) then
+ write(11,"(e7.4)") field_display(ilocnum+ieoff)
+ else
+ ! format spezifier has problems w/ very small values
+ !write(11,"(i10,1x,e7.4)") ireorder(ibool_number),field_display(ilocnum+ieoff)
+ write(11,*) ireorder(ibool_number),field_display(ilocnum+ieoff)
+ endif
endif
- else
- if(USE_OPENDX) then
- write(11,"(e7.4)") field_display(ilocnum+ieoff)
- else
- write(11,"(i10,1x,e7.4)") ireorder(ibool_number),field_display(ilocnum+ieoff)
- endif
endif
- endif
- mask_point(ibool_number) = .true.
+ mask_point(ibool_number) = .true.
+ enddo
enddo
- enddo
-! define OpenDX field
- if(USE_OPENDX) then
- write(11,*) 'attribute "dep" string "positions"'
- write(11,*) 'object "irregular positions irregular connections" class field'
- write(11,*) 'component "positions" value 1'
- write(11,*) 'component "connections" value 2'
- write(11,*) 'component "data" value 3'
- write(11,*) 'end'
+ ! define OpenDX field
+ if(USE_OPENDX) then
+ write(11,*) 'attribute "dep" string "positions"'
+ write(11,*) 'object "irregular positions irregular connections" class field'
+ write(11,*) 'component "positions" value 1'
+ write(11,*) 'component "connections" value 2'
+ write(11,*) 'component "data" value 3'
+ write(11,*) 'end'
+ endif
+
+ ! end of test for GMT format
endif
-! end of test for GMT format
- endif
+ if(.not. UNIQUE_FILE) close(11)
- if(.not. UNIQUE_FILE) close(11)
-
-! end of loop and test on all the time steps for all the movie images
- endif
+ ! end of loop and test on all the time steps for all the movie images
enddo
if(UNIQUE_FILE) close(11)
@@ -785,57 +836,39 @@
if(USE_GMT) print *,'GMT files are stored in ', trim(OUTPUT_FILES), '/gmt_*.xyz'
print *
- end subroutine create_movie_AVS_DX
+ end program xcreate_movie_AVS_DX
+
!
!=====================================================================
!
- subroutine read_AVS_DX_parameters(NEX_XI_out,NEX_ETA_out, &
- NSTEP_out,NTSTEP_BETWEEN_FRAMES_out, &
- NCHUNKS_out,MOVIE_SURFACE_out, &
- NPROCTOT_out,NEX_PER_PROC_XI_out,NEX_PER_PROC_ETA_out)
+ subroutine read_AVS_DX_parameters()
use constants
use shared_parameters
implicit none
-! parameters deduced from parameters read from file
- integer :: NPROCTOT_out,NCHUNKS_out
-
- integer :: NEX_PER_PROC_XI_out,NEX_PER_PROC_ETA_out
- integer :: NEX_XI_out,NEX_ETA_out
-
- logical :: MOVIE_SURFACE_out
- integer :: NSTEP_out,NTSTEP_BETWEEN_FRAMES_out
-
print *
print *,'reading parameter file'
print *
-! read the parameter file and compute additional parameters
+ ! read the parameter file and compute additional parameters
call read_compute_parameters()
-!
- if(MOVIE_COARSE) stop 'create_movie_AVS_DX does not work with MOVIE_COARSE'
- ! return values
- NPROCTOT_out = NPROCTOT
- NCHUNKS_out = NCHUNKS
+ ! checks
+ if(.not. MOVIE_SURFACE) stop 'movie frames were not saved by the solver'
- NEX_PER_PROC_XI_out = NEX_PER_PROC_XI
- NEX_PER_PROC_ETA_out = NEX_PER_PROC_ETA
- NEX_XI_out = NEX_XI
- NEX_ETA_out = NEX_ETA
+! if(MOVIE_COARSE) stop 'create_movie_AVS_DX does not work with MOVIE_COARSE'
- MOVIE_SURFACE_out = MOVIE_SURFACE
- NSTEP_out = NSTEP
- NTSTEP_BETWEEN_FRAMES_out = NTSTEP_BETWEEN_FRAMES
-
end subroutine read_AVS_DX_parameters
-! ------------------------------------------------------------------
+!
+!=====================================================================
+!
+
subroutine get_global_AVS(nspec,xp,yp,zp,iglob,loc,ifseg,nglob,npointot)
! this routine MUST be in double precision to avoid sensitivity
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/create_movie_GMT_global.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/create_movie_GMT_global.f90 2013-08-31 00:33:04 UTC (rev 22750)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/create_movie_GMT_global.f90 2013-09-02 10:52:04 UTC (rev 22751)
@@ -64,36 +64,37 @@
!---------------------
- integer i,j,it
- integer it1,it2
- integer ispec
+ integer :: i,j,it
+ integer :: it1,it2
+ integer :: ispec
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: x,y,z,displn
- real(kind=CUSTOM_REAL) xcoord,ycoord,zcoord,rval,thetaval,phival
- real(kind=CUSTOM_REAL) RRval,rhoval
- real(kind=CUSTOM_REAL) displx,disply,displz
- real(kind=CUSTOM_REAL) normal_x,normal_y,normal_z
- real(kind=CUSTOM_REAL) thetahat_x,thetahat_y,thetahat_z
- real(kind=CUSTOM_REAL) phihat_x,phihat_y
+ real(kind=CUSTOM_REAL) :: xcoord,ycoord,zcoord,rval,thetaval,phival
+ real(kind=CUSTOM_REAL) :: RRval,rhoval
+ real(kind=CUSTOM_REAL) :: displx,disply,displz
+ real(kind=CUSTOM_REAL) :: normal_x,normal_y,normal_z
+ real(kind=CUSTOM_REAL) :: thetahat_x,thetahat_y,thetahat_z
+ real(kind=CUSTOM_REAL) :: phihat_x,phihat_y
! to average maxima over past few steps
- double precision min_field_current,max_field_current,max_absol,max_average
+ double precision :: min_field_current,max_field_current,max_absol,max_average
double precision,dimension(:),allocatable :: max_history
integer :: nmax_history,imax
- real disp,lat,long
- integer nframes,iframe,USE_COMPONENT
+ real :: disp,lat,long
+ integer :: nframes,iframe,USE_COMPONENT
- character(len=150) outputname
+ character(len=150) :: outputname
- integer iproc,ipoin
+ integer :: iproc,ipoin
! for sorting routine
- integer npointot,ilocnum,ielm,ieoff,ispecloc,NIT
+ integer :: npointot,ilocnum,ielm,ieoff,ispecloc
+ integer :: NIT
double precision, dimension(:), allocatable :: xp,yp,zp,field_display
! for dynamic memory allocation
- integer ierror
+ integer :: ierror
! movie files stored by solver
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
@@ -104,11 +105,13 @@
real(kind=CUSTOM_REAL) :: LAT_SOURCE,LON_SOURCE,DEP_SOURCE
real(kind=CUSTOM_REAL) :: dist_lon,dist_lat,distance,mute_factor,val
- character(len=256) line
real(kind=CUSTOM_REAL) :: cmt_hdur,cmt_t_shift,t0,hdur
real(kind=CUSTOM_REAL) :: xmesh,ymesh,zmesh
integer :: istamp1,istamp2
+ character(len=256) :: line
+ character(len=150) :: CMTSOLUTION
+
! ************** PROGRAM STARTS HERE **************
print *
@@ -124,9 +127,35 @@
if(.not. MOVIE_SURFACE) stop 'movie frames were not saved by the solver'
+ ! user input
+ print *,'movie frames have been saved every ',NTSTEP_BETWEEN_FRAMES,' time steps'
print *
+
+ print *,'enter first time step of movie (e.g. 1)'
+ read(5,*) it1
+
+ print *,'enter last time step of movie (e.g. ',NSTEP,'or -1 for all)'
+ 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'
+
+ print *,'enter output ascii (F) or binary (T)'
+ read(5,*) OUTPUT_BINARY
+
+ print *
+ print *,'--------'
+ print *
+
+ ! checks options
+ if( it2 == -1 ) it2 = NSTEP
+
+
+ print *
print *,'There are ',NPROCTOT,' slices numbered from 0 to ',NPROCTOT-1
print *
+
if(MOVIE_COARSE) then
! note:
! nex_per_proc_xi*nex_per_proc_eta = nex_xi*nex_eta/nproc = nspec2d_top(iregion_crust_mantle)
@@ -138,9 +167,10 @@
ilocnum = NGLLX*NGLLY*NEX_PER_PROC_XI*NEX_PER_PROC_ETA
NIT = 1
endif
+
print *
print *,'Allocating arrays for reading data of size ',ilocnum*NPROCTOT,'=', &
- 6*ilocnum*NPROCTOT*CUSTOM_REAL/1000000,'MB'
+ 6*ilocnum*NPROCTOT*CUSTOM_REAL/1024./1024.,'MB'
print *
! allocates movie arrays
@@ -174,28 +204,7 @@
allocate(displn(NGLLX,NGLLY),stat=ierror)
if(ierror /= 0) stop 'error while allocating displn'
- print *,'movie frames have been saved every ',NTSTEP_BETWEEN_FRAMES,' time steps'
- print *
- ! user input
- print *,'--------'
- print *,'enter first time step of movie (e.g. 1)'
- read(5,*) it1
-
- print *,'enter last time step of movie (e.g. ',NSTEP,'or -1 for all)'
- 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'
-
- print *,'enter output ascii (F) or binary (T)'
- read(5,*) OUTPUT_BINARY
- print *,'--------'
-
- ! checks options
- if( it2 == -1 ) it2 = NSTEP
-
print *
print *,'looping from ',it1,' to ',it2,' every ',NTSTEP_BETWEEN_FRAMES,' time steps'
@@ -221,8 +230,8 @@
print *
- print *,'Allocating 4 outputdata arrays of size 4*CUSTOM_REAL',npointot,'=', &
- 4*npointot*CUSTOM_REAL/1000000,' MB'
+ print *,'Allocating 4 outputdata arrays of size 4 * CUSTOM_REAL *',npointot,'=', &
+ 4*npointot*CUSTOM_REAL/1024.0/1024.0,' MB'
print *
allocate(xp(npointot),stat=ierror)
@@ -264,28 +273,30 @@
cmt_t_shift = 0.0
cmt_hdur = 0.0
+ call get_value_string(CMTSOLUTION, 'solver.CMTSOLUTION','DATA/CMTSOLUTION')
+
! reads in source lat/lon
- open(22,file="DATA/CMTSOLUTION",status='old',action='read',iostat=ierror )
+ open(IIN,file=trim(CMTSOLUTION),status='old',action='read',iostat=ierror )
if( ierror == 0 ) then
! skip first line, event name,timeshift,half duration
- read(22,*,iostat=ierror ) line ! PDE line
- read(22,*,iostat=ierror ) line ! event name
+ read(IIN,*,iostat=ierror ) line ! PDE line
+ read(IIN,*,iostat=ierror ) line ! event name
! timeshift
- read(22,'(a256)',iostat=ierror ) line
+ read(IIN,'(a256)',iostat=ierror ) line
if( ierror == 0 ) read(line(12:len_trim(line)),*) cmt_t_shift
! halfduration
- read(22,'(a256)',iostat=ierror ) line
+ read(IIN,'(a256)',iostat=ierror ) line
if( ierror == 0 ) read(line(15:len_trim(line)),*) cmt_hdur
! latitude
- read(22,'(a256)',iostat=ierror ) line
+ read(IIN,'(a256)',iostat=ierror ) line
if( ierror == 0 ) read(line(10:len_trim(line)),*) LAT_SOURCE
! longitude
- read(22,'(a256)',iostat=ierror ) line
+ read(IIN,'(a256)',iostat=ierror ) line
if( ierror == 0 ) read(line(11:len_trim(line)),*) LON_SOURCE
! depth
- read(22,'(a256)',iostat=ierror ) line
+ read(IIN,'(a256)',iostat=ierror ) line
if( ierror == 0 ) read(line(7:len_trim(line)),*) DEP_SOURCE
- close(22)
+ close(IIN)
endif
! effective half duration in movie runs
hdur = sqrt( cmt_hdur**2 + HDUR_MOVIE**2)
@@ -328,537 +339,567 @@
!--- ****** read data saved by solver ******
+ ! movie point locations
+ outputname = "/moviedata_xyz.bin"
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//trim(outputname), &
+ status='old',action='read',form='unformatted',iostat=ierror)
+ if(ierror /= 0 ) then
+ print*,'error opening file: ',trim(OUTPUT_FILES)//trim(outputname)
+ stop 'error opening moviedata file'
+ endif
+
+ ! 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
+ close(IOUT)
+
! --------------------------------------
+
istamp1 = 0
istamp2 = 0
iframe = 0
! loop on all the time steps in the range entered
do it = it1,it2
- ! check if time step corresponds to a movie frame
- if(mod(it,NTSTEP_BETWEEN_FRAMES) == 0) then
- iframe = iframe + 1
+ ! check if time step corresponds to a movie frame
+ if(mod(it,NTSTEP_BETWEEN_FRAMES) /= 0) cycle
- ! read all the elements from the same file
- write(outputname,"('OUTPUT_FILES/moviedata',i6.6)") it
- open(unit=IOUT,file=outputname,status='old',form='unformatted')
+ iframe = iframe + 1
- print *
- print *,'reading snapshot time step ',it,' out of ',NSTEP,' file ',outputname
- !print *
+ print *
+ print *,'reading snapshot time step ',it,' out of ',NSTEP
+ print *
- ! 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
+ ! read all the elements from the same file
+ write(outputname,"('/moviedata',i6.6)") it
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//trim(outputname), &
+ status='old',action='read',form='unformatted',iostat=ierror)
+ if(ierror /= 0 ) then
+ print*,'error opening file: ',trim(OUTPUT_FILES)//trim(outputname)
+ stop 'error opening moviedata file'
+ endif
- ! reads in associated values (displacement or velocity..)
- read(IOUT) store_val_ux
- read(IOUT) store_val_uy
- read(IOUT) store_val_uz
+ ! reads in associated values (displacement or velocity..)
+ read(IOUT) store_val_ux
+ read(IOUT) store_val_uy
+ read(IOUT) store_val_uz
- close(IOUT)
- !print *, 'finished reading ',outputname
+ close(IOUT)
- ! mutes source region
- if( MUTE_SOURCE ) then
- ! initialize factor
- mute_factor = 1.0
+ print *,' file ',trim(OUTPUT_FILES)//trim(outputname)
+ print *
+ !debug
+ print *," data ux min/max: ",minval(store_val_ux),maxval(store_val_ux)
+ print *," data uy min/max: ",minval(store_val_uy),maxval(store_val_uy)
+ print *," data uz min/max: ",minval(store_val_uz),maxval(store_val_uz)
+ print *
- print*,'simulation time: ',(it-1)*DT - t0,'(s)'
+ ! mutes source region
+ if( MUTE_SOURCE ) then
+ ! initialize factor
+ mute_factor = 1.0
- ! muting radius grows/shrinks with time
- if( (it-1)*DT - t0 > STARTTIME_TO_MUTE ) then
+ print*,'simulation time: ',(it-1)*DT - t0,'(s)'
- ! approximate wavefront travel distance in degrees
- ! (~3.5 km/s wave speed for surface waves)
- distance = 3.5 * ((it-1)*DT-t0) / 6371.0 * RADIANS_TO_DEGREES
+ ! muting radius grows/shrinks with time
+ if( (it-1)*DT - t0 > STARTTIME_TO_MUTE ) then
- ! approximate distance to source (in degrees)
- ! (shrinks if waves travel back from antipode)
- !do while ( distance > 360. )
- ! distance = distance - 360.
- !enddo
- ! waves are back at origin, no source tapering anymore
- if( distance > 360.0 ) distance = 0.0
- ! shrinks when waves reached antipode
- !if( distance > 180. ) distance = 360. - distance
- ! shrinks when waves reached half-way to antipode
- if( distance > 90.0 ) distance = 90.0 - distance
+ ! approximate wavefront travel distance in degrees
+ ! (~3.5 km/s wave speed for surface waves)
+ distance = 3.5 * ((it-1)*DT-t0) / 6371.0 * RADIANS_TO_DEGREES
- ! limit size around source (in degrees)
- if( distance < 0.0 ) distance = 0.0
- if( distance > 80.0 ) distance = 80.0
+ ! approximate distance to source (in degrees)
+ ! (shrinks if waves travel back from antipode)
+ !do while ( distance > 360. )
+ ! distance = distance - 360.
+ !enddo
+ ! waves are back at origin, no source tapering anymore
+ if( distance > 360.0 ) distance = 0.0
+ ! shrinks when waves reached antipode
+ !if( distance > 180. ) distance = 360. - distance
+ ! shrinks when waves reached half-way to antipode
+ if( distance > 90.0 ) distance = 90.0 - distance
- print*,'muting radius: ',0.7 * distance,'(degrees)'
+ ! limit size around source (in degrees)
+ if( distance < 0.0 ) distance = 0.0
+ if( distance > 80.0 ) distance = 80.0
- ! new radius of mute area (in rad)
- RADIUS_TO_MUTE = 0.7 * distance * DEGREES_TO_RADIANS
- else
- ! mute_factor used at the beginning for scaling displacement values
- if( STARTTIME_TO_MUTE > TINYVAL ) then
- ! mute factor 1: no masking out
- ! 0: masks out values (within mute radius)
- ! linear scaling between [0,1]:
- ! from 0 (simulation time equal to zero )
- ! to 1 (simulation time equals starttime_to_mute)
- mute_factor = 1.0 - ( STARTTIME_TO_MUTE - ((it-1)*DT-t0) ) / (STARTTIME_TO_MUTE+t0)
- ! threshold value for mute_factor
- if( mute_factor < TINYVAL ) mute_factor = TINYVAL
- if( mute_factor > 1.0 ) mute_factor = 1.0
- endif
- endif
+ print*,'muting radius: ',0.7 * distance,'(degrees)'
+
+ ! new radius of mute area (in rad)
+ RADIUS_TO_MUTE = 0.7 * distance * DEGREES_TO_RADIANS
+ else
+ ! mute_factor used at the beginning for scaling displacement values
+ if( STARTTIME_TO_MUTE > TINYVAL ) then
+ ! mute factor 1: no masking out
+ ! 0: masks out values (within mute radius)
+ ! linear scaling between [0,1]:
+ ! from 0 (simulation time equal to zero )
+ ! to 1 (simulation time equals starttime_to_mute)
+ mute_factor = 1.0 - ( STARTTIME_TO_MUTE - ((it-1)*DT-t0) ) / (STARTTIME_TO_MUTE+t0)
+ ! threshold value for mute_factor
+ if( mute_factor < TINYVAL ) mute_factor = TINYVAL
+ if( mute_factor > 1.0 ) mute_factor = 1.0
endif
+ endif
+ endif
- ! clear number of elements kept
- ispec = 0
+ ! clear number of elements kept
+ ispec = 0
- ! read points for all the slices
- print *,'Converting to geo-coordinates'
- do iproc = 0,NPROCTOT-1
- ! reset point number
- ipoin = 0
- do ispecloc = 1,NEX_PER_PROC_XI*NEX_PER_PROC_ETA
- do j = 1,NGLLY,NIT
- do i = 1,NGLLX,NIT
- ipoin = ipoin + 1
+ print *,'Converting to geo-coordinates'
- ! coordinates actually contain r theta phi
- xcoord = store_val_x(ipoin,iproc)
- ycoord = store_val_y(ipoin,iproc)
- zcoord = store_val_z(ipoin,iproc)
+ ! read points for all the slices
+ do iproc = 0,NPROCTOT-1
- displx = store_val_ux(ipoin,iproc)
- disply = store_val_uy(ipoin,iproc)
- displz = store_val_uz(ipoin,iproc)
+ ! reset point number
+ ipoin = 0
- ! 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)
+ do ispecloc = 1,NEX_PER_PROC_XI*NEX_PER_PROC_ETA
- ! save the results for this element
- x(i,j) = xcoord
- y(i,j) = ycoord
- z(i,j) = zcoord
+ ! gets element values
+ do j = 1,NGLLY,NIT
+ do i = 1,NGLLX,NIT
+ ipoin = ipoin + 1
- ! 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
+ ! coordinates actually contain r theta phi
+ xcoord = store_val_x(ipoin,iproc)
+ ycoord = store_val_y(ipoin,iproc)
+ zcoord = store_val_z(ipoin,iproc)
- displn(i,j) = displx*normal_x + disply*normal_y + displz*normal_z
+ displx = store_val_ux(ipoin,iproc)
+ disply = store_val_uy(ipoin,iproc)
+ displz = store_val_uz(ipoin,iproc)
- else if(USE_COMPONENT == 2) then
+ ! 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 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
+ ! save the results for this element
+ x(i,j) = xcoord
+ y(i,j) = ycoord
+ z(i,j) = zcoord
- displn(i,j) = - (displx*thetahat_x + disply*thetahat_y + displz*thetahat_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
- else if(USE_COMPONENT == 3) then
+ displn(i,j) = displx*normal_x + disply*normal_y + displz*normal_z
- ! 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
+ else if(USE_COMPONENT == 2) then
- displn(i,j) = displx*phihat_x + disply*phihat_y
- endif
+ ! 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)
- ! mute values
- if( MUTE_SOURCE ) then
+ else if(USE_COMPONENT == 3) then
- ! distance in colatitude (in rad)
- ! note: this mixes geocentric (point location) and geographic (source location) coordinates;
- ! since we only need approximate distances here,
- ! this should be fine for the muting region
- dist_lat = thetaval - LAT_SOURCE
+ ! 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
- ! distance in longitude (in rad)
- ! checks source longitude range
- if( LON_SOURCE - RADIUS_TO_MUTE < -PI .or. LON_SOURCE + RADIUS_TO_MUTE > PI ) then
- ! source close to 180. longitudes, shifts range to [0, 2PI]
- if( phival < 0.0 ) phival = phival + TWO_PI
- if( LON_SOURCE < 0.0 ) then
- dist_lon = phival - (LON_SOURCE + TWO_PI)
- else
- dist_lon = phival - LON_SOURCE
- endif
- else
- ! source well between range to [-PI, PI]
- ! shifts phival to be like LON_SOURCE between [-PI,PI]
- if( phival > PI ) phival = phival - TWO_PI
- if( phival < -PI ) phival = phival + TWO_PI
+ displn(i,j) = displx*phihat_x + disply*phihat_y
+ endif
- dist_lon = phival - LON_SOURCE
- endif
- ! distance of point to source (in rad)
- distance = sqrt(dist_lat**2 + dist_lon**2)
+ ! mute values
+ if( MUTE_SOURCE ) then
- ! mutes source region values
- if ( distance < RADIUS_TO_MUTE ) then
- ! muting takes account of the event time
- if( (it-1)*DT-t0 > STARTTIME_TO_MUTE ) then
- ! wavefield will be tapered to mask out noise in source area
- ! factor from 0 to 1
- mute_factor = ( 0.5*(1.0 - cos(distance/RADIUS_TO_MUTE*PI)) )**6
- ! factor from 0.01 to 1
- mute_factor = mute_factor * 0.99 + 0.01
- displn(i,j) = displn(i,j) * mute_factor
- else
- ! wavefield will initially be scaled down to avoid noise being amplified at beginning
- displn(i,j) = displn(i,j) * mute_factor
- endif
- endif
+ ! distance in colatitude (in rad)
+ ! note: this mixes geocentric (point location) and geographic (source location) coordinates;
+ ! since we only need approximate distances here,
+ ! this should be fine for the muting region
+ dist_lat = thetaval - LAT_SOURCE
+ ! distance in longitude (in rad)
+ ! checks source longitude range
+ if( LON_SOURCE - RADIUS_TO_MUTE < -PI .or. LON_SOURCE + RADIUS_TO_MUTE > PI ) then
+ ! source close to 180. longitudes, shifts range to [0, 2PI]
+ if( phival < 0.0 ) phival = phival + TWO_PI
+ if( LON_SOURCE < 0.0 ) then
+ dist_lon = phival - (LON_SOURCE + TWO_PI)
+ else
+ dist_lon = phival - LON_SOURCE
endif
+ else
+ ! source well between range to [-PI, PI]
+ ! shifts phival to be like LON_SOURCE between [-PI,PI]
+ if( phival > PI ) phival = phival - TWO_PI
+ if( phival < -PI ) phival = phival + TWO_PI
- enddo !i
- enddo !j
+ dist_lon = phival - LON_SOURCE
+ endif
+ ! distance of point to source (in rad)
+ distance = sqrt(dist_lat**2 + dist_lon**2)
- ispec = ispec + 1
- if(MOVIE_COARSE) then
- ielm = ispec-1
- else
- ielm = (NGLLX-1)*(NGLLY-1)*(ispec-1)
- endif
- do j = 1,NGLLY-NIT
- do i = 1,NGLLX-NIT
- if(MOVIE_COARSE) then
- ieoff = ielm+1
+ ! mutes source region values
+ if ( distance < RADIUS_TO_MUTE ) then
+ ! muting takes account of the event time
+ if( (it-1)*DT-t0 > STARTTIME_TO_MUTE ) then
+ ! wavefield will be tapered to mask out noise in source area
+ ! factor from 0 to 1
+ mute_factor = ( 0.5*(1.0 - cos(distance/RADIUS_TO_MUTE*PI)) )**6
+ ! factor from 0.01 to 1
+ mute_factor = mute_factor * 0.99 + 0.01
+ displn(i,j) = displn(i,j) * mute_factor
else
- ieoff = (ielm+(i-1)+(j-1)*(NGLLX-1))+1
+ ! wavefield will initially be scaled down to avoid noise being amplified at beginning
+ displn(i,j) = displn(i,j) * mute_factor
endif
+ endif
-! for movie_coarse e.g. x(i,j) is defined at x(1,1), x(1,NGLLY), x(NGLLX,1) and x(NGLLX,NGLLY)
-! be aware that for the cubed sphere, the mapping changes for different chunks,
-! i.e. e.g. x(1,1) and x(5,5) flip left and right sides of the elements in geographical coordinates
- if(MOVIE_COARSE) then
- if(NCHUNKS == 6) then
- ! chunks mapped such that element corners increase in long/lat
- select case (iproc/NPROC+1)
- case(CHUNK_AB)
- xp(ieoff) = dble(x(1,NGLLY))
- yp(ieoff) = dble(y(1,NGLLY))
- zp(ieoff) = dble(z(1,NGLLY))
- field_display(ieoff) = dble(displn(1,NGLLY))
- case(CHUNK_AB_ANTIPODE)
- xp(ieoff) = dble(x(1,1))
- yp(ieoff) = dble(y(1,1))
- zp(ieoff) = dble(z(1,1))
- field_display(ieoff) = dble(displn(1,1))
- case(CHUNK_AC)
- xp(ieoff) = dble(x(1,NGLLY))
- yp(ieoff) = dble(y(1,NGLLY))
- zp(ieoff) = dble(z(1,NGLLY))
- field_display(ieoff) = dble(displn(1,NGLLY))
- case(CHUNK_AC_ANTIPODE)
- xp(ieoff) = dble(x(1,1))
- yp(ieoff) = dble(y(1,1))
- zp(ieoff) = dble(z(1,1))
- field_display(ieoff) = dble(displn(1,1))
- case(CHUNK_BC)
- xp(ieoff) = dble(x(1,NGLLY))
- yp(ieoff) = dble(y(1,NGLLY))
- zp(ieoff) = dble(z(1,NGLLY))
- field_display(ieoff) = dble(displn(1,NGLLY))
- case(CHUNK_BC_ANTIPODE)
- xp(ieoff) = dble(x(NGLLX,NGLLY))
- yp(ieoff) = dble(y(NGLLX,NGLLY))
- zp(ieoff) = dble(z(NGLLX,NGLLY))
- field_display(ieoff) = dble(displn(NGLLX,NGLLY))
- case default
- stop 'incorrect chunk number'
- end select
- else
+ endif
+
+ enddo !i
+ enddo !j
+
+ ispec = ispec + 1
+ if(MOVIE_COARSE) then
+ ielm = ispec-1
+ else
+ ielm = (NGLLX-1)*(NGLLY-1)*(ispec-1)
+ endif
+
+ do j = 1,NGLLY-NIT
+ do i = 1,NGLLX-NIT
+ ! offset (starts at 1)
+ if(MOVIE_COARSE) then
+ ieoff = ielm+1
+ else
+ ieoff = (ielm+(i-1)+(j-1)*(NGLLX-1))+1
+ endif
+
+ ! for movie_coarse e.g. x(i,j) is defined at x(1,1), x(1,NGLLY), x(NGLLX,1) and x(NGLLX,NGLLY)
+ ! be aware that for the cubed sphere, the mapping changes for different chunks,
+ ! i.e. e.g. x(1,1) and x(5,5) flip left and right sides of the elements in geographical coordinates
+ if(MOVIE_COARSE) then
+ if(NCHUNKS == 6) then
+ ! chunks mapped such that element corners increase in long/lat
+ select case (iproc/NPROC+1)
+ case(CHUNK_AB)
+ xp(ieoff) = dble(x(1,NGLLY))
+ yp(ieoff) = dble(y(1,NGLLY))
+ zp(ieoff) = dble(z(1,NGLLY))
+ field_display(ieoff) = dble(displn(1,NGLLY))
+ case(CHUNK_AB_ANTIPODE)
xp(ieoff) = dble(x(1,1))
yp(ieoff) = dble(y(1,1))
zp(ieoff) = dble(z(1,1))
field_display(ieoff) = dble(displn(1,1))
- endif ! NCHUNKS
- else
- xp(ieoff) = dble(x(i,j))
- yp(ieoff) = dble(y(i,j))
- zp(ieoff) = dble(z(i,j))
- field_display(ieoff) = dble(displn(i,j))
- endif ! MOVIE_COARSE
+ case(CHUNK_AC)
+ xp(ieoff) = dble(x(1,NGLLY))
+ yp(ieoff) = dble(y(1,NGLLY))
+ zp(ieoff) = dble(z(1,NGLLY))
+ field_display(ieoff) = dble(displn(1,NGLLY))
+ case(CHUNK_AC_ANTIPODE)
+ xp(ieoff) = dble(x(1,1))
+ yp(ieoff) = dble(y(1,1))
+ zp(ieoff) = dble(z(1,1))
+ field_display(ieoff) = dble(displn(1,1))
+ case(CHUNK_BC)
+ xp(ieoff) = dble(x(1,NGLLY))
+ yp(ieoff) = dble(y(1,NGLLY))
+ zp(ieoff) = dble(z(1,NGLLY))
+ field_display(ieoff) = dble(displn(1,NGLLY))
+ case(CHUNK_BC_ANTIPODE)
+ xp(ieoff) = dble(x(NGLLX,NGLLY))
+ yp(ieoff) = dble(y(NGLLX,NGLLY))
+ zp(ieoff) = dble(z(NGLLX,NGLLY))
+ field_display(ieoff) = dble(displn(NGLLX,NGLLY))
+ case default
+ stop 'incorrect chunk number'
+ end select
+ else
+ ! takes lower-left point only
+ xp(ieoff) = dble(x(1,1))
+ yp(ieoff) = dble(y(1,1))
+ zp(ieoff) = dble(z(1,1))
+ field_display(ieoff) = dble(displn(1,1))
+ endif ! NCHUNKS
+ else
+ xp(ieoff) = dble(x(i,j))
+ yp(ieoff) = dble(y(i,j))
+ zp(ieoff) = dble(z(i,j))
+ field_display(ieoff) = dble(displn(i,j))
+ endif ! MOVIE_COARSE
- ! determines noth / sourth pole index for stamping maximum values
- if( USE_AVERAGED_MAXIMUM .and. NORMALIZE_VALUES ) then
- xmesh = xp(ieoff)
- ymesh = yp(ieoff)
- zmesh = zp(ieoff)
- if(zmesh > -SMALL_VAL_ANGLE .and. zmesh <= ZERO) zmesh = -SMALL_VAL_ANGLE
- if(zmesh < SMALL_VAL_ANGLE .and. zmesh >= ZERO) zmesh = SMALL_VAL_ANGLE
- thetaval = atan2(sqrt(xmesh*xmesh+ymesh*ymesh),zmesh)
- ! thetaval between 0 and PI / 2
- !print*,'thetaval:',thetaval * 180. / PI
- ! close to north pole
- if( thetaval >= 0.495 * PI ) istamp1 = ieoff
- ! close to south pole
- if( thetaval <= 0.01 ) istamp2 = ieoff
- endif
+ ! determines noth / sourth pole index for stamping maximum values
+ if( USE_AVERAGED_MAXIMUM .and. NORMALIZE_VALUES ) then
+ xmesh = xp(ieoff)
+ ymesh = yp(ieoff)
+ zmesh = zp(ieoff)
+ if(zmesh > -SMALL_VAL_ANGLE .and. zmesh <= ZERO) zmesh = -SMALL_VAL_ANGLE
+ if(zmesh < SMALL_VAL_ANGLE .and. zmesh >= ZERO) zmesh = SMALL_VAL_ANGLE
+ thetaval = atan2(sqrt(xmesh*xmesh+ymesh*ymesh),zmesh)
+ ! thetaval between 0 and PI / 2
+ !print*,'thetaval:',thetaval * 180. / PI
+ ! close to north pole
+ if( thetaval >= 0.495 * PI ) istamp1 = ieoff
+ ! close to south pole
+ if( thetaval <= 0.01 ) istamp2 = ieoff
+ endif
- enddo !i
- enddo !j
+ enddo !i
+ enddo !j
- enddo !ispec
+ enddo !ispec
- enddo !nproc
+ enddo !nproc
- ! compute min and max of data value to normalize
- min_field_current = minval(field_display(:))
- max_field_current = maxval(field_display(:))
+ ! compute min and max of data value to normalize
+ min_field_current = minval(field_display(:))
+ max_field_current = maxval(field_display(:))
- ! print minimum and maximum amplitude in current snapshot
- print *
- print *,'minimum amplitude in current snapshot = ',min_field_current
- print *,'maximum amplitude in current snapshot = ',max_field_current
+ ! print minimum and maximum amplitude in current snapshot
+ print *
+ print *,'minimum amplitude in current snapshot = ',min_field_current
+ print *,'maximum amplitude in current snapshot = ',max_field_current
- ! takes average over last few snapshots available and uses it
- ! to normalize field values
- if( USE_AVERAGED_MAXIMUM ) then
+ ! takes average over last few snapshots available and uses it
+ ! to normalize field values
+ if( USE_AVERAGED_MAXIMUM ) then
- ! (average) maximum between positive and negative values
- max_absol = (abs(min_field_current)+abs(max_field_current))/2.0
+ ! (average) maximum between positive and negative values
+ max_absol = (abs(min_field_current)+abs(max_field_current))/2.0
- ! stores last few maxima
- ! index between 1 and nmax_history
- imax = mod(iframe-1,nmax_history) + 1
- max_history( imax ) = max_absol
+ ! stores last few maxima
+ ! index between 1 and nmax_history
+ imax = mod(iframe-1,nmax_history) + 1
+ max_history( imax ) = max_absol
- ! average over history
- max_average = sum( max_history )
- if( iframe < nmax_history ) then
- ! history not filled yet, only average over available entries
- max_average = max_average / iframe
- else
- ! average over all history entries
- max_average = max_average / nmax_history
- endif
+ ! average over history
+ max_average = sum( max_history )
+ if( iframe < nmax_history ) then
+ ! history not filled yet, only average over available entries
+ max_average = max_average / iframe
+ else
+ ! average over all history entries
+ max_average = max_average / nmax_history
+ endif
- print *,'maximum amplitude averaged in current snapshot = ',max_absol
- print *,'maximum amplitude over averaged last snapshots = ',max_average
+ print *,'maximum amplitude averaged in current snapshot = ',max_absol
+ print *,'maximum amplitude over averaged last snapshots = ',max_average
- ! thresholds positive & negative maximum values
- where( field_display(:) > max_absol ) field_display = max_absol
- where( field_display(:) < - max_absol ) field_display = -max_absol
+ ! thresholds positive & negative maximum values
+ where( field_display(:) > max_absol ) field_display = max_absol
+ where( field_display(:) < - max_absol ) field_display = -max_absol
- ! sets new maxima for decaying wavefield
- ! this should avoid flickering when normalizing wavefields
- if( NORMALIZE_VALUES ) then
- ! checks stamp indices for maximum values
- if( istamp1 == 0 ) istamp1 = ieoff
- if( istamp2 == 0 ) istamp2 = ieoff-1
- !print *, 'stamp: ',istamp1,istamp2
+ ! sets new maxima for decaying wavefield
+ ! this should avoid flickering when normalizing wavefields
+ if( NORMALIZE_VALUES ) then
+ ! checks stamp indices for maximum values
+ if( istamp1 == 0 ) istamp1 = ieoff
+ if( istamp2 == 0 ) istamp2 = ieoff-1
+ !print *, 'stamp: ',istamp1,istamp2
- if( max_absol < max_average ) then
- ! distance (in degree) of surface waves travelled
- distance = 3.5 * ((it-1)*DT-t0) / 6371.0 * RADIANS_TO_DEGREES
- if( distance > 10.0 .and. distance <= 20.0 ) then
- ! smooth transition between 10 and 20 degrees
- ! sets positive and negative maximum
- field_display(istamp1) = + max_absol + (max_average-max_absol) * (distance - 10.0)/10.0
- field_display(istamp2) = - ( max_absol + (max_average-max_absol) * (distance - 10.0)/10.0 )
- else if( distance > 20.0 ) then
- ! sets positive and negative maximum
- field_display(istamp1) = + max_average
- field_display(istamp2) = - max_average
- endif
- else
- ! thresholds positive & negative maximum values
- where( field_display(:) > max_average ) field_display = max_average
- where( field_display(:) < - max_average ) field_display = -max_average
- ! sets positive and negative maximum
- field_display(istamp1) = + max_average
- field_display(istamp2) = - max_average
- endif
- ! updates current wavefield maxima
- min_field_current = minval(field_display(:))
- max_field_current = maxval(field_display(:))
- max_absol = (abs(min_field_current)+abs(max_field_current))/2.0
+ if( max_absol < max_average ) then
+ ! distance (in degree) of surface waves travelled
+ distance = 3.5 * ((it-1)*DT-t0) / 6371.0 * RADIANS_TO_DEGREES
+ if( distance > 10.0 .and. distance <= 20.0 ) then
+ ! smooth transition between 10 and 20 degrees
+ ! sets positive and negative maximum
+ field_display(istamp1) = + max_absol + (max_average-max_absol) * (distance - 10.0)/10.0
+ field_display(istamp2) = - ( max_absol + (max_average-max_absol) * (distance - 10.0)/10.0 )
+ else if( distance > 20.0 ) then
+ ! sets positive and negative maximum
+ field_display(istamp1) = + max_average
+ field_display(istamp2) = - max_average
endif
-
- ! scales field values up to match average
- if( abs(max_absol) > TINYVAL) &
- field_display = field_display * max_average / max_absol
-
- ! thresholds after scaling positive & negative maximum values
+ else
+ ! thresholds positive & negative maximum values
where( field_display(:) > max_average ) field_display = max_average
where( field_display(:) < - max_average ) field_display = -max_average
+ ! sets positive and negative maximum
+ field_display(istamp1) = + max_average
+ field_display(istamp2) = - max_average
+ endif
+ ! updates current wavefield maxima
+ min_field_current = minval(field_display(:))
+ max_field_current = maxval(field_display(:))
+ max_absol = (abs(min_field_current)+abs(max_field_current))/2.0
+ endif
- ! normalizes field values
- if( NORMALIZE_VALUES ) then
- if( MUTE_SOURCE ) then
- ! checks if source wavefield kicked in
- if( (it-1)*DT - t0 > STARTTIME_TO_MUTE ) then
- ! wavefield should be visible at surface now
- ! normalizes wavefield
- if( abs(max_average) > TINYVAL ) field_display = field_display / max_average
- else
- ! no wavefield yet assumed
+ ! scales field values up to match average
+ if( abs(max_absol) > TINYVAL) &
+ field_display = field_display * max_average / max_absol
- ! we set two single field values (last in array)
- ! to: +/- 100 * max_average
- ! to avoid further amplifying when
- ! a normalization routine is used for rendering images
- ! (which for example is the case for shakemovies)
- if( STARTTIME_TO_MUTE > TINYVAL ) then
- ! with additional scale factor:
- ! linear scaling between [0,1]:
- ! from 0 (simulation time equal to -t0 )
- ! to 1 (simulation time equals starttime_to_mute)
- mute_factor = 1.0 - ( STARTTIME_TO_MUTE - ((it-1)*DT-t0) ) / (STARTTIME_TO_MUTE+t0)
- ! takes complement and shifts scale to (1,100)
- ! thus, mute factor is 100 at simulation start and 1.0 at starttime_to_mute
- mute_factor = abs(1.0 - mute_factor) * 99.0 + 1.0
- ! positive and negative maximum reach average when wavefield appears
- val = mute_factor * max_average
- else
- ! uses a constant factor
- val = 100.0 * max_average
- endif
- ! positive and negative maximum
- field_display(istamp1) = + val
- field_display(istamp2) = - val
- if( abs(max_average) > TINYVAL ) field_display = field_display / val
- endif
+ ! thresholds after scaling positive & negative maximum values
+ where( field_display(:) > max_average ) field_display = max_average
+ where( field_display(:) < - max_average ) field_display = -max_average
+
+ ! normalizes field values
+ if( NORMALIZE_VALUES ) then
+ if( MUTE_SOURCE ) then
+ ! checks if source wavefield kicked in
+ if( (it-1)*DT - t0 > STARTTIME_TO_MUTE ) then
+ ! wavefield should be visible at surface now
+ ! normalizes wavefield
+ if( abs(max_average) > TINYVAL ) field_display = field_display / max_average
+ else
+ ! no wavefield yet assumed
+
+ ! we set two single field values (last in array)
+ ! to: +/- 100 * max_average
+ ! to avoid further amplifying when
+ ! a normalization routine is used for rendering images
+ ! (which for example is the case for shakemovies)
+ if( STARTTIME_TO_MUTE > TINYVAL ) then
+ ! with additional scale factor:
+ ! linear scaling between [0,1]:
+ ! from 0 (simulation time equal to -t0 )
+ ! to 1 (simulation time equals starttime_to_mute)
+ mute_factor = 1.0 - ( STARTTIME_TO_MUTE - ((it-1)*DT-t0) ) / (STARTTIME_TO_MUTE+t0)
+ ! takes complement and shifts scale to (1,100)
+ ! thus, mute factor is 100 at simulation start and 1.0 at starttime_to_mute
+ mute_factor = abs(1.0 - mute_factor) * 99.0 + 1.0
+ ! positive and negative maximum reach average when wavefield appears
+ val = mute_factor * max_average
else
- ! no source to mute
- ! normalizes wavefield
- if( abs(max_average) > TINYVAL ) field_display = field_display / max_average
+ ! uses a constant factor
+ val = 100.0 * max_average
endif
+ ! positive and negative maximum
+ field_display(istamp1) = + val
+ field_display(istamp2) = - val
+ if( abs(max_average) > TINYVAL ) field_display = field_display / val
endif
+ else
+ ! no source to mute
+ ! normalizes wavefield
+ if( abs(max_average) > TINYVAL ) field_display = field_display / max_average
endif
+ endif
+ endif
- print *
- print *,'initial number of points (with multiples) was ',npointot
- print *,'final number of points is ',ieoff
+ print *
+ print *,'initial number of points (with multiples) was ',npointot
+ print *,'final number of points is ',ieoff
- !--- ****** create GMT file ******
+ !--- ****** create GMT file ******
- ! create file name and open file
- if(OUTPUT_BINARY) then
- if(USE_COMPONENT == 1) then
- write(outputname,"('bin_movie_',i6.6,'.d')") it
- else if(USE_COMPONENT == 2) then
- write(outputname,"('bin_movie_',i6.6,'.N')") it
- else if(USE_COMPONENT == 3) then
- write(outputname,"('bin_movie_',i6.6,'.E')") it
- endif
- open(unit=11,file='OUTPUT_FILES/'//trim(outputname),status='unknown', &
- form='unformatted',action='write')
- if(iframe == 1) open(unit=12,file='OUTPUT_FILES/bin_movie.xy', &
- status='unknown',form='unformatted',action='write')
+ ! create file name and open file
+ if(OUTPUT_BINARY) then
+ if(USE_COMPONENT == 1) then
+ write(outputname,"('/bin_movie_',i6.6,'.d')") it
+ else if(USE_COMPONENT == 2) then
+ write(outputname,"('/bin_movie_',i6.6,'.N')") it
+ else if(USE_COMPONENT == 3) then
+ write(outputname,"('/bin_movie_',i6.6,'.E')") it
+ endif
+ open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown', &
+ form='unformatted',action='write')
+ if(iframe == 1) open(unit=12,file=trim(OUTPUT_FILES)//'/bin_movie.xy', &
+ status='unknown',form='unformatted',action='write')
+ else
+ if(USE_COMPONENT == 1) then
+ write(outputname,"('/ascii_movie_',i6.6,'.d')") it
+ else if(USE_COMPONENT == 2) then
+ write(outputname,"('/ascii_movie_',i6.6,'.N')") it
+ else if(USE_COMPONENT == 3) then
+ write(outputname,"('/ascii_movie_',i6.6,'.E')") it
+ endif
+ open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown', &
+ action='write')
+ if(iframe == 1) open(unit=12,file=trim(OUTPUT_FILES)//'/ascii_movie.xy', &
+ status='unknown',action='write')
+ endif
+ ! clear number of elements kept
+ ispec = 0
+
+ ! read points for all the slices
+ print *,'Writing output',outputname
+ do iproc = 0,NPROCTOT-1
+
+ ! reset point number
+ ipoin = 0
+
+ do ispecloc = 1,NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ ispec = ispec + 1
+ if(MOVIE_COARSE) then
+ ielm = ispec - 1
else
- if(USE_COMPONENT == 1) then
- write(outputname,"('ascii_movie_',i6.6,'.d')") it
- else if(USE_COMPONENT == 2) then
- write(outputname,"('ascii_movie_',i6.6,'.N')") it
- else if(USE_COMPONENT == 3) then
- write(outputname,"('ascii_movie_',i6.6,'.E')") it
- endif
- open(unit=11,file='OUTPUT_FILES/'//trim(outputname),status='unknown', &
- action='write')
- if(iframe == 1) open(unit=12,file='OUTPUT_FILES/ascii_movie.xy', &
- status='unknown',action='write')
+ ielm = (NGLLX-1)*(NGLLY-1)*(ispec-1)
endif
- ! clear number of elements kept
- ispec = 0
- ! read points for all the slices
- print *,'Writing output',outputname
- do iproc = 0,NPROCTOT-1
-
- ! reset point number
- ipoin = 0
-
- do ispecloc = 1,NEX_PER_PROC_XI*NEX_PER_PROC_ETA
- ispec = ispec + 1
+ do j = 1,NGLLY-NIT
+ do i = 1,NGLLX-NIT
if(MOVIE_COARSE) then
- ielm = ispec - 1
+ ieoff = ielm + 1
else
- ielm = (NGLLX-1)*(NGLLY-1)*(ispec-1)
+ ieoff = (ielm+(i-1)+(j-1)*(NGLLX-1))+1
endif
- do j = 1,NGLLY-NIT
- do i = 1,NGLLX-NIT
- if(MOVIE_COARSE) then
- ieoff = ielm + 1
- else
- ieoff = (ielm+(i-1)+(j-1)*(NGLLX-1))+1
- endif
+ ! point position
+ if(iframe == 1) then
+ ! gets cartesian coordinates
+ xcoord = sngl(xp(ieoff))
+ ycoord = sngl(yp(ieoff))
+ zcoord = sngl(zp(ieoff))
- ! point position
- if(iframe == 1) then
- ! gets cartesian coordinates
- xcoord = sngl(xp(ieoff))
- ycoord = sngl(yp(ieoff))
- zcoord = sngl(zp(ieoff))
+ ! location latitude/longitude (with geocentric colatitude theta )
+ call xyz_2_rthetaphi(xcoord,ycoord,zcoord,rval,thetaval,phival)
- ! location latitude/longitude (with geocentric colatitude theta )
- call xyz_2_rthetaphi(xcoord,ycoord,zcoord,rval,thetaval,phival)
+ ! converts the geocentric colatitude to a geographic colatitude
+ if(.not. ASSUME_PERFECT_SPHERE) then
+ thetaval = PI_OVER_TWO - &
+ datan(1.006760466d0*dcos(dble(thetaval))/dmax1(TINYVAL,dble(sin(thetaval))))
+ endif
- ! converts the geocentric colatitude to a geographic colatitude
- if(.not. ASSUME_PERFECT_SPHERE) then
- thetaval = PI_OVER_TWO - &
- datan(1.006760466d0*dcos(dble(thetaval))/dmax1(TINYVAL,dble(sin(thetaval))))
- endif
+ ! gets geographic latitude and longitude in degrees
+ lat = sngl(90.d0 - thetaval*RADIANS_TO_DEGREES)
+ long = sngl(phival*RADIANS_TO_DEGREES)
+ if(long > 180.0) long = long-360.0
+ endif
- ! gets geographic latitude and longitude in degrees
- lat = sngl(90.d0 - thetaval*RADIANS_TO_DEGREES)
- long = sngl(phival*RADIANS_TO_DEGREES)
- if(long > 180.0) long = long-360.0
- endif
+ ! displacement
+ disp = sngl(field_display(ieoff))
- ! displacement
- disp = sngl(field_display(ieoff))
+ ! writes displacement and latitude/longitude to corresponding files
+ if(OUTPUT_BINARY) then
+ write(11) disp
+ if(iframe == 1) write(12) long,lat
+ else
+ write(11,*) disp
+ if(iframe == 1) write(12,*) long,lat
+ endif
- ! writes displacement and latitude/longitude to corresponding files
- if(OUTPUT_BINARY) then
- write(11) disp
- if(iframe == 1) write(12) long,lat
- else
- write(11,*) disp
- if(iframe == 1) write(12,*) long,lat
- endif
+ enddo !i
+ enddo !j
+ enddo !ispecloc
+ enddo !iproc
+ close(11)
+ if(iframe == 1) close(12)
- enddo !i
- enddo !j
- enddo !ispecloc
- enddo !iproc
- close(11)
- if(iframe == 1) close(12)
-
-! end of loop and test on all the time steps for all the movie images
- endif
+ ! end of loop and test on all the time steps for all the movie images
enddo
print *,'done creating movie'
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_stability.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_stability.f90 2013-08-31 00:33:04 UTC (rev 22750)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/check_stability.f90 2013-09-02 10:52:04 UTC (rev 22751)
@@ -299,7 +299,7 @@
weekday_name(day_of_week_remote),month_name(mon_remote),day_remote,year_remote,hr_remote,minutes_remote
endif
- if (it < 100) then
+ if (it_run < 100) then
write(IMAIN,*) '************************************************************'
write(IMAIN,*) '**** BEWARE: the above time estimates are not reliable'
write(IMAIN,*) '**** because fewer than 100 iterations have been performed'
@@ -403,9 +403,18 @@
data month_name /'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'/
data weekday_name /'Sun', 'Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat'/
+ ! information about the current run only
+ it_run = it - it_begin + 1
+ nstep_run = it_end - it_begin + 1
+
! write time stamp file to give information about progression of simulation
- write(outputname,"('/timestamp',i6.6)") it
+ if(SIMULATION_TYPE == 1) then
+ write(outputname,"('/timestamp_forward',i6.6)") it
+ else
+ write(outputname,"('/timestamp_backward',i6.6)") it
+ endif
+ ! file output
open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='unknown',action='write')
write(IOUT,*) 'Time step # ',it
@@ -428,9 +437,6 @@
if( NUMBER_OF_RUNS > 1 .and. NUMBER_OF_THIS_RUN < NUMBER_OF_RUNS ) then
! this is in the case of restart files, when a given run consists of several partial runs
- ! information about the current run only
- it_run = it - it_begin + 1
- nstep_run = it_end - it_begin + 1
write(IOUT,*) 'Time steps done for this run = ',it_run,' out of ',nstep_run
write(IOUT,*) 'Time steps done in total = ',it,' out of ',NSTEP
write(IOUT,*) 'Time steps remaining for this run = ',it_end - it
@@ -451,7 +457,7 @@
write(IOUT,*) 'We have done ',sngl(100.d0*dble(it)/dble(NSTEP)),'% of that'
write(IOUT,*)
- if(it < NSTEP) then
+ if(it < it_end) then
write(IOUT,"(' The run will finish approximately on (in local time): ',a3,' ',a3,' ',i2.2,', ',i4.4,' ',i2.2,':',i2.2)") &
weekday_name(day_of_week),month_name(mon),day,year,hr,minutes
@@ -472,7 +478,7 @@
day_remote,year_remote,hr_remote,minutes_remote
endif
- if(it < 100) then
+ if(it_run < 100) then
write(IOUT,*)
write(IOUT,*) '************************************************************'
write(IOUT,*) '**** BEWARE: the above time estimates are not reliable'
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic_calling_routine.F90 2013-08-31 00:33:04 UTC (rev 22750)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_acoustic_calling_routine.F90 2013-09-02 10:52:04 UTC (rev 22751)
@@ -45,9 +45,6 @@
integer :: iphase
logical :: phase_is_inner
- ! checks
- if( SIMULATION_TYPE == 3 ) return
-
! compute internal forces in the fluid region
if(CUSTOM_REAL == SIZE_REAL) then
time = sngl((dble(it-1)*DT-t0)*scale_t_inv)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-08-31 00:33:04 UTC (rev 22750)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-09-02 10:52:04 UTC (rev 22751)
@@ -48,9 +48,6 @@
integer :: iphase
logical :: phase_is_inner
- ! checks
- if( SIMULATION_TYPE == 3 ) return
-
! ****************************************************
! big loop over all spectral elements in the solid
! ****************************************************
@@ -152,7 +149,9 @@
! add the sources
! add adjoint sources
- if (SIMULATION_TYPE == 2 ) then
+ ! note: this must remain here even when SIMULATION_TYPE == 3 because it applies to array
+ ! accel_crust_mantle rather than b_accel_crust_mantle
+ if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3 ) then
if( nadj_rec_local > 0 ) call compute_add_sources_adjoint()
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90 2013-08-31 00:33:04 UTC (rev 22750)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/finalize_simulation.f90 2013-09-02 10:52:04 UTC (rev 22751)
@@ -232,13 +232,9 @@
endif
! movies
- if(MOVIE_SURFACE .or. NOISE_TOMOGRAPHY /= 0 ) then
- deallocate(store_val_x,store_val_y,store_val_z, &
- store_val_ux,store_val_uy,store_val_uz)
- if (MOVIE_SURFACE) then
- deallocate(store_val_x_all,store_val_y_all,store_val_z_all, &
- store_val_ux_all,store_val_uy_all,store_val_uz_all)
- endif
+ if (MOVIE_SURFACE) then
+ deallocate(store_val_ux,store_val_uy,store_val_uz)
+ deallocate(store_val_ux_all,store_val_uy_all,store_val_uz_all)
endif
if(MOVIE_VOLUME) then
deallocate(nu_3dmovie)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90 2013-08-31 00:33:04 UTC (rev 22750)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90 2013-09-02 10:52:04 UTC (rev 22751)
@@ -574,13 +574,11 @@
call gather_all_i(ispec_selected_source_subset,NSOURCES_SUBSET_current_size, &
ispec_selected_source_all,NSOURCES_SUBSET_current_size,NPROCTOT_VAL)
- ! daniel debug
- !print*,'rank',myrank,'ispec:',ispec_selected_source_subset(:),'all:',ispec_selected_source_all(:,:)
-
! checks that the gather operation went well
if(myrank == 0) then
if(minval(ispec_selected_source_all(:,:)) <= 0) then
- print*,'error ispec all:',ispec_selected_source_all(:,:)
+ print*,'error ispec all:',NPROCTOT_VAL,NSOURCES_SUBSET_current_size
+ print*,ispec_selected_source_all(:,:)
call exit_MPI(myrank,'gather operation failed for source')
endif
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2013-08-31 00:33:04 UTC (rev 22750)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2013-09-02 10:52:04 UTC (rev 22751)
@@ -49,7 +49,7 @@
call prepare_timerun_convert_coord()
! allocate files to save movies
- ! for noise tomography, store_val_x/y/z/ux/uy/uz needed for 'surface movie'
+ ! for noise tomography, number of movie points (nmovie_points) needed for 'surface movie'
if( MOVIE_SURFACE .or. NOISE_TOMOGRAPHY /= 0 ) then
call prepare_timerun_movie_surface()
endif
@@ -496,39 +496,59 @@
! local parameters
integer :: ier
+ ! user output
if(myrank == 0 ) then
write(IMAIN,*) "preparing movie surface."
write(IMAIN,*)
call flush_IMAIN()
endif
-
- if(MOVIE_COARSE .and. NOISE_TOMOGRAPHY ==0) then ! only output corners !for noise tomography, must NOT be coarse
- nmovie_points = 2 * 2 * NSPEC2D_TOP(IREGION_CRUST_MANTLE)
- if(NGLLX /= NGLLY) &
+ ! only output corners
+ ! note: for noise tomography, must NOT be coarse (have to be saved on all gll points)
+ if( MOVIE_COARSE .and. NOISE_TOMOGRAPHY == 0 ) then
+ ! checks setup
+ if(NGLLX /= NGLLY) &
call exit_MPI(myrank,'MOVIE_COARSE together with MOVIE_SURFACE requires NGLLX=NGLLY')
- NIT = NGLLX - 1
+ ! number of points
+ nmovie_points = 2 * 2 * NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+ NIT = NGLLX - 1
else
- nmovie_points = NGLLX * NGLLY * NSPEC2D_TOP(IREGION_CRUST_MANTLE)
- NIT = 1
+ ! number of points
+ nmovie_points = NGLLX * NGLLY * NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+ NIT = 1
endif
- allocate(store_val_x(nmovie_points), &
- store_val_y(nmovie_points), &
- store_val_z(nmovie_points), &
- store_val_ux(nmovie_points), &
- store_val_uy(nmovie_points), &
- store_val_uz(nmovie_points),stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating movie surface arrays')
- if (MOVIE_SURFACE) then ! those arrays are not neccessary for noise tomography, so only allocate them in MOVIE_SURFACE case
- allocate(store_val_x_all(nmovie_points,0:NPROCTOT_VAL-1), &
- store_val_y_all(nmovie_points,0:NPROCTOT_VAL-1), &
- store_val_z_all(nmovie_points,0:NPROCTOT_VAL-1), &
- store_val_ux_all(nmovie_points,0:NPROCTOT_VAL-1), &
- store_val_uy_all(nmovie_points,0:NPROCTOT_VAL-1), &
- store_val_uz_all(nmovie_points,0:NPROCTOT_VAL-1),stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating movie surface all arrays')
+ ! checks exact number of points nmovie_points
+ call movie_surface_count_points()
+
+ ! those arrays are not neccessary for noise tomography, so only allocate them in MOVIE_SURFACE case
+ if( MOVIE_SURFACE ) then
+ ! writes out movie point locations to file
+ call write_movie_surface_mesh()
+
+ ! allocates movie surface arrays for wavefield values
+ allocate(store_val_ux(nmovie_points), &
+ store_val_uy(nmovie_points), &
+ store_val_uz(nmovie_points),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating movie surface arrays')
+
+ ! allocates arrays for gathering wavefield values
+ if( myrank == 0 ) then
+ ! only master needs full arrays
+ allocate(store_val_ux_all(nmovie_points,0:NPROCTOT_VAL-1), &
+ store_val_uy_all(nmovie_points,0:NPROCTOT_VAL-1), &
+ store_val_uz_all(nmovie_points,0:NPROCTOT_VAL-1),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating movie surface all arrays')
+ else
+ ! slave processes only need dummy arrays
+ allocate(store_val_ux_all(1,1), &
+ store_val_uy_all(1,1), &
+ store_val_uz_all(1,1),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating movie surface all arrays')
+ endif
endif
+
+ ! user output
if(myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) 'Movie surface:'
@@ -1522,11 +1542,11 @@
endif
allocate(noise_sourcearray(NDIM,NGLLX,NGLLY,NGLLZ,NSTEP), &
- normal_x_noise(nmovie_points), &
- normal_y_noise(nmovie_points), &
- normal_z_noise(nmovie_points), &
- mask_noise(nmovie_points), &
- noise_surface_movie(NDIM,NGLLX,NGLLY,NSPEC_TOP),stat=ier)
+ normal_x_noise(nmovie_points), &
+ normal_y_noise(nmovie_points), &
+ normal_z_noise(nmovie_points), &
+ mask_noise(nmovie_points), &
+ noise_surface_movie(NDIM,NGLLX,NGLLY,NSPEC_TOP),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating noise arrays')
noise_sourcearray(:,:,:,:,:) = 0._CUSTOM_REAL
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90 2013-08-31 00:33:04 UTC (rev 22750)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90 2013-09-02 10:52:04 UTC (rev 22751)
@@ -713,14 +713,9 @@
! to save movie frames
integer :: nmovie_points,NIT
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
- store_val_x,store_val_y,store_val_z, &
- store_val_ux,store_val_uy,store_val_uz
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_ux,store_val_uy,store_val_uz
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: store_val_ux_all,store_val_uy_all,store_val_uz_all
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
- store_val_x_all,store_val_y_all,store_val_z_all, &
- store_val_ux_all,store_val_uy_all,store_val_uz_all
-
! to save movie volume
double precision :: scalingval
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_output.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_output.f90 2013-08-31 00:33:04 UTC (rev 22750)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_output.f90 2013-09-02 10:52:04 UTC (rev 22751)
@@ -52,8 +52,11 @@
! gets resulting array values onto CPU
if( GPU_MODE ) then
! transfers whole fields
- call transfer_displ_cm_from_device(NDIM*NGLOB_CRUST_MANTLE,displ_crust_mantle,Mesh_pointer)
- call transfer_veloc_cm_from_device(NDIM*NGLOB_CRUST_MANTLE,veloc_crust_mantle,Mesh_pointer)
+ if(MOVIE_VOLUME_TYPE == 5) then
+ call transfer_displ_cm_from_device(NDIM*NGLOB_CRUST_MANTLE,displ_crust_mantle,Mesh_pointer)
+ else
+ call transfer_veloc_cm_from_device(NDIM*NGLOB_CRUST_MANTLE,veloc_crust_mantle,Mesh_pointer)
+ endif
endif
! save velocity here to avoid static offset on displacement for movies
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_surface.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_surface.f90 2013-08-31 00:33:04 UTC (rev 22750)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_surface.f90 2013-09-02 10:52:04 UTC (rev 22751)
@@ -25,6 +25,135 @@
!
!=====================================================================
+
+ subroutine movie_surface_count_points()
+
+ use specfem_par
+ use specfem_par_crustmantle,only: NSPEC_TOP
+ use specfem_par_movie,only: NIT,nmovie_points
+
+ implicit none
+
+ ! local parameters
+ integer :: ipoin,ispec2D,i,j,k,npoin
+
+ ! gets number of points on surface mesh
+ ipoin = 0
+ do ispec2D = 1, NSPEC_TOP ! NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+ k = NGLLZ
+ ! loop on all the points inside the element
+ do j = 1,NGLLY,NIT
+ do i = 1,NGLLX,NIT
+ ipoin = ipoin + 1
+ enddo
+ enddo
+ enddo
+ npoin = ipoin
+
+ ! checks
+ if( npoin /= nmovie_points ) then
+ print*,'error: movie points collected ',npoin,'not equal to calculated :',nmovie_points
+ call exit_mpi(myrank,'error confusing number of movie points')
+ endif
+
+ end subroutine movie_surface_count_points
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine write_movie_surface_mesh()
+
+! writes out movie point locations to file
+
+ use specfem_par
+ use specfem_par_crustmantle
+ use specfem_par_movie
+
+ implicit none
+
+ ! local parameters
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_x,store_val_y,store_val_z
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: store_val_x_all,store_val_y_all,store_val_z_all
+ integer :: ipoin,ispec2D,ispec,i,j,k,ier,iglob,npoin
+ character(len=150) :: outputname
+
+ ! allocates movie surface arrays
+ allocate(store_val_x(nmovie_points), &
+ store_val_y(nmovie_points), &
+ store_val_z(nmovie_points),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating movie surface location arrays')
+
+ ! allocates arrays for gathering movie point locations
+ if( myrank == 0 ) then
+ ! only master needs full arrays
+ allocate(store_val_x_all(nmovie_points,0:NPROCTOT_VAL-1), &
+ store_val_y_all(nmovie_points,0:NPROCTOT_VAL-1), &
+ store_val_z_all(nmovie_points,0:NPROCTOT_VAL-1),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating movie surface all arrays')
+ else
+ ! slave processes only need dummy arrays
+ allocate(store_val_x_all(1,1), &
+ store_val_y_all(1,1), &
+ store_val_z_all(1,1),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating movie surface all arrays')
+ endif
+
+ ! gets coordinates of surface mesh
+ ipoin = 0
+ do ispec2D = 1, NSPEC_TOP ! NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+ ispec = ibelm_top_crust_mantle(ispec2D)
+ ! in case of global, NCHUNKS_VAL == 6 simulations, be aware that for
+ ! the cubed sphere, the mapping changes for different chunks,
+ ! i.e. e.g. x(1,1) and x(5,5) flip left and right sides of the elements in geographical coordinates.
+ ! for future consideration, like in create_movie_GMT_global.f90 ...
+ k = NGLLZ
+ ! loop on all the points inside the element
+ do j = 1,NGLLY,NIT
+ do i = 1,NGLLX,NIT
+ ipoin = ipoin + 1
+ ! stores values
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+ store_val_x(ipoin) = xstore_crust_mantle(iglob) ! <- radius r (normalized)
+ store_val_y(ipoin) = ystore_crust_mantle(iglob) ! <- colatitude theta (in radian)
+ store_val_z(ipoin) = zstore_crust_mantle(iglob) ! <- longitude phi (in radian)
+ enddo
+ enddo
+ enddo
+ npoin = ipoin
+ if( npoin /= nmovie_points ) call exit_mpi(myrank,'error number of movie points not equal to nmovie_points')
+
+ ! gather info on master proc
+ call gather_all_cr(store_val_x,nmovie_points,store_val_x_all,nmovie_points,NPROCTOT_VAL)
+ call gather_all_cr(store_val_y,nmovie_points,store_val_y_all,nmovie_points,NPROCTOT_VAL)
+ call gather_all_cr(store_val_z,nmovie_points,store_val_z_all,nmovie_points,NPROCTOT_VAL)
+
+ ! save movie data locations to disk in home directory
+ if(myrank == 0) then
+
+ ! outputs movie point locations to moviedata_xyz.bin file
+ outputname = "/moviedata_xyz.bin"
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//trim(outputname), &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening moviedata_xyz.bin file')
+
+ ! point coordinates
+ ! (given as r theta phi for geocentric coordinate system)
+ write(IOUT) store_val_x_all
+ write(IOUT) store_val_y_all
+ write(IOUT) store_val_z_all
+ close(IOUT)
+ endif
+
+ deallocate(store_val_x_all,store_val_y_all,store_val_z_all)
+ deallocate(store_val_x,store_val_y,store_val_z)
+
+ end subroutine write_movie_surface_mesh
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine write_movie_surface()
use specfem_par
@@ -39,7 +168,7 @@
! by default: save velocity here to avoid static offset on displacement for movies
- ! get coordinates of surface mesh and surface displacement
+ ! gets coordinates of surface mesh and surface displacement
ipoin = 0
do ispec2D = 1, NSPEC_TOP ! NSPEC2D_TOP(IREGION_CRUST_MANTLE)
ispec = ibelm_top_crust_mantle(ispec2D)
@@ -55,9 +184,8 @@
do i = 1,NGLLX,NIT
ipoin = ipoin + 1
iglob = ibool_crust_mantle(i,j,k,ispec)
- store_val_x(ipoin) = xstore_crust_mantle(iglob) ! <- radius r (normalized)
- store_val_y(ipoin) = ystore_crust_mantle(iglob) ! <- colatitude theta (in radian)
- store_val_z(ipoin) = zstore_crust_mantle(iglob) ! <- longitude phi (in radian)
+
+ ! wavefield values
if(MOVIE_VOLUME_TYPE == 5) then
! stores displacement
store_val_ux(ipoin) = displ_crust_mantle(1,iglob)*scale_displ
@@ -72,17 +200,13 @@
enddo
enddo
-
enddo
! gather info on master proc
- ispec = nmovie_points
- call gather_all_cr(store_val_x,ispec,store_val_x_all,ispec,NPROCTOT_VAL)
- call gather_all_cr(store_val_y,ispec,store_val_y_all,ispec,NPROCTOT_VAL)
- call gather_all_cr(store_val_z,ispec,store_val_z_all,ispec,NPROCTOT_VAL)
- call gather_all_cr(store_val_ux,ispec,store_val_ux_all,ispec,NPROCTOT_VAL)
- call gather_all_cr(store_val_uy,ispec,store_val_uy_all,ispec,NPROCTOT_VAL)
- call gather_all_cr(store_val_uz,ispec,store_val_uz_all,ispec,NPROCTOT_VAL)
+ ! wavefield
+ call gather_all_cr(store_val_ux,nmovie_points,store_val_ux_all,nmovie_points,NPROCTOT_VAL)
+ call gather_all_cr(store_val_uy,nmovie_points,store_val_uy_all,nmovie_points,NPROCTOT_VAL)
+ call gather_all_cr(store_val_uz,nmovie_points,store_val_uz_all,nmovie_points,NPROCTOT_VAL)
! save movie data to disk in home directory
if(myrank == 0) then
@@ -91,12 +215,13 @@
status='unknown',form='unformatted',action='write',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening moviedata file')
- write(IOUT) store_val_x_all
- write(IOUT) store_val_y_all
- write(IOUT) store_val_z_all
+ ! -> find movie point locations stored in file moviedata_xyz.bin
+
+ ! wavefield values
write(IOUT) store_val_ux_all
write(IOUT) store_val_uy_all
write(IOUT) store_val_uz_all
+
close(IOUT)
endif
More information about the CIG-COMMITS
mailing list