[cig-commits] r17262 - in seismo/3D/ADJOINT_TOMO/iterate_adj/pangu: model_vp_vs model_vp_vs/src smooth smooth/src sum_kernel sum_kernel/src
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Tue Oct 12 15:46:52 PDT 2010
Author: danielpeter
Date: 2010-10-12 15:46:52 -0700 (Tue, 12 Oct 2010)
New Revision: 17262
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/go_globe_pbs.bash
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/readme_globe
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/run_globe.bash
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/Makefile_globe
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_iso.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso_iso.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/go_globe_pbs.bash
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/readme_globe
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/run_globe.bash
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/Makefile_globe
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/smooth_sem_globe.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/go_globe_pbs.bash
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/readme_globe
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/run_globe.bash
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/Makefile_globe
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/sum_kernels_globe.f90
Log:
adds routines for global inversions to iterate_adj package for sum_kernel,smooth and model_vp_vs; see readme_globe files in corresponding directories for more infos
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/go_globe_pbs.bash
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/go_globe_pbs.bash (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/go_globe_pbs.bash 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,38 @@
+#!/bin/bash
+#PBS -S /bin/bash
+
+## job name and output file
+#PBS -N model_update
+#PBS -j oe
+#PBS -o OUTPUT_FILES/job.o
+
+###########################################################
+# USER PARAMETERS
+
+## 144 CPUs ( 18*8 ), walltime 5 hour
+#PBS -l nodes=18:ppn=8,walltime=5:00:00
+
+numnodes=144
+
+# model update percentage 3%
+percentage=0.03
+
+
+###########################################################
+
+# (takes about 10 min...)
+
+echo `date`
+
+cd $PBS_O_WORKDIR
+
+# obtain lsf job information
+cat $PBS_NODEFILE > OUTPUT_FILES/compute_nodes
+echo "$PBS_JOBID" > OUTPUT_FILES/jobid
+
+# steepest descent step
+mpiexec -np $numnodes $PWD/src/add_model_globe_iso $percentage
+
+
+echo "done successfully"
+echo `date`
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/go_globe_pbs.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/readme_globe
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/readme_globe (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/readme_globe 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,75 @@
+----------------------------
+readme
+----------------------------
+
+add_model_globe_iso:
+
+ takes ISOTROPIC model and ISOTROPIC kernels and
+ makes a steepest descent model update
+
+ the new models vp & vs are then computed from the old ones (in INPU_MODEL/), assuming
+ that the gradients/kernels relate relative model perturbations dln(m/m0) to traveltime/multitaper measurements.
+
+ setup:
+
+ input:
+ - step_fac : step length to update the models, f.e. 0.03 for plusminus 3%
+
+ - INPUT_MODEL/ contains:
+ proc000***_reg1_vs.bin & proc000***_reg1_vp.bin etc.
+ - INPUT_GRADIENT/ contains:
+ proc000***_reg1_alpha_kernel_smooth.bin & proc000***_reg1_beta_kernel_smooth.bin etc.
+ - topo/ contains:
+ proc000***_reg1_solver_data_1.bin
+
+ new models are stored in
+ - OUTPUT_MODEL/ as proc000***_reg1_vp_new.bin and proc000***_reg1_vs_new.bin
+
+
+ usage: > ./add_model_globe_iso 0.3
+
+
+
+add_model_globe_tiso_iso:
+
+ takes TRANSVERSE ISOTROPIC model and ISOTROPIC kernels and
+ makes a steepest descent model update
+
+add_model_globe_tiso:
+
+ takes TRANSVERSE ISOTROPIC model and TRANSVERSE ISOTROPIC kernels and
+ makes a steepest descent model update
+
+
+
+steps:
+------
+
+0. copy file SPECFEM3D_GLOBE/constants.h
+ over to: src/constants_globe.h
+
+ copy file SPECFEM3D_GLOBE/precision.h
+ over to: src/precision_globe.h
+
+ copy file SPECFEM3D_GLOBE/OUTPUT_FILES/values_from_mesher.h
+ over to: src/values_from_mesher_globe.h
+ ( when generated for a kernel run ...)
+
+1. compile add_model_globe_iso:
+ > cd src/
+ > make -f Makefile_globe
+
+2. setup input model links:
+ > ln -s /my_input_model/proc***.bin INPUT_MODEL/
+
+3. setup gradient directory:
+ > ln -s /my_summed_smoothed_event_kernel_directory INPUT_GRADIENT
+
+4. run job:
+
+ modify step_length in go_globe_pbs.bash,
+ then
+ > run_globe.bash
+
+
+
\ No newline at end of file
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/run_globe.bash
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/run_globe.bash (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/run_globe.bash 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,16 @@
+#!/bin/bash
+
+# cleans outputs
+mkdir -p OUTPUT_FILES
+rm -f OUTPUT_FILES/*
+rm -f OUTPUT_MODEL/*.*
+
+# compiles executable
+cd src/
+make -f Makefile_globe
+cd ..
+
+# submits job
+date
+qsub go_globe_pbs.bash
+
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/run_globe.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/Makefile_globe
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/Makefile_globe (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/Makefile_globe 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,22 @@
+FFLAGS=-O3 -assume byterecl
+
+all: add_model_globe_iso \
+ add_model_globe_tiso \
+ add_model_globe_tiso_iso
+
+add_model_globe_iso: add_model_globe_iso.f90
+ mpif90 -o add_model_globe_iso $(FFLAGS) add_model_globe_iso.f90 gll_library.f90 exit_mpi.f90
+
+add_model_globe_tiso: add_model_globe_tiso.f90
+ mpif90 -o add_model_globe_tiso $(FFLAGS) add_model_globe_tiso.f90 gll_library.f90 exit_mpi.f90
+
+add_model_globe_tiso_iso: add_model_globe_tiso_iso.f90
+ mpif90 -o add_model_globe_tiso_iso $(FFLAGS) add_model_globe_tiso_iso.f90 gll_library.f90 exit_mpi.f90
+
+
+
+clean:
+ rm -f add_model_globe_iso \
+ add_model_globe_tiso \
+ add_model_globe_tiso_iso \
+ *.o *.mod
\ No newline at end of file
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_iso.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_iso.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_iso.f90 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,661 @@
+! add_model_globe_iso
+!
+! this program can be used to update ISOTROPIC model files with
+! (smoothed & summed) event kernels.
+! the kernels are given for isotropic parameters (alpha,beta,rho) or ( bulk_c,beta,rho).
+!
+! the algorithm uses a steepest descent method with a step length
+! determined by the given maximum update percentage.
+!
+! input:
+! - step_fac : step length to update the models, f.e. 0.03 for plusminus 3%
+!
+! setup:
+!
+!- INPUT_MODEL/ contains:
+! proc000***_reg1_vs.bin &
+! proc000***_reg1_vp.bin &
+! proc000***_reg1_rho.bin
+!
+!- INPUT_GRADIENT/ contains:
+! proc000***_reg1_bulk_c_kernel_smooth.bin &
+! proc000***_reg1_bulk_beta_kernel_smooth.bin &
+! proc000***_reg1_rho_kernel_smooth.bin
+! or
+! proc000***_reg1_alpha_kernel_smooth.bin &
+! proc000***_reg1_beta_kernel_smooth.bin &
+! proc000***_reg1_rho_kernel_smooth.bin
+!
+!- topo/ contains:
+! proc000***_reg1_solver_data_1.bin
+!
+! new models are stored in
+!- OUTPUT_MODEL/ as
+! proc000***_reg1_vp_new.bin and
+! proc000***_reg1_vs_new.bin and
+! proc000***_reg1_rho_new.bin and
+!
+! USAGE: e.g. ./add_model_globe_iso 0.3
+
+program add_model
+
+ implicit none
+
+ include 'mpif.h'
+ include 'constants_globe.h'
+ include 'precision_globe.h'
+ include 'values_from_mesher_globe.h'
+
+ ! ======================================================
+
+ ! USER PARAMETERS
+
+ ! by default, this algorithm uses (bulk,bulk_beta,rho) kernels to update vp,vs,rho
+ ! if you prefer using (alpha,beta,rho) kernels, set this flag to true
+ logical, parameter :: USE_ALPHA_BETA_RHO = .false.
+
+ ! ignore rho kernel, but use density perturbations as a scaling of Vs perturbations
+ logical, parameter :: USE_RHO_SCALING = .false.
+
+ ! in case of rho scaling, specifies density scaling factor with shear perturbations
+ ! see e.g. Montagner & Anderson (1989), Panning & Romanowicz (2006)
+ real(kind=CUSTOM_REAL),parameter :: RHO_SCALING = 0.33_CUSTOM_REAL
+
+ ! ======================================================
+
+ integer, parameter :: NSPEC = NSPEC_CRUST_MANTLE
+ integer, parameter :: NGLOB = NGLOB_CRUST_MANTLE
+
+ character(len=150) :: sline, m_file, fname
+
+ ! model values, model updates and kernels
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ model_vp, model_vs, model_rho, &
+ model_dA, model_dB, model_dR, &
+ kernel_a, kernel_b, kernel_rho, total_model
+
+ ! steepest descent lengths
+ real(kind=CUSTOM_REAL) :: step_fac,step_length
+
+ integer :: nfile, myrank, sizeprocs, ier
+ integer :: i, j, k,ispec, iglob, ishell, n, it, j1, ib, npts_sem, ios
+
+ ! statistics
+ real(kind=CUSTOM_REAL) :: min_vp,min_vs,max_vp,max_vs, &
+ min_rho,max_rho,max,minmax(4),vs_sum,vp_sum,rho_sum
+ character(len=150) :: s_step_fac
+
+ ! volume
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: x, y, z
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: volume_spec
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: volume_glob,volume_glob_sum
+
+ ! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: xigll, wxgll
+ double precision, dimension(NGLLY) :: yigll, wygll
+ double precision, dimension(NGLLZ) :: zigll, wzgll
+
+ ! array with all the weights in the cube
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+
+ ! jacobian
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: jacobian
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl, &
+ jacobianl,volumel
+
+ ! calculates volume associated with (global) point
+ logical,parameter:: VOLUME_GLOBALPOINT = .true.
+ real(kind=CUSTOM_REAL) :: kernel_integral_alpha,kernel_integral_beta,kernel_integral_rho
+ real(kind=CUSTOM_REAL) :: integral_alpha_sum,integral_beta_sum,integral_rho_sum
+
+ ! ============ program starts here =====================
+
+ ! initialize the MPI communicator and start the NPROCTOT MPI processes
+ call MPI_INIT(ier)
+ call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
+ call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
+
+ if( sizeprocs /= nchunks_val*nproc_xi_val*nproc_eta_val ) then
+ print*,'sizeprocs:',sizeprocs,nchunks_val,nproc_xi_val,nproc_eta_val
+ call exit_mpi('error number sizeprocs')
+ endif
+
+ model_vp = 0.0_CUSTOM_REAL
+ model_vs = 0.0_CUSTOM_REAL
+ model_rho = 0.0_CUSTOM_REAL
+ model_dA = 0.0_CUSTOM_REAL
+ model_dB = 0.0_CUSTOM_REAL
+ model_dR = 0.0_CUSTOM_REAL
+ kernel_a = 0.0_CUSTOM_REAL
+ kernel_b = 0.0_CUSTOM_REAL
+ kernel_rho = 0.0_CUSTOM_REAL
+
+ !--------------------------------------------------------
+
+ ! subjective step length to multiply to the gradient
+ ! e.g. step_fac = 0.03
+
+ call getarg(1,s_step_fac)
+
+ if (trim(s_step_fac) == '') then
+ call exit_MPI(myrank,'Usage: add_model_globe_iso step_factor')
+ endif
+
+ ! read in parameter information
+ read(s_step_fac,*) step_fac
+ if( abs(step_fac) < 1.e-10) then
+ print*,'error: step factor ',step_fac
+ call exit_MPI(myrank,'error step factor')
+ endif
+
+ if (myrank == 0) then
+ print*,'defaults'
+ print*,' NPROC_XI , NPROC_ETA: ',nproc_xi_val,nproc_eta_val
+ print*,' NCHUNKS: ',nchunks_val
+ print*
+ print*,'model update for vs & vp & rho'
+ print*,' step_fac = ',step_fac
+ print*
+ if( USE_ALPHA_BETA_RHO ) then
+ print*,'kernel parameterization: (alpha,beta,rho)'
+ else
+ print*,'kernel parameterization: (bulk,beta,rho)'
+ endif
+ print*
+ if( USE_RHO_SCALING ) then
+ print*,'scaling rho perturbations'
+ print*
+ endif
+
+ open(20,file='step_fac',status='unknown',action='write')
+ write(20,'(1e24.12)') step_fac
+ close(20)
+ endif
+
+ !-----------------------------------------------------
+ ! reads in current vp & vs & rho model
+ ! vp model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vp.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_vp(:,:,:,1:nspec)
+ close(12)
+
+ ! vs model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vs.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_vs(:,:,:,1:nspec)
+ close(12)
+
+ ! rho model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_rho.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_rho(:,:,:,1:nspec)
+ close(12)
+
+ ! statistics
+ call mpi_reduce(minval(model_vp),min_vp,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_vp),max_vp,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_vs),min_vs,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_vs),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_rho),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_rho),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'initial models:'
+ print*,' vs min/max: ',min_vs,max_vs
+ print*,' vp min/max: ',min_vp,max_vp
+ print*,' rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+ ! reads in smoothed (& summed) event kernel
+ if( USE_ALPHA_BETA_RHO ) then
+ ! reads in alpha kernel
+ fname = 'alpha_kernel_smooth'
+ else
+ ! reads in bulk_c kernel
+ fname = 'bulk_c_kernel_smooth'
+ endif
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) kernel_a(:,:,:,1:nspec)
+ close(12)
+
+ ! beta kernel
+ if( USE_ALPHA_BETA_RHO ) then
+ ! reads in beta kernel
+ fname = 'beta_kernel_smooth'
+ else
+ ! reads in bulk_beta kernel
+ fname = 'bulk_beta_kernel_smooth'
+ endif
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) kernel_b(:,:,:,1:nspec)
+ close(12)
+
+ ! rho kernel
+ if( USE_RHO_SCALING ) then
+
+ ! uses scaling relation with shear perturbations
+ kernel_rho(:,:,:,:) = RHO_SCALING * kernel_b(:,:,:,:)
+
+ else
+
+ ! uses rho kernel
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_rho_kernel_smooth.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) kernel_rho(:,:,:,1:nspec)
+ close(12)
+ endif
+
+ ! statistics
+ call mpi_reduce(minval(kernel_a),min_vp,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(kernel_a),max_vp,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(kernel_b),min_vs,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(kernel_b),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(kernel_rho),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(kernel_rho),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'initial kernels:'
+ print*,' beta min/max: ',min_vs,max_vs
+ print*,' alpha min/max: ',min_vp,max_vp
+ print*,' rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+
+ !--------------------------------------------------------
+
+ if( VOLUME_GLOBALPOINT ) then
+ ! computes volume element associated with points
+ ! GLL points
+ wgll_cube = 0.0
+ call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ wgll_cube(i,j,k) = wxgll(i)*wygll(j)*wzgll(k)
+ enddo
+ enddo
+ enddo
+
+ ! global addressing
+ write(m_file,'(a,i6.6,a)') 'topo/proc',myrank,'_reg1_solver_data_2.bin'
+ open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(11) x(1:nglob)
+ read(11) y(1:nglob)
+ read(11) z(1:nglob)
+ read(11) ibool(:,:,:,1:nspec)
+ close(11)
+
+ ! builds jacobian
+ write(m_file,'(a,i6.6,a)') 'topo/proc',myrank,'_reg1_solver_data_1.bin'
+ open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(11) xix
+ read(11) xiy
+ read(11) xiz
+ read(11) etax
+ read(11) etay
+ read(11) etaz
+ read(11) gammax
+ read(11) gammay
+ read(11) gammaz
+ close(11)
+
+ jacobian = 0.0
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ ! gets derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ ! computes the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+ jacobian(i,j,k,ispec) = jacobianl
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! volume associated with global point
+ volume_glob = 0.0
+ kernel_integral_alpha = 0._CUSTOM_REAL
+ kernel_integral_beta = 0._CUSTOM_REAL
+ kernel_integral_rho = 0._CUSTOM_REAL
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+ if( iglob == 0 ) then
+ print*,'iglob zero',i,j,k,ispec
+ print*
+ print*,'ibool:',ispec
+ print*,ibool(:,:,:,ispec)
+ print*
+ call exit_MPI('error ibool')
+ endif
+
+ ! volume associated with GLL point
+ volumel = jacobian(i,j,k,ispec)*wgll_cube(i,j,k)
+
+ ! kernel integration: for each element
+ kernel_integral_alpha = kernel_integral_alpha &
+ + volumel * kernel_a(i,j,k,ispec)
+
+ kernel_integral_beta = kernel_integral_beta &
+ + volumel * kernel_b(i,j,k,ispec)
+
+ kernel_integral_rho = kernel_integral_rho &
+ + volumel * kernel_rho(i,j,k,ispec)
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! statistics
+ ! kernel integration: for whole volume
+ call mpi_reduce(kernel_integral_alpha,integral_alpha_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(kernel_integral_beta,integral_beta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(kernel_integral_rho,integral_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'integral kernels:'
+ print*,' beta : ',integral_beta_sum
+ print*,' a : ',integral_alpha_sum
+ print*,' rho : ',integral_rho_sum
+ print*
+ endif
+
+ endif
+
+ ! calculates gradient
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ ! gradient in direction of steepest descent
+
+ ! for vp
+ model_dA(i,j,k,ispec) = - kernel_a(i,j,k,ispec)
+
+ ! for shear
+ model_dB(i,j,k,ispec) = - kernel_b(i,j,k,ispec)
+
+ ! for rho
+ model_dR(i,j,k,ispec) = - kernel_rho(i,j,k,ispec)
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! statistics
+ call mpi_reduce(minval(model_dA),min_vp,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dA),max_vp,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_dB),min_vs,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dB),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_dR),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dR),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'initial gradients:'
+ print*,' beta min/max : ',min_vs,max_vs
+ print*,' a min/max: ',min_vp,max_vp
+ print*,' rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+ ! master determines step length based on maximum gradient value (either vp or vs)
+ if( myrank == 0 ) then
+ minmax(1) = abs(min_vs)
+ minmax(2) = abs(max_vs)
+
+ minmax(3) = abs(min_vp)
+ minmax(4) = abs(max_vp)
+
+ max = maxval(minmax)
+ step_length = step_fac/max
+
+ print*,' step length : ',step_length,max
+ print*
+ endif
+ call mpi_bcast(step_length,1,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+
+
+ ! gradient length
+ max_vp = sum( model_dA * model_dA )
+ max_vs = sum( model_dB * model_dB )
+ max_rho = sum( model_dR * model_dR )
+ max_vp = sqrt(max_vp)
+ max_vs = sqrt(max_vs)
+ max_rho = sqrt(max_rho)
+
+ call mpi_reduce(max_vp,vp_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(max_vs,vs_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(max_rho,rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,' initial beta length : ',vs_sum
+ print*,' initial a length: ',vp_sum
+ print*,' initial rho length: ',rho_sum
+ print*
+ endif
+
+ ! model updates:
+ ! multiply model updates by a subjective factor that will change the step
+
+ model_dA(:,:,:,:) = step_length * model_dA(:,:,:,:)
+ model_dB(:,:,:,:) = step_length * model_dB(:,:,:,:)
+ model_dR(:,:,:,:) = step_length * model_dR(:,:,:,:)
+
+
+ ! statistics
+ call mpi_reduce(minval(model_dA),min_vp,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dA),max_vp,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_dB),min_vs,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dB),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_dR),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dR),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'scaled gradients:'
+ print*,' beta min/max : ',min_vs,max_vs
+ print*,' a min/max: ',min_vp,max_vp
+ print*,' rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+ !-----------------------------------------------------
+ ! new model:
+ ! computes new model values for alpha, beta and rho
+ ! and stores new model files
+
+ ! S wavespeed model
+ total_model = 0._CUSTOM_REAL
+ total_model = model_vs * exp( model_dB )
+
+ call mpi_reduce(maxval(total_model),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vs,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ fname = 'vs_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) total_model(:,:,:,1:nspec)
+ close(12)
+
+ ! Rho density model
+ total_model = 0._CUSTOM_REAL
+ total_model = model_rho * exp( model_dR )
+
+ call mpi_reduce(maxval(total_model),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ fname = 'rho_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) total_model(:,:,:,1:nspec)
+ close(12)
+
+ ! P wavespeed model
+ total_model = 0._CUSTOM_REAL
+ if( USE_ALPHA_BETA_RHO ) then
+ ! new vp values use alpha model update
+ total_model = model_vp * exp( model_dA )
+ else
+ ! new vp values use bulk model update:
+ ! this is based on vp_new = sqrt( bulk_new**2 + 4/3 vs_new**2 )
+ total_model = sqrt( model_vp**2 * exp(2.0*model_dA) + FOUR_THIRDS * model_vs**2 *( &
+ exp(2.0*model_dB) - exp(2.0*model_dA) ) )
+ endif
+ call mpi_reduce(maxval(total_model),max_vp,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vp,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ fname = 'vp_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) total_model(:,:,:,1:nspec)
+ close(12)
+
+
+ if( myrank == 0 ) then
+ print*,'new models:'
+ print*,' vs min/max: ',min_vs,max_vs
+ print*,' vp min/max: ',min_vp,max_vp
+ print*,' rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+
+ ! P & S & rho model update
+ ! stores the linear approximations of the model updates
+
+ ! stores relative Vp model perturbations
+ where( model_vp /= 0.0 ) total_model = ( total_model - model_vp) / model_vp
+
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvpvp.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model(:,:,:,1:nspec)
+ close(12)
+ call mpi_reduce(maxval(total_model),max_vp,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vp,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! relative S model perturbations
+ total_model = 0._CUSTOM_REAL
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ if ( model_vs(i,j,k,ispec) > 1.e-10 ) then
+ total_model(i,j,k,ispec) = ( model_vs(i,j,k,ispec)*exp(model_dB(i,j,k,ispec)) &
+ - model_vs(i,j,k,ispec) ) / model_vs(i,j,k,ispec)
+ endif
+ enddo
+ enddo
+ enddo
+ enddo
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvsvs.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model(:,:,:,1:nspec)
+ close(12)
+ call mpi_reduce(maxval(total_model),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vs,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! relative rho model perturbations
+ total_model = 0._CUSTOM_REAL
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ if ( model_rho(i,j,k,ispec) > 1.e-10 ) then
+ total_model(i,j,k,ispec) = ( model_rho(i,j,k,ispec)*exp(model_dR(i,j,k,ispec)) &
+ - model_rho(i,j,k,ispec) ) / model_rho(i,j,k,ispec)
+ endif
+ enddo
+ enddo
+ enddo
+ enddo
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_drhorho.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model(:,:,:,1:nspec)
+ close(12)
+ call mpi_reduce(maxval(total_model),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+
+ if( myrank == 0 ) then
+ print*,'relative Vs & Vp update:'
+ print*,' dvs/vs min/max: ',min_vs,max_vs
+ print*,' dvp/vp min/max: ',min_vp,max_vp
+ print*,' drho/rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+
+ !-----------------------------------------------------
+
+ ! stop all the MPI processes, and exit
+ call MPI_FINALIZE(ier)
+
+end program add_model
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso.f90 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,1052 @@
+! add_model_globe_tiso
+!
+! this program can be used to update TRANSVERSE ISOTROPIC model files
+! based on smoothed event kernels.
+! the kernels are given for tranverse isotropic parameters (bulk_c,bulk_betav,bulk_betah,eta).
+!
+! the algorithm uses a steepest descent method with a step length
+! determined by the given maximum update percentage.
+!
+! input:
+! - step_fac : step length to update the models, f.e. 0.03 for plusminus 3%
+!
+! setup:
+!
+!- INPUT_MODEL/ contains:
+! proc000***_reg1_vsv.bin &
+! proc000***_reg1_vsh.bin &
+! proc000***_reg1_vpv.bin &
+! proc000***_reg1_vph.bin &
+! proc000***_reg1_eta.bin &
+! proc000***_reg1_rho.bin
+!
+!- INPUT_GRADIENT/ contains:
+! proc000***_reg1_bulk_c_kernel_smooth.bin &
+! proc000***_reg1_bulk_betav_kernel_smooth.bin &
+! proc000***_reg1_bulk_betah_kernel_smooth.bin &
+! proc000***_reg1_eta_kernel_smooth.bin
+!
+!- topo/ contains:
+! proc000***_reg1_solver_data_1.bin
+!
+! new models are stored in
+!- OUTPUT_MODEL/ as
+! proc000***_reg1_vpv_new.bin &
+! proc000***_reg1_vph_new.bin &
+! proc000***_reg1_vsv_new.bin &
+! proc000***_reg1_vsh_new.bin &
+! proc000***_reg1_eta_new.bin &
+! proc000***_reg1_rho_new.bin
+!
+! USAGE: ./add_model_globe_tiso 0.3
+
+module model_update
+
+ include 'mpif.h'
+ include 'constants_globe.h'
+ include 'precision_globe.h'
+ include 'values_from_mesher_globe.h'
+
+ ! ======================================================
+
+ ! density scaling factor with shear perturbations
+ ! see e.g. Montagner & Anderson (1989), Panning & Romanowicz (2006)
+ real(kind=CUSTOM_REAL),parameter :: RHO_SCALING = 0.33_CUSTOM_REAL
+
+ ! constraint on eta model
+ real(kind=CUSTOM_REAL),parameter :: LIMIT_ETA_MIN = 0.5_CUSTOM_REAL
+ real(kind=CUSTOM_REAL),parameter :: LIMIT_ETA_MAX = 1.5_CUSTOM_REAL
+
+ ! ======================================================
+
+ integer, parameter :: NSPEC = NSPEC_CRUST_MANTLE
+ integer, parameter :: NGLOB = NGLOB_CRUST_MANTLE
+
+ ! transverse isotropic model files
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ model_vpv,model_vph,model_vsv,model_vsh,model_eta,model_rho
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ model_vpv_new,model_vph_new,model_vsv_new,model_vsh_new,model_eta_new,model_rho_new
+
+ ! model updates
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ model_dbulk,model_dbetah,model_dbetav,model_deta
+
+ ! kernels
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ kernel_bulk,kernel_betav,kernel_betah,kernel_eta
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: total_model
+
+ ! volume
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: x, y, z
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: volume_spec
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: volume_glob,volume_glob_sum
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
+ integer, dimension(NSPEC) :: idoubling
+
+ ! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: xigll, wxgll
+ double precision, dimension(NGLLY) :: yigll, wygll
+ double precision, dimension(NGLLZ) :: zigll, wzgll
+
+ ! array with all the weights in the cube
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+
+ ! gradient vector norm ( v^T * v )
+ real(kind=CUSTOM_REAL) :: norm_bulk,norm_betav,norm_betah,norm_eta
+ real(kind=CUSTOM_REAL) :: norm_bulk_sum,norm_betav_sum, &
+ norm_betah_sum,norm_eta_sum
+
+ ! model update length
+ real(kind=CUSTOM_REAL) :: step_fac,step_length
+
+ real(kind=CUSTOM_REAL) :: min_vpv,min_vph,min_vsv,min_vsh, &
+ max_vpv,max_vph,max_vsv,max_vsh,min_eta,max_eta,min_bulk,max_bulk, &
+ min_rho,max_rho,max,minmax(4)
+
+ real(kind=CUSTOM_REAL) :: betav1,betah1,betav0,betah0,rho1,rho0, &
+ betaiso1,betaiso0,eta1,eta0,alphav1,alphav0,alphah1,alphah0
+ real(kind=CUSTOM_REAL) :: dbetaiso,dbulk
+
+ integer :: nfile, myrank, sizeprocs, ier
+ integer :: i, j, k,ispec, iglob, ishell, n, it, j1, ib, npts_sem, ios
+ character(len=150) :: sline, m_file, fname
+
+
+end module model_update
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+program add_model
+
+ use model_update
+
+ implicit none
+
+ ! ============ program starts here =====================
+
+ ! initializes arrays
+ call initialize()
+
+ ! reads in parameters needed
+ call read_parameters()
+
+ ! reads in current transverse isotropic model files: vpv.. & vsv.. & eta & rho
+ call read_model()
+
+ ! reads in smoothed kernels: bulk, betav, betah, eta
+ call read_kernels()
+
+ ! computes volume element associated with points, calculates kernel integral for statistics
+ call compute_volume()
+
+ ! calculates gradient
+ ! steepest descent method
+ call get_gradient()
+
+ ! compute new model in terms of alpha, beta, eta and rho
+ ! (see also Carl's Latex notes)
+
+ ! model update:
+ ! transverse isotropic update only in layer Moho to 220 (where SPECFEM3D_GLOBE considers TISO)
+ ! everywhere else uses an isotropic update
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ ! initial model values
+ eta0 = model_eta(i,j,k,ispec)
+ betav0 = model_vsv(i,j,k,ispec)
+ betah0 = model_vsh(i,j,k,ispec)
+ rho0 = model_rho(i,j,k,ispec)
+ alphav0 = model_vpv(i,j,k,ispec)
+ alphah0 = model_vph(i,j,k,ispec)
+
+ eta1 = 0._CUSTOM_REAL
+ betav1 = 0._CUSTOM_REAL
+ betah1 = 0._CUSTOM_REAL
+ rho1 = 0._CUSTOM_REAL
+ alphav1 = 0._CUSTOM_REAL
+ alphah1 = 0._CUSTOM_REAL
+
+ ! do not use transverse isotropy except if element is between d220 and Moho
+ if(.not. ( idoubling(ispec)==IFLAG_220_80 .or. idoubling(ispec)==IFLAG_80_MOHO) ) then
+
+ ! isotropic model update
+
+ ! no eta perturbation, since eta = 1 in isotropic media
+ eta1 = eta0
+
+ ! shear values
+ ! isotropic kernel K_beta = K_betav + K_betah
+ ! with same scaling step_length the model update dbeta_iso = dbetav + dbetah
+ ! note:
+ ! this step length can be twice as big as that given by the input
+ dbetaiso = model_dbetav(i,j,k,ispec) + model_dbetah(i,j,k,ispec)
+ betav1 = betav0 * exp( dbetaiso )
+ betah1 = betah0 * exp( dbetaiso )
+ ! note: betah is probably not really used in isotropic layers
+ ! (see SPECFEM3D_GLOBE/get_model.f90)
+
+ ! density: uses scaling relation with isotropic shear perturbations
+ ! dln rho = RHO_SCALING * dln betaiso
+ rho1 = rho0 * exp( RHO_SCALING * dbetaiso )
+
+ ! alpha values
+ dbulk = model_dbulk(i,j,k,ispec)
+ alphav1 = sqrt( alphav0**2 * exp(2.0*dbulk) + FOUR_THIRDS * betav0**2 * ( &
+ exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
+ alphah1 = sqrt( alphah0**2 * exp(2.0*dbulk) + FOUR_THIRDS * betah0**2 * ( &
+ exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
+ ! note: alphah probably not used in SPECFEM3D_GLOBE
+
+ else
+
+ ! transverse isotropic model update
+
+ ! eta value : limits updated values for eta range constraint
+ eta1 = eta0 * exp( model_deta(i,j,k,ispec) )
+ if( eta1 < LIMIT_ETA_MIN ) eta1 = LIMIT_ETA_MIN
+ if( eta1 > LIMIT_ETA_MAX ) eta1 = LIMIT_ETA_MAX
+
+ ! shear values
+ betav1 = betav0 * exp( model_dbetav(i,j,k,ispec) )
+ betah1 = betah0 * exp( model_dbetah(i,j,k,ispec) )
+
+ ! density: uses scaling relation with Voigt average of shear perturbations
+ betaiso0 = sqrt( ( 2.0 * betav0**2 + betah0**2 ) / 3.0 )
+ betaiso1 = sqrt( ( 2.0 * betav1**2 + betah1**2 ) / 3.0 )
+ dbetaiso = log( betaiso1 / betaiso0 )
+ rho1 = rho0 * exp( RHO_SCALING * dbetaiso )
+
+ ! alpha values
+ dbulk = model_dbulk(i,j,k,ispec)
+ alphav1 = sqrt( alphav0**2 * exp(2.0*dbulk) &
+ + FOUR_THIRDS * betav0**2 * ( &
+ exp(2.0*model_dbetav(i,j,k,ispec)) - exp(2.0*dbulk) ) )
+ alphah1 = sqrt( alphah0**2 * exp(2.0*dbulk) &
+ + FOUR_THIRDS * betah0**2 * ( &
+ exp(2.0*model_dbetah(i,j,k,ispec)) - exp(2.0*dbulk) ) )
+
+ endif
+
+
+ ! stores new model values
+ model_vpv_new(i,j,k,ispec) = alphav1
+ model_vph_new(i,j,k,ispec) = alphah1
+ model_vsv_new(i,j,k,ispec) = betav1
+ model_vsh_new(i,j,k,ispec) = betah1
+ model_eta_new(i,j,k,ispec) = eta1
+ model_rho_new(i,j,k,ispec) = rho1
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! stores new model in files
+ call store_new_model()
+
+ ! stores relative model perturbations
+ call store_perturbations()
+
+ ! stop all the MPI processes, and exit
+ call MPI_FINALIZE(ier)
+
+end program add_model
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine initialize()
+
+! initializes arrays
+
+ use model_update
+ implicit none
+
+ ! initialize the MPI communicator and start the NPROCTOT MPI processes
+ call MPI_INIT(ier)
+ call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
+ call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
+
+ if( sizeprocs /= nchunks_val*nproc_xi_val*nproc_eta_val ) then
+ print*,'sizeprocs:',sizeprocs,nchunks_val,nproc_xi_val,nproc_eta_val
+ call exit_mpi('error number sizeprocs')
+ endif
+
+ ! model
+ model_vpv = 0.0_CUSTOM_REAL
+ model_vph = 0.0_CUSTOM_REAL
+ model_vsv = 0.0_CUSTOM_REAL
+ model_vsh = 0.0_CUSTOM_REAL
+ model_eta = 0.0_CUSTOM_REAL
+ model_rho = 0.0_CUSTOM_REAL
+
+ model_vpv_new = 0.0_CUSTOM_REAL
+ model_vph_new = 0.0_CUSTOM_REAL
+ model_vsv_new = 0.0_CUSTOM_REAL
+ model_vsh_new = 0.0_CUSTOM_REAL
+ model_eta_new = 0.0_CUSTOM_REAL
+ model_rho_new = 0.0_CUSTOM_REAL
+
+ ! model updates
+ model_dbulk = 0.0_CUSTOM_REAL
+ model_dbetah = 0.0_CUSTOM_REAL
+ model_dbetav = 0.0_CUSTOM_REAL
+ model_deta = 0.0_CUSTOM_REAL
+
+ ! gradients
+ kernel_bulk = 0.0_CUSTOM_REAL
+ kernel_betav = 0.0_CUSTOM_REAL
+ kernel_betah = 0.0_CUSTOM_REAL
+ kernel_eta = 0.0_CUSTOM_REAL
+
+end subroutine initialize
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine read_parameters()
+
+! reads in parameters needed
+
+ use model_update
+ implicit none
+ character(len=150) :: s_step_fac
+
+ ! subjective step length to multiply to the gradient
+ !step_fac = 0.03
+
+ call getarg(1,s_step_fac)
+
+ if (trim(s_step_fac) == '') then
+ call exit_MPI(myrank,'Usage: add_model_globe_tiso step_factor')
+ endif
+
+ ! read in parameter information
+ read(s_step_fac,*) step_fac
+ !if( abs(step_fac) < 1.e-10) then
+ ! print*,'error: step factor ',step_fac
+ ! call exit_MPI(myrank,'error step factor')
+ !endif
+
+ if (myrank == 0) then
+ print*,'defaults'
+ print*,' NPROC_XI , NPROC_ETA: ',nproc_xi_val,nproc_eta_val
+ print*,' NCHUNKS: ',nchunks_val
+ print*
+ print*,'model update for vsv,vsh,vpv,vph,eta,rho:'
+ print*,' step_fac = ',step_fac
+ print*
+
+ endif
+
+
+end subroutine read_parameters
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine read_model()
+
+! reads in current transverse isotropic model: vpv.. & vsv.. & eta & rho
+
+ use model_update
+ implicit none
+
+ ! vpv model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vpv.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_vpv(:,:,:,1:nspec)
+ close(12)
+
+ ! vph model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vph.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_vph(:,:,:,1:nspec)
+ close(12)
+
+ ! vsv model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vsv.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_vsv(:,:,:,1:nspec)
+ close(12)
+
+ ! vsh model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vsh.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_vsh(:,:,:,1:nspec)
+ close(12)
+
+ ! eta model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_eta.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_eta(:,:,:,1:nspec)
+ close(12)
+
+ ! rho model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_rho.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_rho(:,:,:,1:nspec)
+ close(12)
+
+ ! statistics
+ call mpi_reduce(minval(model_vpv),min_vpv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_vpv),max_vpv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_vph),min_vph,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_vph),max_vph,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_vsv),min_vsv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_vsv),max_vsv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_vsh),min_vsh,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_vsh),max_vsh,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_eta),min_eta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_eta),max_eta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_rho),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_rho),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'initial models:'
+ print*,' vpv min/max: ',min_vpv,max_vpv
+ print*,' vph min/max: ',min_vph,max_vph
+ print*,' vsv min/max: ',min_vsv,max_vsv
+ print*,' vsh min/max: ',min_vsh,max_vsh
+ print*,' eta min/max: ',min_eta,max_eta
+ print*,' rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+end subroutine read_model
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine read_kernels()
+
+! reads in smoothed kernels: bulk, betav, betah, eta
+
+ use model_update
+ implicit none
+
+ ! bulk kernel
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_bulk_c_kernel_smooth.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) kernel_bulk(:,:,:,1:nspec)
+ close(12)
+
+ ! betav kernel
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_bulk_betav_kernel_smooth.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) kernel_betav(:,:,:,1:nspec)
+ close(12)
+
+ ! betah kernel
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_bulk_betah_kernel_smooth.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) kernel_betah(:,:,:,1:nspec)
+ close(12)
+
+ ! eta kernel
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_eta_kernel_smooth.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) kernel_eta(:,:,:,1:nspec)
+ close(12)
+
+
+ ! statistics
+ call mpi_reduce(minval(kernel_bulk),min_bulk,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(kernel_bulk),max_bulk,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(kernel_betah),min_vsh,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(kernel_betah),max_vsh,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(kernel_betav),min_vsv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(kernel_betav),max_vsv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(kernel_eta),min_eta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(kernel_eta),max_eta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'initial kernels:'
+ print*,' bulk min/max : ',min_bulk,max_bulk
+ print*,' betav min/max: ',min_vsv,max_vsv
+ print*,' betah min/max: ',min_vsh,max_vsh
+ print*,' eta min/max : ',min_eta,max_eta
+ print*
+ endif
+
+end subroutine read_kernels
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine compute_volume()
+
+! computes volume element associated with points
+
+ use model_update
+ implicit none
+ ! jacobian
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: jacobian
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl, &
+ jacobianl,volumel
+ ! integration values
+ real(kind=CUSTOM_REAL) :: integral_bulk_sum,integral_betav_sum, &
+ integral_betah_sum,integral_eta_sum
+ real(kind=CUSTOM_REAL) :: kernel_integral_bulk,kernel_integral_betav,&
+ kernel_integral_betah,kernel_integral_eta
+
+ ! GLL points
+ wgll_cube = 0.0
+ call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ wgll_cube(i,j,k) = wxgll(i)*wygll(j)*wzgll(k)
+ enddo
+ enddo
+ enddo
+
+ ! global addressing
+ write(m_file,'(a,i6.6,a)') 'topo/proc',myrank,'_reg1_solver_data_2.bin'
+ open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(11) x(1:nglob)
+ read(11) y(1:nglob)
+ read(11) z(1:nglob)
+ read(11) ibool(:,:,:,1:nspec)
+ read(11) idoubling(1:nspec)
+ close(11)
+
+ ! builds jacobian
+ write(m_file,'(a,i6.6,a)') 'topo/proc',myrank,'_reg1_solver_data_1.bin'
+ open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(11) xix
+ read(11) xiy
+ read(11) xiz
+ read(11) etax
+ read(11) etay
+ read(11) etaz
+ read(11) gammax
+ read(11) gammay
+ read(11) gammaz
+ close(11)
+
+ jacobian = 0.0
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ ! gets derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ ! computes the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+ jacobian(i,j,k,ispec) = jacobianl
+
+ !if( abs(jacobianl) < 1.e-8 ) then
+ ! print*,'rank ',myrank,'jacobian: ',jacobianl,i,j,k,wgll_cube(i,j,k)
+ !endif
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! volume associated with global point
+ volume_glob = 0.0
+ kernel_integral_bulk = 0._CUSTOM_REAL
+ kernel_integral_betav = 0._CUSTOM_REAL
+ kernel_integral_betah = 0._CUSTOM_REAL
+ kernel_integral_eta = 0._CUSTOM_REAL
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+ if( iglob == 0 ) then
+ print*,'iglob zero',i,j,k,ispec
+ print*
+ print*,'ibool:',ispec
+ print*,ibool(:,:,:,ispec)
+ print*
+ call exit_MPI('error ibool')
+ endif
+
+ ! volume associated with GLL point
+ volumel = jacobian(i,j,k,ispec)*wgll_cube(i,j,k)
+
+ ! kernel integration: for each element
+ kernel_integral_bulk = kernel_integral_bulk &
+ + volumel * kernel_bulk(i,j,k,ispec)
+
+ kernel_integral_betav = kernel_integral_betav &
+ + volumel * kernel_betav(i,j,k,ispec)
+
+ kernel_integral_betah = kernel_integral_betah &
+ + volumel * kernel_betah(i,j,k,ispec)
+
+ kernel_integral_eta = kernel_integral_eta &
+ + volumel * kernel_eta(i,j,k,ispec)
+
+ ! gradient vector norm sqrt( v^T * v )
+ norm_bulk = norm_bulk + kernel_bulk(i,j,k,ispec)**2
+ norm_betav = norm_betav + kernel_betav(i,j,k,ispec)**2
+ norm_betah = norm_betah + kernel_betah(i,j,k,ispec)**2
+ norm_eta = norm_eta + kernel_eta(i,j,k,ispec)**2
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! statistics
+ ! kernel integration: for whole volume
+ call mpi_reduce(kernel_integral_bulk,integral_bulk_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(kernel_integral_betav,integral_betav_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(kernel_integral_betah,integral_betah_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(kernel_integral_eta,integral_eta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'integral kernels:'
+ print*,' bulk : ',integral_bulk_sum
+ print*,' betav : ',integral_betav_sum
+ print*,' betah : ',integral_betah_sum
+ print*,' eta : ',integral_eta_sum
+ print*
+ endif
+
+ ! norms: for whole volume
+ call mpi_reduce(norm_bulk,norm_bulk_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_betav,norm_betav_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_betah,norm_betah_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_eta,norm_eta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+
+ norm_bulk = sqrt(norm_bulk_sum)
+ norm_betav = sqrt(norm_betav_sum)
+ norm_betah = sqrt(norm_betah_sum)
+ norm_eta = sqrt(norm_eta_sum)
+
+
+ if( myrank == 0 ) then
+ print*,'norm kernels:'
+ print*,' bulk : ',norm_bulk
+ print*,' betav : ',norm_betav
+ print*,' betah : ',norm_betah
+ print*,' eta : ',norm_eta
+ print*
+ endif
+
+end subroutine compute_volume
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine get_gradient()
+
+! calculates gradient by steepest descent method
+
+ use model_update
+ implicit none
+
+ ! gradient in negative direction for steepest descent
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ ! for bulk
+ model_dbulk(i,j,k,ispec) = - kernel_bulk(i,j,k,ispec)
+
+ ! for shear
+ model_dbetav(i,j,k,ispec) = - kernel_betav(i,j,k,ispec)
+ model_dbetah(i,j,k,ispec) = - kernel_betah(i,j,k,ispec)
+
+ ! for eta
+ model_deta(i,j,k,ispec) = - kernel_eta(i,j,k,ispec)
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! stores model_dbulk, ... arrays
+ call store_kernel_updates()
+
+ ! statistics
+ call mpi_reduce(minval(model_dbulk),min_bulk,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dbulk),max_bulk,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_dbetav),min_vsv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dbetav),max_vsv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_dbetah),min_vsh,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dbetah),max_vsh,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_deta),min_eta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_deta),max_eta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'initial gradients:'
+ print*,' bulk min/max : ',min_bulk,max_bulk
+ print*,' betav min/max: ',min_vsv,max_vsv
+ print*,' betah min/max: ',min_vsh,max_vsh
+ print*,' eta min/max : ',min_eta,max_eta
+ print*
+ endif
+
+ ! determines step length based on maximum gradient value (either vsv or vsh)
+ if( myrank == 0 ) then
+ ! maximum gradient values
+ minmax(1) = abs(min_vsv)
+ minmax(2) = abs(max_vsv)
+ minmax(3) = abs(min_vsh)
+ minmax(4) = abs(max_vsh)
+
+ ! chooses step length such that it becomes the desired, given step factor as inputted
+ max = maxval(minmax)
+ step_length = step_fac/max
+
+ print*,' step length : ',step_length,max
+ print*
+ endif
+ call mpi_bcast(step_length,1,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+
+
+ ! gradient length sqrt( v^T * v )
+ norm_bulk = sum( model_dbulk * model_dbulk )
+ norm_betav = sum( model_dbetav * model_dbetav )
+ norm_betah = sum( model_dbetah * model_dbetah )
+ norm_eta = sum( model_deta * model_deta )
+
+ call mpi_reduce(norm_bulk,norm_bulk_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_betav,norm_betav_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_betah,norm_betah_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_eta,norm_eta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+
+ norm_bulk = sqrt(norm_bulk_sum)
+ norm_betav = sqrt(norm_betav_sum)
+ norm_betah = sqrt(norm_betah_sum)
+ norm_eta = sqrt(norm_eta_sum)
+
+ if( myrank == 0 ) then
+ print*,'norm model updates:'
+ print*,' bulk : ',norm_bulk
+ print*,' betav: ',norm_betav
+ print*,' betah: ',norm_betah
+ print*,' eta : ',norm_eta
+ print*
+ endif
+
+ ! multiply model updates by a subjective factor that will change the step
+ model_dbulk(:,:,:,:) = step_length * model_dbulk(:,:,:,:)
+ model_dbetav(:,:,:,:) = step_length * model_dbetav(:,:,:,:)
+ model_dbetah(:,:,:,:) = step_length * model_dbetah(:,:,:,:)
+ model_deta(:,:,:,:) = step_length * model_deta(:,:,:,:)
+
+
+ ! statistics
+ call mpi_reduce(minval(model_dbulk),min_bulk,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dbulk),max_bulk,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_dbetav),min_vsv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dbetav),max_vsv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_dbetah),min_vsh,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dbetah),max_vsh,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_deta),min_eta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_deta),max_eta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'scaled gradients:'
+ print*,' bulk min/max : ',min_bulk,max_bulk
+ print*,' betav min/max: ',min_vsv,max_vsv
+ print*,' betah min/max: ',min_vsh,max_vsh
+ print*,' eta min/max : ',min_eta,max_eta
+ print*
+ endif
+
+end subroutine get_gradient
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine store_kernel_updates()
+
+! file output for new model
+
+ use model_update
+ implicit none
+
+ ! kernel updates
+ fname = 'dbulk_c'
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_dbulk
+ close(12)
+
+ fname = 'dbetav'
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_dbetav
+ close(12)
+
+ fname = 'dbetah'
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_dbetah
+ close(12)
+
+ fname = 'deta'
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_deta
+ close(12)
+
+end subroutine store_kernel_updates
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine store_new_model()
+
+! file output for new model
+
+ use model_update
+ implicit none
+
+ ! vpv model
+ call mpi_reduce(maxval(model_vpv_new),max_vpv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_vpv_new),min_vpv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'vpv_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_vpv_new
+ close(12)
+
+ ! vph model
+ call mpi_reduce(maxval(model_vph_new),max_vph,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_vph_new),min_vph,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'vph_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_vph_new
+ close(12)
+
+ ! vsv model
+ call mpi_reduce(maxval(model_vsv_new),max_vsv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_vsv_new),min_vsv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'vsv_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_vsv_new
+ close(12)
+
+ ! vsh model
+ call mpi_reduce(maxval(model_vsh_new),max_vsh,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_vsh_new),min_vsh,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'vsh_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_vsh_new
+ close(12)
+
+ ! eta model
+ call mpi_reduce(maxval(model_eta_new),max_eta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_eta_new),min_eta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'eta_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_eta_new
+ close(12)
+
+ ! rho model
+ call mpi_reduce(maxval(model_rho_new),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_rho_new),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'rho_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_rho_new
+ close(12)
+
+
+ if( myrank == 0 ) then
+ print*,'new models:'
+ print*,' vpv min/max: ',min_vpv,max_vpv
+ print*,' vph min/max: ',min_vph,max_vph
+ print*,' vsv min/max: ',min_vsv,max_vsv
+ print*,' vsh min/max: ',min_vsh,max_vsh
+ print*,' eta min/max: ',min_eta,max_eta
+ print*,' rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+
+end subroutine store_new_model
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine store_perturbations()
+
+! file output for new model
+
+ use model_update
+ implicit none
+
+ ! vpv relative perturbations
+ ! logarithmic perturbation: log( v_new) - log( v_old) = log( v_new / v_old )
+ total_model = 0.0_CUSTOM_REAL
+ where( model_vpv /= 0.0 ) total_model = log( model_vpv_new / model_vpv)
+ ! or
+ ! linear approximation: (v_new - v_old) / v_old
+ !where( model_vpv /= 0.0 ) total_model = ( model_vpv_new - model_vpv) / model_vpv
+
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvpvvpv.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_vpv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vpv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! vph relative perturbations
+ total_model = 0.0_CUSTOM_REAL
+ where( model_vph /= 0.0 ) total_model = log( model_vph_new / model_vph)
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvphvph.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_vph,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vph,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! vsv relative perturbations
+ total_model = 0.0_CUSTOM_REAL
+ where( model_vsv /= 0.0 ) total_model = log( model_vsv_new / model_vsv)
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvsvvsv.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_vsv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vsv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! vsh relative perturbations
+ total_model = 0.0_CUSTOM_REAL
+ where( model_vsh /= 0.0 ) total_model = log( model_vsh_new / model_vsh)
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvshvsh.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_vsh,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vsh,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! eta relative perturbations
+ total_model = 0.0_CUSTOM_REAL
+ where( model_eta /= 0.0 ) total_model = log( model_eta_new / model_eta)
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_detaeta.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_eta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_eta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! rho relative model perturbations
+ total_model = 0.0_CUSTOM_REAL
+ where( model_rho /= 0.0 ) total_model = log( model_rho_new / model_rho)
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_drhorho.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'relative update:'
+ print*,' dvpv/vpv min/max: ',min_vpv,max_vpv
+ print*,' dvph/vph min/max: ',min_vph,max_vph
+ print*,' dvsv/vsv min/max: ',min_vsv,max_vsv
+ print*,' dvsh/vsh min/max: ',min_vsh,max_vsh
+ print*,' deta/eta min/max: ',min_eta,max_eta
+ print*,' drho/rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+end subroutine store_perturbations
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso_iso.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso_iso.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso_iso.f90 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,1005 @@
+! add_model_globe_tiso_iso
+!
+! this program can be used to update TRANSVERSE ISOTROPIC model files
+! based on smoothed, summed, ISOTROPIC event kernels.
+! the kernels are given for isotropic parameters (bulk_c,bulk_beta,rho) or (alpha,beta,rho).
+!
+! the algorithm uses a steepest descent method with a step length
+! determined by the given maximum update percentage.
+!
+! input:
+! - step_fac : step length to update the models, f.e. 0.03 for plusminus 3%
+!
+! setup:
+!
+!- INPUT_MODEL/ contains:
+! proc000***_reg1_vsv.bin &
+! proc000***_reg1_vsh.bin &
+! proc000***_reg1_vpv.bin &
+! proc000***_reg1_vph.bin &
+! proc000***_reg1_eta.bin &
+! proc000***_reg1_rho.bin
+!
+!- INPUT_GRADIENT/ contains:
+! proc000***_reg1_bulk_c_kernel_smooth.bin &
+! proc000***_reg1_bulk_beta_kernel_smooth.bin &
+! proc000***_reg1_rho_kernel_smooth.bin
+! or
+! proc000***_reg1_alpha_kernel_smooth.bin &
+! proc000***_reg1_beta_kernel_smooth.bin &
+! proc000***_reg1_rho_kernel_smooth.bin
+!
+!- topo/ contains:
+! proc000***_reg1_solver_data_1.bin
+!
+! new models are stored in
+!- OUTPUT_MODEL/ as
+! proc000***_reg1_vpv_new.bin &
+! proc000***_reg1_vph_new.bin &
+! proc000***_reg1_vsv_new.bin &
+! proc000***_reg1_vsh_new.bin &
+! proc000***_reg1_eta_new.bin &
+! proc000***_reg1_rho_new.bin
+!
+! USAGE: ./add_model_globe_tiso_iso 0.3
+
+module model_update_iso
+
+ include 'mpif.h'
+ include 'constants_globe.h'
+ include 'precision_globe.h'
+ include 'values_from_mesher_globe.h'
+
+ ! ======================================================
+
+ ! USER PARAMETERS
+
+ ! by default, this algorithm uses (bulk_c,bulk_beta,rho) kernels to update vpv,vph,vsv,vsh,rho,eta
+ ! if you prefer using (alpha,beta,rho) kernels, set this flag to true
+ logical, parameter :: USE_ALPHA_BETA_RHO = .false.
+
+ ! ignore rho kernel, but use density perturbations as a scaling of Vs perturbations
+ logical, parameter :: USE_RHO_SCALING = .false.
+
+ ! in case of rho scaling, specifies density scaling factor with shear perturbations
+ ! see e.g. Montagner & Anderson (1989), Panning & Romanowicz (2006)
+ real(kind=CUSTOM_REAL),parameter :: RHO_SCALING = 0.33_CUSTOM_REAL
+
+ ! ======================================================
+
+ integer, parameter :: NSPEC = NSPEC_CRUST_MANTLE
+ integer, parameter :: NGLOB = NGLOB_CRUST_MANTLE
+
+ ! transverse isotropic model files
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ model_vpv,model_vph,model_vsv,model_vsh,model_eta,model_rho
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ model_vpv_new,model_vph_new,model_vsv_new,model_vsh_new,model_eta_new,model_rho_new
+
+ ! model updates
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ model_dA,model_dB,model_dR
+
+ ! kernels
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ kernel_a,kernel_b,kernel_rho
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: total_model
+
+ ! volume
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: x, y, z
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: volume_spec
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: volume_glob,volume_glob_sum
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
+ integer, dimension(NSPEC) :: idoubling
+
+ ! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: xigll, wxgll
+ double precision, dimension(NGLLY) :: yigll, wygll
+ double precision, dimension(NGLLZ) :: zigll, wzgll
+
+ ! array with all the weights in the cube
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+
+ ! gradient vector norm ( v^T * v )
+ real(kind=CUSTOM_REAL) :: norm_a,norm_beta,norm_rho
+ real(kind=CUSTOM_REAL) :: norm_a_sum,norm_beta_sum,norm_rho_sum
+
+ ! model update length
+ real(kind=CUSTOM_REAL) :: step_fac,step_length
+
+ real(kind=CUSTOM_REAL) :: min_vpv,min_vph,min_vsv,min_vsh, &
+ max_vpv,max_vph,max_vsv,max_vsh,min_b,max_b,min_a,max_a, &
+ min_rho,max_rho,min_eta,max_eta,max,minmax(4)
+
+ real(kind=CUSTOM_REAL) :: betav1,betah1,betav0,betah0,rho1,rho0, &
+ betaiso1,betaiso0,eta1,eta0,alphav1,alphav0,alphah1,alphah0
+ real(kind=CUSTOM_REAL) :: dbetaiso,dbulk
+
+ integer :: nfile, myrank, sizeprocs, ier
+ integer :: i, j, k,ispec, iglob, ishell, n, it, j1, ib, npts_sem, ios
+ character(len=150) :: sline, m_file, fname
+
+
+end module model_update_iso
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+program add_model
+
+ use model_update_iso
+
+ implicit none
+
+ ! ============ program starts here =====================
+
+ ! initializes arrays
+ call initialize()
+
+ ! reads in parameters needed
+ call read_parameters()
+
+ ! reads in current transverse isotropic model files: vpv.. & vsv.. & eta & rho
+ call read_model()
+
+ ! reads in smoothed kernels: bulk, beta, rho
+ call read_kernels()
+
+ ! computes volume element associated with points, calculates kernel integral for statistics
+ call compute_volume()
+
+ ! calculates gradient
+ ! steepest descent method
+ call get_gradient()
+
+ ! compute new model in terms of alpha, beta, eta and rho
+ ! (see also Carl's Latex notes)
+
+ ! model update:
+ ! transverse isotropic update only in layer Moho to 220 (where SPECFEM3D_GLOBE considers TISO)
+ ! everywhere else uses an isotropic update
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ ! initial model values
+ eta0 = model_eta(i,j,k,ispec)
+ betav0 = model_vsv(i,j,k,ispec)
+ betah0 = model_vsh(i,j,k,ispec)
+ rho0 = model_rho(i,j,k,ispec)
+ alphav0 = model_vpv(i,j,k,ispec)
+ alphah0 = model_vph(i,j,k,ispec)
+
+ eta1 = 0._CUSTOM_REAL
+ betav1 = 0._CUSTOM_REAL
+ betah1 = 0._CUSTOM_REAL
+ rho1 = 0._CUSTOM_REAL
+ alphav1 = 0._CUSTOM_REAL
+ alphah1 = 0._CUSTOM_REAL
+
+ ! isotropic model update
+
+ ! no eta perturbation for isotropic updates, takes the same as initial values
+ eta1 = eta0
+
+ ! shear values
+ dbetaiso = model_dB(i,j,k,ispec)
+ betav1 = betav0 * exp( dbetaiso )
+ betah1 = betah0 * exp( dbetaiso )
+ ! note: betah is probably not really used in isotropic layers
+ ! (see SPECFEM3D_GLOBE/get_model.f90)
+
+ ! density
+ rho1 = rho0 * exp( model_dR(i,j,k,ispec) )
+
+ ! alpha values
+ dbulk = model_dA(i,j,k,ispec)
+ alphav1 = sqrt( alphav0**2 * exp(2.0*dbulk) + FOUR_THIRDS * betav0**2 * ( &
+ exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
+ alphah1 = sqrt( alphah0**2 * exp(2.0*dbulk) + FOUR_THIRDS * betah0**2 * ( &
+ exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
+ ! note: alphah probably not used in SPECFEM3D_GLOBE
+
+ ! stores new model values
+ model_vpv_new(i,j,k,ispec) = alphav1
+ model_vph_new(i,j,k,ispec) = alphah1
+ model_vsv_new(i,j,k,ispec) = betav1
+ model_vsh_new(i,j,k,ispec) = betah1
+ model_eta_new(i,j,k,ispec) = eta1
+ model_rho_new(i,j,k,ispec) = rho1
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! stores new model in files
+ call store_new_model()
+
+ ! stores relative model perturbations
+ call store_perturbations()
+
+ ! stop all the MPI processes, and exit
+ call MPI_FINALIZE(ier)
+
+end program add_model
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine initialize()
+
+! initializes arrays
+
+ use model_update_iso
+ implicit none
+
+ ! initialize the MPI communicator and start the NPROCTOT MPI processes
+ call MPI_INIT(ier)
+ call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
+ call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
+
+ if( sizeprocs /= nchunks_val*nproc_xi_val*nproc_eta_val ) then
+ print*,'sizeprocs:',sizeprocs,nchunks_val,nproc_xi_val,nproc_eta_val
+ call exit_mpi('error number sizeprocs')
+ endif
+
+ ! model
+ model_vpv = 0.0_CUSTOM_REAL
+ model_vph = 0.0_CUSTOM_REAL
+ model_vsv = 0.0_CUSTOM_REAL
+ model_vsh = 0.0_CUSTOM_REAL
+ model_eta = 0.0_CUSTOM_REAL
+ model_rho = 0.0_CUSTOM_REAL
+
+ model_vpv_new = 0.0_CUSTOM_REAL
+ model_vph_new = 0.0_CUSTOM_REAL
+ model_vsv_new = 0.0_CUSTOM_REAL
+ model_vsh_new = 0.0_CUSTOM_REAL
+ model_eta_new = 0.0_CUSTOM_REAL
+ model_rho_new = 0.0_CUSTOM_REAL
+
+ ! model updates
+ model_dA = 0.0_CUSTOM_REAL
+ model_dB = 0.0_CUSTOM_REAL
+ model_dR = 0.0_CUSTOM_REAL
+
+ ! gradients
+ kernel_a = 0.0_CUSTOM_REAL
+ kernel_b = 0.0_CUSTOM_REAL
+ kernel_rho = 0.0_CUSTOM_REAL
+
+end subroutine initialize
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine read_parameters()
+
+! reads in parameters needed
+
+ use model_update_iso
+ implicit none
+ character(len=150) :: s_step_fac
+
+ ! subjective step length to multiply to the gradient
+ !step_fac = 0.03
+
+ call getarg(1,s_step_fac)
+
+ if (trim(s_step_fac) == '') then
+ call exit_MPI(myrank,'Usage: add_model_globe_tiso step_factor')
+ endif
+
+ ! read in parameter information
+ read(s_step_fac,*) step_fac
+ !if( abs(step_fac) < 1.e-10) then
+ ! print*,'error: step factor ',step_fac
+ ! call exit_MPI(myrank,'error step factor')
+ !endif
+
+ if (myrank == 0) then
+ print*,'defaults'
+ print*,' NPROC_XI , NPROC_ETA: ',nproc_xi_val,nproc_eta_val
+ print*,' NCHUNKS: ',nchunks_val
+ print*
+ print*,'model update for vsv,vsh,vpv,vph,eta,rho:'
+ print*,' step_fac = ',step_fac
+ print*
+ if( USE_ALPHA_BETA_RHO ) then
+ print*,'kernel parameterization: (alpha,beta,rho)'
+ else
+ print*,'kernel parameterization: (bulk,beta,rho)'
+ endif
+ print*
+ if( USE_RHO_SCALING ) then
+ print*,'scaling rho perturbations'
+ print*
+ endif
+
+ endif
+
+
+end subroutine read_parameters
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine read_model()
+
+! reads in current transverse isotropic model: vpv.. & vsv.. & eta & rho
+
+ use model_update_iso
+ implicit none
+
+ ! vpv model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vpv.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_vpv(:,:,:,1:nspec)
+ close(12)
+
+ ! vph model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vph.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_vph(:,:,:,1:nspec)
+ close(12)
+
+ ! vsv model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vsv.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_vsv(:,:,:,1:nspec)
+ close(12)
+
+ ! vsh model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vsh.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_vsh(:,:,:,1:nspec)
+ close(12)
+
+ ! eta model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_eta.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_eta(:,:,:,1:nspec)
+ close(12)
+
+ ! rho model
+ write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_rho.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) model_rho(:,:,:,1:nspec)
+ close(12)
+
+ ! statistics
+ call mpi_reduce(minval(model_vpv),min_vpv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_vpv),max_vpv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_vph),min_vph,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_vph),max_vph,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_vsv),min_vsv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_vsv),max_vsv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_vsh),min_vsh,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_vsh),max_vsh,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_eta),min_eta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_eta),max_eta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_rho),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_rho),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'initial models:'
+ print*,' vpv min/max: ',min_vpv,max_vpv
+ print*,' vph min/max: ',min_vph,max_vph
+ print*,' vsv min/max: ',min_vsv,max_vsv
+ print*,' vsh min/max: ',min_vsh,max_vsh
+ print*,' eta min/max: ',min_eta,max_eta
+ print*,' rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+end subroutine read_model
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine read_kernels()
+
+! reads in smoothed kernels: bulk, betav, betah, eta
+
+ use model_update_iso
+ implicit none
+
+ ! bulk kernel
+ if( USE_ALPHA_BETA_RHO ) then
+ ! reads in alpha kernel
+ fname = 'alpha_kernel_smooth'
+ else
+ ! reads in bulk_c kernel
+ fname = 'bulk_c_kernel_smooth'
+ endif
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) kernel_a(:,:,:,1:nspec)
+ close(12)
+
+ ! shear kernel
+ if( USE_ALPHA_BETA_RHO ) then
+ ! reads in beta kernel
+ fname = 'beta_kernel_smooth'
+ else
+ ! reads in bulk_beta kernel
+ fname = 'bulk_beta_kernel_smooth'
+ endif
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) kernel_b(:,:,:,1:nspec)
+ close(12)
+
+ ! rho kernel
+ if( USE_RHO_SCALING ) then
+
+ ! uses scaling relation with shear perturbations
+ kernel_rho(:,:,:,:) = RHO_SCALING * kernel_b(:,:,:,:)
+
+ else
+
+ ! uses rho kernel
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_rho_kernel_smooth.bin'
+ open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(12) kernel_rho(:,:,:,1:nspec)
+ close(12)
+ endif
+
+
+ ! statistics
+ call mpi_reduce(minval(kernel_a),min_a,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(kernel_a),max_a,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(kernel_b),min_b,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(kernel_b),max_b,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(kernel_rho),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(kernel_rho),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'initial kernels:'
+ print*,' a min/max : ',min_a,max_a
+ print*,' beta min/max: ',min_b,max_b
+ print*,' rho min/max : ',min_rho,max_rho
+ print*
+ endif
+
+end subroutine read_kernels
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine compute_volume()
+
+! computes volume element associated with points
+
+ use model_update_iso
+ implicit none
+ ! jacobian
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: jacobian
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl, &
+ jacobianl,volumel
+ ! integration values
+ real(kind=CUSTOM_REAL) :: integral_a_sum,integral_b_sum,integral_rho_sum
+ real(kind=CUSTOM_REAL) :: kernel_integral_a,kernel_integral_b,kernel_integral_rho
+
+ ! GLL points
+ wgll_cube = 0.0
+ call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ wgll_cube(i,j,k) = wxgll(i)*wygll(j)*wzgll(k)
+ enddo
+ enddo
+ enddo
+
+ ! global addressing
+ write(m_file,'(a,i6.6,a)') 'topo/proc',myrank,'_reg1_solver_data_2.bin'
+ open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(11) x(1:nglob)
+ read(11) y(1:nglob)
+ read(11) z(1:nglob)
+ read(11) ibool(:,:,:,1:nspec)
+ read(11) idoubling(1:nspec)
+ close(11)
+
+ ! builds jacobian
+ write(m_file,'(a,i6.6,a)') 'topo/proc',myrank,'_reg1_solver_data_1.bin'
+ open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',trim(m_file)
+ call exit_mpi('file not found')
+ endif
+ read(11) xix
+ read(11) xiy
+ read(11) xiz
+ read(11) etax
+ read(11) etay
+ read(11) etaz
+ read(11) gammax
+ read(11) gammay
+ read(11) gammaz
+ close(11)
+
+ jacobian = 0.0
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ ! gets derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ ! computes the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+ jacobian(i,j,k,ispec) = jacobianl
+
+ !if( abs(jacobianl) < 1.e-8 ) then
+ ! print*,'rank ',myrank,'jacobian: ',jacobianl,i,j,k,wgll_cube(i,j,k)
+ !endif
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! volume associated with global point
+ volume_glob = 0.0
+ kernel_integral_a = 0._CUSTOM_REAL
+ kernel_integral_b = 0._CUSTOM_REAL
+ kernel_integral_rho = 0._CUSTOM_REAL
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+ if( iglob == 0 ) then
+ print*,'iglob zero',i,j,k,ispec
+ print*
+ print*,'ibool:',ispec
+ print*,ibool(:,:,:,ispec)
+ print*
+ call exit_MPI('error ibool')
+ endif
+
+ ! volume associated with GLL point
+ volumel = jacobian(i,j,k,ispec)*wgll_cube(i,j,k)
+
+ ! kernel integration: for each element
+ kernel_integral_a = kernel_integral_a &
+ + volumel * kernel_a(i,j,k,ispec)
+
+ kernel_integral_b = kernel_integral_b &
+ + volumel * kernel_b(i,j,k,ispec)
+
+ kernel_integral_rho = kernel_integral_rho &
+ + volumel * kernel_rho(i,j,k,ispec)
+
+ ! gradient vector norm sqrt( v^T * v )
+ norm_a = norm_a + kernel_a(i,j,k,ispec)**2
+ norm_beta = norm_beta + kernel_b(i,j,k,ispec)**2
+ norm_rho = norm_rho + kernel_rho(i,j,k,ispec)**2
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! statistics
+ ! kernel integration: for whole volume
+ call mpi_reduce(kernel_integral_a,integral_a_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(kernel_integral_b,integral_b_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(kernel_integral_rho,integral_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'integral kernels:'
+ print*,' a : ',integral_a_sum
+ print*,' beta : ',integral_b_sum
+ print*,' rho : ',integral_rho_sum
+ print*
+ endif
+
+ ! norms: for whole volume
+ call mpi_reduce(norm_a,norm_a_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_beta,norm_beta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_rho,norm_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+
+ norm_a = sqrt(norm_a_sum)
+ norm_beta = sqrt(norm_beta_sum)
+ norm_rho = sqrt(norm_rho_sum)
+
+
+ if( myrank == 0 ) then
+ print*,'norm kernels:'
+ print*,' a : ',norm_a
+ print*,' beta : ',norm_beta
+ print*,' rho : ',norm_rho
+ print*
+ endif
+
+end subroutine compute_volume
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine get_gradient()
+
+! calculates gradient by steepest descent method
+
+ use model_update_iso
+ implicit none
+
+ ! gradient in negative direction for steepest descent
+ do ispec = 1, NSPEC
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ ! for bulk
+ model_dA(i,j,k,ispec) = - kernel_a(i,j,k,ispec)
+
+ ! for shear
+ model_dB(i,j,k,ispec) = - kernel_b(i,j,k,ispec)
+
+ ! for rho
+ model_dR(i,j,k,ispec) = - kernel_rho(i,j,k,ispec)
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! stores model_dbulk, ... arrays
+ call store_kernel_updates()
+
+ ! statistics
+ call mpi_reduce(minval(model_dA),min_a,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dA),max_a,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_dB),min_b,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dB),max_b,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_dR),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dR),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'initial gradients:'
+ print*,' a min/max : ',min_a,max_a
+ print*,' beta min/max: ',min_b,max_b
+ print*,' rho min/max : ',min_rho,max_rho
+ print*
+ endif
+
+ ! determines step length based on maximum gradient value (either shear or bulk)
+ if( myrank == 0 ) then
+ ! maximum gradient values
+ minmax(1) = abs(min_b)
+ minmax(2) = abs(max_b)
+ minmax(3) = abs(min_a)
+ minmax(4) = abs(max_a)
+
+ ! chooses step length such that it becomes the desired, given step factor as inputted
+ max = maxval(minmax)
+ step_length = step_fac/max
+
+ print*,' step length : ',step_length,max
+ print*
+ endif
+ call mpi_bcast(step_length,1,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+
+
+ ! gradient length sqrt( v^T * v )
+ norm_a = sum( model_dA * model_dA )
+ norm_beta = sum( model_dB * model_dB )
+ norm_rho = sum( model_dR * model_dR )
+
+ call mpi_reduce(norm_a,norm_a_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_beta,norm_beta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(norm_rho,norm_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+
+ norm_a = sqrt(norm_a_sum)
+ norm_beta = sqrt(norm_beta_sum)
+ norm_rho = sqrt(norm_rho_sum)
+
+ if( myrank == 0 ) then
+ print*,'norm model updates:'
+ print*,' a : ',norm_a
+ print*,' beta: ',norm_beta
+ print*,' rho : ',norm_rho
+ print*
+ endif
+
+ ! multiply model updates by a subjective factor that will change the step
+ model_dA(:,:,:,:) = step_length * model_dA(:,:,:,:)
+ model_dB(:,:,:,:) = step_length * model_dB(:,:,:,:)
+ model_dR(:,:,:,:) = step_length * model_dR(:,:,:,:)
+
+
+ ! statistics
+ call mpi_reduce(minval(model_dA),min_a,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dA),max_a,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_dB),min_b,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dB),max_b,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ call mpi_reduce(minval(model_dR),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(maxval(model_dR),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'scaled gradients:'
+ print*,' a min/max : ',min_a,max_a
+ print*,' beta min/max: ',min_b,max_b
+ print*,' rho min/max : ',min_rho,max_rho
+ print*
+ endif
+
+end subroutine get_gradient
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine store_kernel_updates()
+
+! file output for new model
+
+ use model_update_iso
+ implicit none
+
+ ! kernel updates
+ fname = 'dA'
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_dA
+ close(12)
+
+ fname = 'dB'
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_dB
+ close(12)
+
+ fname = 'dR'
+ write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_dR
+ close(12)
+
+end subroutine store_kernel_updates
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine store_new_model()
+
+! file output for new model
+
+ use model_update_iso
+ implicit none
+
+ ! vpv model
+ call mpi_reduce(maxval(model_vpv_new),max_vpv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_vpv_new),min_vpv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'vpv_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_vpv_new
+ close(12)
+
+ ! vph model
+ call mpi_reduce(maxval(model_vph_new),max_vph,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_vph_new),min_vph,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'vph_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_vph_new
+ close(12)
+
+ ! vsv model
+ call mpi_reduce(maxval(model_vsv_new),max_vsv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_vsv_new),min_vsv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'vsv_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_vsv_new
+ close(12)
+
+ ! vsh model
+ call mpi_reduce(maxval(model_vsh_new),max_vsh,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_vsh_new),min_vsh,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'vsh_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_vsh_new
+ close(12)
+
+ ! eta model
+ call mpi_reduce(maxval(model_eta_new),max_eta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_eta_new),min_eta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'eta_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_eta_new
+ close(12)
+
+ ! rho model
+ call mpi_reduce(maxval(model_rho_new),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(model_rho_new),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+ fname = 'rho_new'
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
+ open(12,file=trim(m_file),form='unformatted')
+ write(12) model_rho_new
+ close(12)
+
+
+ if( myrank == 0 ) then
+ print*,'new models:'
+ print*,' vpv min/max: ',min_vpv,max_vpv
+ print*,' vph min/max: ',min_vph,max_vph
+ print*,' vsv min/max: ',min_vsv,max_vsv
+ print*,' vsh min/max: ',min_vsh,max_vsh
+ print*,' eta min/max: ',min_eta,max_eta
+ print*,' rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+
+end subroutine store_new_model
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine store_perturbations()
+
+! file output for new model
+
+ use model_update_iso
+ implicit none
+
+ ! vpv relative perturbations
+ ! logarithmic perturbation: log( v_new) - log( v_old) = log( v_new / v_old )
+ total_model = 0.0_CUSTOM_REAL
+ where( model_vpv /= 0.0 ) total_model = log( model_vpv_new / model_vpv)
+ ! or
+ ! linear approximation: (v_new - v_old) / v_old
+ !where( model_vpv /= 0.0 ) total_model = ( model_vpv_new - model_vpv) / model_vpv
+
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvpvvpv.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_vpv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vpv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! vph relative perturbations
+ total_model = 0.0_CUSTOM_REAL
+ where( model_vph /= 0.0 ) total_model = log( model_vph_new / model_vph)
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvphvph.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_vph,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vph,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! vsv relative perturbations
+ total_model = 0.0_CUSTOM_REAL
+ where( model_vsv /= 0.0 ) total_model = log( model_vsv_new / model_vsv)
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvsvvsv.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_vsv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vsv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! vsh relative perturbations
+ total_model = 0.0_CUSTOM_REAL
+ where( model_vsh /= 0.0 ) total_model = log( model_vsh_new / model_vsh)
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvshvsh.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_vsh,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_vsh,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! eta relative perturbations
+ total_model = 0.0_CUSTOM_REAL
+ where( model_eta /= 0.0 ) total_model = log( model_eta_new / model_eta)
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_detaeta.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_eta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_eta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ ! rho relative model perturbations
+ total_model = 0.0_CUSTOM_REAL
+ where( model_rho /= 0.0 ) total_model = log( model_rho_new / model_rho)
+ write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_drhorho.bin'
+ open(12,file=trim(m_file),form='unformatted',action='write')
+ write(12) total_model
+ close(12)
+ call mpi_reduce(maxval(total_model),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(minval(total_model),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
+
+ if( myrank == 0 ) then
+ print*,'relative update:'
+ print*,' dvpv/vpv min/max: ',min_vpv,max_vpv
+ print*,' dvph/vph min/max: ',min_vph,max_vph
+ print*,' dvsv/vsv min/max: ',min_vsv,max_vsv
+ print*,' dvsh/vsh min/max: ',min_vsh,max_vsh
+ print*,' deta/eta min/max: ',min_eta,max_eta
+ print*,' drho/rho min/max: ',min_rho,max_rho
+ print*
+ endif
+
+end subroutine store_perturbations
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/go_globe_pbs.bash
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/go_globe_pbs.bash (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/go_globe_pbs.bash 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,59 @@
+#!/bin/bash
+#PBS -S /bin/bash
+
+## job name and output file
+#PBS -N smooth_kernel
+#PBS -j oe
+#PBS -o OUTPUT_FILES/job_bulk_c.o
+
+###########################################################
+# USER PARAMETERS
+
+## 144 CPUs ( 18*8 ), walltime 10 hour
+#PBS -l nodes=18:ppn=8,walltime=10:00:00
+
+numnodes=144
+
+## horizontal: sigmah (in km)
+## vertical: sigmav (in km)
+## e.g. period 40 s: wavelength lambda = 40 s * 4 km/s ~ 160 km
+sigmah=160
+sigmav=40
+
+# kernel to smooth
+kernel=bulk_c_kernel
+
+###########################################################
+
+# (takes about: 1h 48 min for sigmah=160/sigmav=40 ...)
+# 25 min for sigmah=40/sigmav=40 ... )
+
+
+echo `date`
+echo "kernel: $kernel"
+cd $PBS_O_WORKDIR
+
+# obtain lsf job information
+cat $PBS_NODEFILE > OUTPUT_FILES/compute_nodes
+echo "$PBS_JOBID" > OUTPUT_FILES/jobid
+
+
+# the neigboring points within 3*sigma+element_size are used for the smoothing
+# the 3D Gaussian function is defined as
+# G(x,y,z) = 1/(sqrt(2*pi)*sigma)**3 * exp[-r**2/(2*sigma**2)]
+# which is parallel to the 1D Gaussian function
+
+# the horizontal smoothing radius would be sigma_h:
+# 160 km (~about wavelength 4km * 40 s )
+# the vertical smoothing radius would be sigma_v:
+# 40 km ( ~about quater wavelength 4km * 40 s / 4 )
+
+# usage: smooth_sem_globe sigma_h sigma_v kernel_name kernel_input_dir/ topo_dir/
+
+echo "kernel smoothing: $kernel"
+mpiexec -np $numnodes $PWD/src/smooth_sem_globe $sigmah $sigmav $kernel $PWD/OUTPUT_SUM $PWD/topo
+echo
+
+
+echo "done successfully"
+echo `date`
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/go_globe_pbs.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/readme_globe
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/readme_globe (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/readme_globe 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,62 @@
+-----------------------------------------------
+readme
+-----------------------------------------------
+
+
+smooth_sem_globe:
+
+ reads in GLOBAL mesh and kernels, in parallel on all processors, and
+ smooths by a three-dimensional gaussian
+
+1. create symbolic link 'topo' to link to topography files:
+ > ln -s ../topo
+ (make sure ../topo links to files proc_**_save_data_1.bin and proc_**_save_data_2.bin)
+
+2. copy file SPECFEM3D_GLOBE/constants.h
+ over to: src/constants_globe.h
+
+ copy file SPECFEM3D_GLOBE/precision.h
+ over to: src/precision_globe.h
+
+ copy file SPECFEM3D_GLOBE/OUTPUT_FILES/values_from_mesher.h
+ over to: src/values_from_mesher_globe.h
+ ( when generated for a kernel run ...)
+
+ copy file SPECFEM3D_GLOBE/config.h
+ over to: src/config_globe.h
+
+
+3. edit file 'smooth_sem_globe.f90':
+ set your
+ NPROC_XI
+ NPROC_ETA
+ NCHUNKS
+ element_size
+ accordingly
+
+
+then, compile with:
+ > cd src/
+ > make -f Makefile_globe
+
+
+4. create symbolic link directory 'OUTPUT_SUM/'
+ which holds the summed kernel files proc**_alpha_kernel.bin, etc.as outputted by sum_kernel/ executable
+
+
+ run with:
+ ---------
+ > ./run_globe.bash
+
+ (will run for about 2 hours...)
+
+ which creates smoothed kernel partitions "_smooth.bin" in the same directory
+ as the input ones.
+
+ check size of smoothing in file go_globe_pbs.bash:
+
+ - example: for periods 20s - 100s: minimum wavelength lambda ~ 4 km/s * 20 s = 80 km
+ then smoothing sizes become
+ horizontal = lambda
+ vertical = lambda/4
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/run_globe.bash
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/run_globe.bash (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/run_globe.bash 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,32 @@
+#!/bin/bash
+
+# compiles executable in src/ directory
+current_pwd=$PWD
+cd src/
+make -f Makefile_globe
+cd ../
+
+# cleans up outputs
+rm -f OUTPUT_FILES/*
+
+# runs job
+date
+
+# creates job scripts from template:
+template=go_globe_pbs.bash
+
+for tag in "bulk_c" "bulk_beta" "rho"
+do
+
+ echo "kernel: $tag"
+
+ # prepares job script from template
+ sed -e "s:bulk_c:$tag:" $template > go_globe_pbs.${tag}.bash
+
+ # smooths: bulk_c, bulk_beta and rho kernels
+ qsub go_globe_pbs.${tag}.bash
+ sleep 2
+
+done
+
+
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/run_globe.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/Makefile_globe
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/Makefile_globe (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/Makefile_globe 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,10 @@
+FFLAGS=-O3 -assume byterecl #-g -traceback
+
+all: smooth_sem_globe
+
+smooth_sem_globe: smooth_sem_globe.f90
+ mpif90 -o smooth_sem_globe $(FFLAGS) smooth_sem_globe.f90 smooth_sub.f90 exit_mpi.f90 gll_library_globe.f90
+
+
+clean:
+ rm -f smooth_sem_globe *.o
\ No newline at end of file
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/smooth_sem_globe.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/smooth_sem_globe.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/smooth_sem_globe.f90 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,595 @@
+! smooth_sem_globe
+!
+! this program can be used for smoothing a (summed) event kernel,
+! where it smooths files with a given input kernel name:
+!
+! Usage:
+! ./smooth_sem_globe sigma_h(km) sigma_v(km) kernel_file_name scratch_file_dir scratch_topo_dir
+! e.g.
+! ./smooth_sem_globe 160 10 bulk_c_kernel OUTPUT_SUM/ topo/
+!
+! where:
+! sigma_h - gaussian width for horizontal smoothing (in km)
+! sigma_v - gaussian width for vertical smoothing (in km)
+! kernel_file_name - takes file with this kernel name,
+! e.g. "bulk_c_kernel"
+! scratch_file_dir - directory containing kernel files,
+! e.g. proc***_reg1_bulk_c_kernel.bin
+! topo_dir - directory containing mesh topo files:
+! proc***_solver_data1.bin, proc***_solver_data2.bin
+! outputs:
+! puts the resulting, smoothed kernel files into the same directory as scratch_file_dir/
+! with a file ending "proc***_kernel_smooth.bin"
+
+program smooth_sem_globe
+
+! this is the embarassingly-parallel program that smooths any specfem function (primarily the kernels)
+! that has the dimension of (NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT)
+!
+! notice that it uses the constants_globe.h and precision_globe.h files
+! from the original SPECFEM3D_GLOBE package, and the
+! values_from_mesher_globe.h file from the output of the mesher (or create_header_file),
+! therefore, you need to compile it for your specific case
+!
+! NOTE: smoothing can be different in radial & horizontal directions; mesh is in spherical geometry.
+! algorithm uses vector components in radial/horizontal direction
+
+ implicit none
+ include 'mpif.h'
+ include 'constants_globe.h'
+ include 'precision_globe.h'
+ include 'values_from_mesher_globe.h'
+
+! ======================================================
+ ! USER PARAMETERS
+
+ ! copy from Par_file values
+ integer, parameter :: NPROC_XI = 12
+ integer, parameter :: NPROC_ETA = 12
+ integer, parameter :: NCHUNKS = 1
+
+ ! taken from values_from_mesher.h:
+ ! average size of a spectral element in km = ...
+ ! e.g. nproc 12x12, nex 192: element_size = 52.122262
+ real(kind=CUSTOM_REAL),parameter:: element_size = 52.12262
+
+! ======================================================
+
+ !takes region 1 kernels
+ integer, parameter :: NSPEC_MAX = NSPEC_CRUST_MANTLE_ADJOINT
+ integer, parameter :: NGLOB_MAX = NGLOB_CRUST_MANTLE
+
+ ! only include the neighboring 3 x 3 slices
+ integer, parameter :: NSLICES = 3
+ integer ,parameter :: NSLICES2 = NSLICES * NSLICES
+
+ character(len=256) :: s_sigma_h, s_sigma_v
+ character(len=256) :: kernel_file_name, scratch_topo_dir, scratch_file_dir
+ integer :: sizeprocs,ier,myrank,ichunk, ixi, ieta, iglob
+ integer :: islice(NSLICES2), islice0(NSLICES2), nums
+
+ real(kind=CUSTOM_REAL) :: sigma_h, sigma_h2, sigma_h3, sigma_v, sigma_v2, sigma_v3
+
+ real(kind=CUSTOM_REAL) :: x0, y0, z0, norm, norm_h, norm_v, element_size_m, max_old, max_new
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: factor, exp_val
+
+ character(len=256) :: ks_file, reg_name
+ character(len=256), dimension(NSLICES2) :: x_file, y_file, z_file, j_file, k_file, i_file
+ character(len=256), dimension(NSLICES2) :: solver1_file,solver2_file
+
+ logical :: global_code
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_MAX) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_MAX) :: kernel, kernel_smooth
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_MAX) :: tk, bk, jacobian, xl, yl, zl, xx, yy, zz
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_MAX) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ real(kind=CUSTOM_REAL), dimension(NGLOB_MAX) :: x, y, z
+ real(kind=CUSTOM_REAL), dimension(NSPEC_MAX) :: cx0, cy0, cz0, cx, cy, cz
+
+ integer :: i,ii,j,jj,k,kk,ispec,iproc,ispec2,nspec(NSLICES2),nglob(NSLICES2)
+
+ ! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: xigll, wxgll
+ double precision, dimension(NGLLY) :: yigll, wygll
+ double precision, dimension(NGLLZ) :: zigll, wzgll
+
+ ! array with all the weights in the cube
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+
+ real(kind=CUSTOM_REAL) :: dist_h,dist_v
+ real(kind=CUSTOM_REAL) :: r1,theta1
+
+ ! ============ program starts here =====================
+
+ ! initialize the MPI communicator and start the NPROCTOT MPI processes
+ call MPI_INIT(ier)
+ call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
+ call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
+
+ if (myrank == 0) print*,"smooth:"
+ call mpi_barrier(MPI_COMM_WORLD,ier)
+
+ ! arguments
+ call getarg(1,s_sigma_h)
+ call getarg(2,s_sigma_v)
+ call getarg(3,kernel_file_name)
+ call getarg(4,scratch_file_dir)
+ call getarg(5,scratch_topo_dir)
+
+ if ( trim(s_sigma_h) == '' .or. trim(s_sigma_v) == '' &
+ .or. trim(kernel_file_name) == '' &
+ .or. trim(scratch_file_dir) == '' &
+ .or. trim(scratch_topo_dir) == '') then
+ call exit_MPI(myrank,'Usage: smooth_sem_globe sigma_h(km) sigma_v(km) kernel_file_name scratch_file_dir scratch_topo_dir')
+ endif
+
+ ! read in parameter information
+ read(s_sigma_h,*) sigma_h
+ read(s_sigma_v,*) sigma_v
+
+ ! checks if basin code or global code: global code uses nchunks /= 0
+ if (NCHUNKS == 0) then
+ global_code = .false.
+ call exit_mpi(myrank,'Error nchunks')
+ else
+ global_code = .true.
+ reg_name='_reg1_'
+ endif
+ if (sizeprocs /= NPROC_XI*NPROC_ETA*NCHUNKS) call exit_mpi(myrank,'Error total number of slices')
+
+ ! user output
+ if (myrank == 0) then
+ print*,"defaults:"
+ print*," NPROC_XI , NPROC_ETA: ",NPROC_XI,NPROC_ETA
+ print*," NCHUNKS : ",NCHUNKS
+ print*," element size on surface(km): ",element_size
+ print*," smoothing sigma_h , sigma_v: ",sigma_h,sigma_v
+ endif
+ ! synchronizes
+ call mpi_barrier(MPI_COMM_WORLD,ier)
+
+ ! initializes lengths
+ element_size_m = element_size * 1000 ! e.g. 9 km on the surface, 36 km at CMB
+ if (global_code) element_size_m = element_size_m/R_EARTH
+
+ sigma_h = sigma_h * 1000.0 ! m
+ if (global_code) sigma_h = sigma_h / R_EARTH ! scale
+ sigma_v = sigma_v * 1000.0 ! m
+ if (global_code) sigma_v = sigma_v / R_EARTH ! scale
+
+ sigma_h2 = 2.0 * sigma_h ** 2 ! factor two for gaussian distribution with standard variance sigma
+ sigma_v2 = 2.0 * sigma_v ** 2
+
+ ! search radius
+ sigma_h3 = 3.0 * sigma_h + element_size_m
+ sigma_v3 = 3.0 * sigma_v + element_size_m
+
+ ! theoretic normal value
+ ! (see integral over -inf to +inf of exp[- x*x/(2*sigma) ] = sigma * sqrt(2*pi) )
+
+! note: smoothing is using a gaussian (ellipsoid for sigma_h /= sigma_v),
+! but in spherical coordinates, we use horizontal distance as epicentral distance
+! and vertical distance as radial distance?
+
+! not squared since epicentral distance is taken? values from bk seem to be closer to squared ones...
+ !norm_h = sqrt(2.0*PI) * sigma_h
+ norm_h = 2.0*PI*sigma_h**2
+ norm_v = sqrt(2.0*PI) * sigma_v
+ norm = norm_h * norm_v
+ !norm = (sqrt(2.0*PI) * sigma) ** 3 ! for sigma_h = sigma_v = sigma
+
+
+ ! GLL points weights
+ call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ wgll_cube(i,j,k) = wxgll(i)*wygll(j)*wzgll(k)
+ enddo
+ enddo
+ enddo
+
+ ! ---- figure out the neighboring 8 or 7 slices: (ichunk,ixi,ieta) index start at 0------
+ ichunk = myrank / (NPROC_XI * NPROC_ETA)
+ ieta = (myrank - ichunk * NPROC_XI * NPROC_ETA) / NPROC_XI
+ ixi = myrank - ichunk * NPROC_XI * NPROC_ETA - ieta * NPROC_XI
+
+ ! get the neighboring slices:
+ call get_all_eight_slices(ichunk,ixi,ieta,&
+ islice0(1),islice0(2),islice0(3),islice0(4),islice0(5),islice0(6),islice0(7),islice0(8),&
+ NPROC_XI,NPROC_ETA)
+
+ ! remove the repeated slices (only 8 for corner slices in global case)
+ islice(1) = myrank; j = 1
+ do i = 1, 8
+ if (.not. any(islice(1:i) == islice0(i)) .and. islice0(i) < sizeprocs) then
+ j = j + 1
+ islice(j) = islice0(i)
+ endif
+ enddo
+ nums = j
+
+ if( myrank == 0 ) then
+ print *,'slices:',nums
+ print *,' ',islice(:)
+ print *
+ endif
+
+ ! read in the topology files of the current and neighboring slices
+ do i = 1, nums
+ write(k_file(i),'(a,i6.6,a)') &
+ trim(scratch_file_dir)//'/proc',islice(i),trim(reg_name)//trim(kernel_file_name)//'.bin'
+
+ write(solver1_file(i),'(a,i6.6,a)') &
+ trim(scratch_topo_dir)//'/proc',islice(i),trim(reg_name)//'solver_data_1.bin'
+ write(solver2_file(i),'(a,i6.6,a)') &
+ trim(scratch_topo_dir)//'/proc',islice(i),trim(reg_name)//'solver_data_2.bin'
+
+ nspec(i) = NSPEC_MAX
+ nglob(i) = NGLOB_MAX
+ enddo
+
+ ! point locations
+ open(11,file=solver2_file(1),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening solver2 file')
+
+ read(11) x(1:nglob(1))
+ read(11) y(1:nglob(1))
+ read(11) z(1:nglob(1))
+ read(11) ibool(:,:,:,1:nspec(1))
+ close(11)
+
+ ! jacobian
+ open(11,file=solver1_file(1),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening solver1 file')
+
+ read(11) xix
+ read(11) xiy
+ read(11) xiz
+ read(11) etax
+ read(11) etay
+ read(11) etaz
+ read(11) gammax
+ read(11) gammay
+ read(11) gammaz
+ close(11)
+
+ ! get the location of the center of the elements
+ do ispec = 1, nspec(1)
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+ xl(i,j,k,ispec) = x(iglob)
+ yl(i,j,k,ispec) = y(iglob)
+ zl(i,j,k,ispec) = z(iglob)
+
+ ! build jacobian
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ ! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ jacobian(i,j,k,ispec) = jacobianl
+
+
+
+ enddo
+ enddo
+ enddo
+ cx0(ispec) = (xl(1,1,1,ispec) + xl(NGLLX,NGLLY,NGLLZ,ispec))/2.0
+ cy0(ispec) = (yl(1,1,1,ispec) + yl(NGLLX,NGLLY,NGLLZ,ispec))/2.0
+ cz0(ispec) = (zl(1,1,1,ispec) + zl(NGLLX,NGLLY,NGLLZ,ispec))/2.0
+ enddo
+
+ if (myrank == 0) write(*,*) 'start looping over elements and points for smoothing ...'
+
+ ! smoothed kernel file name
+ write(ks_file,'(a,i6.6,a)') trim(scratch_file_dir)//'/proc',myrank, &
+ trim(reg_name)//trim(kernel_file_name)//'_smooth.bin'
+
+ tk = 0.0
+ bk = 0.0
+ kernel_smooth=0.0
+
+ ! loop over all the slices
+ do iproc = 1, nums
+
+ ! read in the topology, kernel files, calculate center of elements
+ ! point locations
+ ! given in cartesian coordinates
+ open(11,file=solver2_file(iproc),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening slices: solver2 file')
+
+ read(11) x(1:nglob(iproc))
+ read(11) y(1:nglob(iproc))
+ read(11) z(1:nglob(iproc))
+ read(11) ibool(:,:,:,1:nspec(iproc))
+ close(11)
+
+ open(11,file=solver1_file(iproc),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening slices: solver1 file')
+
+ read(11) xix
+ read(11) xiy
+ read(11) xiz
+ read(11) etax
+ read(11) etay
+ read(11) etaz
+ read(11) gammax
+ read(11) gammay
+ read(11) gammaz
+ close(11)
+
+ ! get the location of the center of the elements
+ do ispec = 1, nspec(iproc)
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ ! build jacobian
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ ! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+ jacobian(i,j,k,ispec) = jacobianl
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! kernel file
+ open(11,file=k_file(iproc),status='old',form='unformatted',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening kernel file')
+
+ read(11) kernel(:,:,:,1:nspec(iproc))
+ close(11)
+
+
+ ! get the global maximum value of the original kernel file
+ if (iproc == 1) then
+ call mpi_reduce(maxval(abs(kernel(:,:,:,1:nspec(iproc)))), max_old, 1, &
+ CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+ endif
+
+ ! calculate element center location
+ do ispec2 = 1, nspec(iproc)
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec2)
+ xx(i,j,k,ispec2) = x(iglob)
+ yy(i,j,k,ispec2) = y(iglob)
+ zz(i,j,k,ispec2) = z(iglob)
+ enddo
+ enddo
+ enddo
+ cx(ispec2) = (xx(1,1,1,ispec2) + xx(NGLLX,NGLLZ,NGLLY,ispec2))/2.0
+ cy(ispec2) = (yy(1,1,1,ispec2) + yy(NGLLX,NGLLZ,NGLLY,ispec2))/2.0
+ cz(ispec2) = (zz(1,1,1,ispec2) + zz(NGLLX,NGLLZ,NGLLY,ispec2))/2.0
+ enddo
+
+ ! loop over elements to be smoothed in the current slice
+ do ispec = 1, nspec(1)
+
+ ! --- only double loop over the elements in the search radius ---
+ do ispec2 = 1, nspec(iproc)
+
+ ! calculates horizontal and vertical distance between two element centers
+
+ ! vector approximation
+ call get_distance_vec(dist_h,dist_v,cx0(ispec),cy0(ispec),cz0(ispec),&
+ cx(ispec2),cy(ispec2),cz(ispec2))
+
+ ! note: distances and sigmah, sigmav are normalized by R_EARTH
+
+ ! checks distance between centers of elements
+ if ( dist_h > sigma_h3 .or. abs(dist_v) > sigma_v3 ) cycle
+
+ ! integration factors
+ factor(:,:,:) = jacobian(:,:,:,ispec2) * wgll_cube(:,:,:)
+
+
+ ! loop over GLL points of the elements in current slice (ispec)
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ ! current point (i,j,k,ispec) location, cartesian coordinates
+ x0 = xl(i,j,k,ispec)
+ y0 = yl(i,j,k,ispec)
+ z0 = zl(i,j,k,ispec)
+
+ ! calculate weights based on gaussian smoothing
+ call smoothing_weights_vec(x0,y0,z0,ispec2,sigma_h2,sigma_v2,exp_val,&
+ xx(:,:,:,ispec2),yy(:,:,:,ispec2),zz(:,:,:,ispec2))
+
+ ! adds GLL integration weights
+ exp_val(:,:,:) = exp_val(:,:,:) * factor(:,:,:)
+
+ ! adds contribution of element ispec2 to smoothed kernel values
+ tk(i,j,k,ispec) = tk(i,j,k,ispec) + sum(exp_val(:,:,:) * kernel(:,:,:,ispec2))
+
+ ! normalization, integrated values of gaussian smoothing function
+ bk(i,j,k,ispec) = bk(i,j,k,ispec) + sum(exp_val(:,:,:))
+
+ enddo
+ enddo
+ enddo ! (i,j,k)
+ enddo ! (ispec2)
+ enddo ! (ispec)
+ enddo ! islice
+
+ if (myrank == 0) write(*,*) 'Done with integration ...'
+
+ ! compute the smoothed kernel values
+ do ispec = 1, nspec(1)
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ ! checks the normalization criterion
+ ! e.g. sigma_h 160km, sigma_v 40km:
+ ! norm (not squared sigma_h ) ~ 0.001
+ ! norm ( squared sigma_h) ~ 6.23 * e-5
+ if (abs(bk(i,j,k,ispec) - norm) > 1.e-4 ) then
+ print *, 'Problem norm here --- ', myrank, ispec, i, j, k, bk(i,j,k,ispec), norm
+ !call exit_mpi(myrank, 'Error computing Gaussian function on the grid')
+ endif
+
+ ! normalizes smoothed kernel values by integral value of gaussian weighting
+ kernel_smooth(i,j,k,ispec) = tk(i,j,k,ispec) / bk(i,j,k,ispec)
+
+ enddo
+ enddo
+ enddo
+ enddo
+ if (myrank == 0) write(*,*) ' norm: ',norm
+
+ ! file output
+ open(11,file=trim(ks_file),status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening smoothed kernel file')
+
+ ! Note: output the following instead of kernel_smooth(:,:,:,1:nspec(1)) to create files of the same sizes
+ write(11) kernel_smooth(:,:,:,:)
+ close(11)
+
+ if (myrank == 0) print *,' written:',trim(ks_file)
+
+
+ ! the maximum value for the smoothed kernel
+ call mpi_reduce(maxval(abs(kernel_smooth(:,:,:,1:nspec(1)))), max_new, 1, &
+ CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+
+ if (myrank == 0) then
+ print *
+ print *, 'Maximum data value before smoothing = ', max_old
+ print *, 'Maximum data value after smoothing = ', max_new
+ print *
+ endif
+
+! stop all the MPI processes, and exit
+ call MPI_FINALIZE(ier)
+
+end program smooth_sem_globe
+
+!
+! -----------------------------------------------------------------------------
+!
+ subroutine smoothing_weights_vec(x0,y0,z0,ispec2,sigma_h2,sigma_v2,exp_val,&
+ xx_elem,yy_elem,zz_elem)
+
+ implicit none
+ include "constants_globe.h"
+
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ),intent(out) :: exp_val
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ),intent(in) :: xx_elem, yy_elem, zz_elem
+ real(kind=CUSTOM_REAL),intent(in) :: x0,y0,z0,sigma_h2,sigma_v2
+ integer,intent(in) :: ispec2
+
+ ! local parameters
+ integer :: ii,jj,kk
+ real(kind=CUSTOM_REAL) :: dist_h,dist_v
+ !real(kind=CUSTOM_REAL) :: r0,r1,theta1
+
+ ! >>>>>
+ ! uniform sigma
+ !exp_val(:,:,:) = exp( -((xx(:,:,:,ispec2)-x0)**2+(yy(:,:,:,ispec2)-y0)**2 &
+ ! +(zz(:,:,:,ispec2)-z0)**2 )/(2*sigma2) )*factor(:,:,:)
+
+ ! from basin code smoothing:
+ ! gaussian function
+ !exp_val(:,:,:) = exp( -(xx(:,:,:,ispec2)-x0)**2/(sigma_h2) &
+ ! -(yy(:,:,:,ispec2)-y0)**2/(sigma_h2) &
+ ! -(zz(:,:,:,ispec2)-z0)**2/(sigma_v2) ) * factor(:,:,:)
+ ! >>>>>
+
+ do kk = 1, NGLLZ
+ do jj = 1, NGLLY
+ do ii = 1, NGLLX
+ ! point in second slice
+
+ ! vector approximation:
+ call get_distance_vec(dist_h,dist_v,x0,y0,z0, &
+ xx_elem(ii,jj,kk),yy_elem(ii,jj,kk),zz_elem(ii,jj,kk))
+
+ ! gaussian function
+ exp_val(ii,jj,kk) = exp( - dist_h*dist_h/sigma_h2 &
+ - dist_v*dist_v/sigma_v2 ) ! * factor(ii,jj,kk)
+
+ enddo
+ enddo
+ enddo
+
+ end subroutine smoothing_weights_vec
+
+
+!
+! -----------------------------------------------------------------------------
+!
+
+ subroutine get_distance_vec(dist_h,dist_v,x0,y0,z0,x1,y1,z1)
+
+! returns vector lengths as distances in radial and horizontal direction
+
+ implicit none
+ include "constants_globe.h"
+
+ real(kind=CUSTOM_REAL),intent(out) :: dist_h,dist_v
+ real(kind=CUSTOM_REAL),intent(in) :: x0,y0,z0,x1,y1,z1
+
+ ! local parameters
+ real(kind=CUSTOM_REAL) :: r0,r1,alpha
+ real(kind=CUSTOM_REAL) :: vx,vy,vz
+
+ ! vertical distance
+ r0 = sqrt( x0*x0 + y0*y0 + z0*z0 ) ! length of first position vector
+ r1 = sqrt( x1*x1 + y1*y1 + z1*z1 )
+ dist_v = r1 - r0
+ ! only for flat earth with z in depth: dist_v = sqrt( (cz(ispec2)-cz0(ispec))** 2)
+
+ ! horizontal distance
+ ! length of vector from point 0 to point 1
+ ! assuming small earth curvature (since only for neighboring elements)
+
+ ! scales r0 to have same length as r1
+ alpha = r1 / r0
+ vx = alpha * x0
+ vy = alpha * y0
+ vz = alpha * z0
+
+ ! vector in horizontal between new r0 and r1
+ vx = x1 - vx
+ vy = y1 - vy
+ vz = z1 - vz
+
+ ! distance is vector length
+ dist_h = sqrt( vx*vx + vy*vy + vz*vz )
+
+ end subroutine get_distance_vec
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/go_globe_pbs.bash
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/go_globe_pbs.bash (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/go_globe_pbs.bash 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,30 @@
+#!/bin/bash
+#PBS -S /bin/bash
+
+## job name and output file
+#PBS -N sum_kernel
+#PBS -j oe
+#PBS -o OUTPUT_FILES/job.o
+
+###########################################################
+# USER PARAMETERS
+
+## e.g. 144 CPUs ( 18*8 ), walltime 2 hour
+#PBS -l nodes=18:ppn=8,walltime=2:00:00
+
+numnodes=144
+
+###########################################################
+
+#(takes about 10 minutes...)
+
+echo `date`
+
+cd $PBS_O_WORKDIR
+
+# obtain lsf job information
+cat $PBS_NODEFILE > OUTPUT_FILES/compute_nodes
+echo "$PBS_JOBID" > OUTPUT_FILES/jobid
+
+mpiexec -np $numnodes $PWD/src/sum_kernels_globe
+
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/go_globe_pbs.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/readme_globe
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/readme_globe (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/readme_globe 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,46 @@
+------------------------------------
+readme
+------------------------------------
+
+sum_kernels_globe:
+
+ reads in several event kernels, in parallel on all processors,
+ and sums them all up
+
+
+1. copy file SPECFEM3D_GLOBE/constants.h
+ over to: src/constants_globe.h
+
+ copy file SPECFEM3D_GLOBE/precision.h
+ over to: src/precision_globe.h
+
+ copy file SPECFEM3D_GLOBE/OUTPUT_FILES/values_from_mesher.h
+ over to: src/values_from_mesher_globe.h
+ ( when generated for a kernel run ...)
+
+ copy file SPECFEM3D_GLOBE/config.h
+ over to: src/config_globe.h
+
+2. compile executable 'sum_kernels_globe':
+ > cd src/
+ > make -f Makefile_globe
+
+3. create symbolic link 'INPUT_KERNELS' to kernel directories:
+ > ln -s /my_OUTPUTS_KERNEL INPUT_KERNEL
+ or
+ > cd INPUT_KERNEL/
+ > ln -s /my_OUTPUTS_KERNEL/*.* ./
+
+4. copy each single kernel directory name to file 'kernels_run.globe':
+ > ls -1 INPUT_KERNELS/ > kernels_run.globe
+
+
+
+then, run job:
+
+ on sesame: > run_globe.bash
+
+ (will run for about 10 minutes)
+
+
+ find summed event kernels in directory OUTPUT_SUM/
\ No newline at end of file
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/run_globe.bash
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/run_globe.bash (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/run_globe.bash 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,33 @@
+#!/bin/bash
+
+###########################################################
+# USER PARAMETERS
+
+# set this to your directory which contains all event kernel directories
+outputs_kernel_directory=my_kernels
+
+###########################################################
+
+# clean up
+mkdir -p OUTPUT_FILES
+rm -f OUTPUT_FILES/*
+rm -f OUTPUT_SUM/*
+
+# update kernel links
+cd INPUT_KERNELS/
+rm -rf ./*.*
+ln -s ${outputs_kernel_directory}/*.* ./
+cd ..
+
+# update input kernel list
+ls -1 INPUT_KERNELS/ > kernels_run.globe
+
+# compiles
+cd src/
+make -f Makefile_globe
+cd ..
+
+# runs job
+date
+qsub go_globe_pbs.bash
+
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/run_globe.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/Makefile_globe
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/Makefile_globe (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/Makefile_globe 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,17 @@
+FC=mpif90
+FFLAGS=-O3 -assume byterecl
+
+OBJS = sum_kernels_globe.o exit_mpi.o
+
+all: sum_kernels_globe
+
+sum_kernels_globe: $(OBJS)
+ $(FC) $(FFLAGS) -o $@ $(OBJS)
+
+
+.SUFFIXES: .o .f90
+.f90.o:
+ $(FC) $(FFLAGS) -c $<
+
+clean:
+ rm -rf *.o *~ sum_kernels_globe
\ No newline at end of file
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/sum_kernels_globe.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/sum_kernels_globe.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/sum_kernel/src/sum_kernels_globe.f90 2010-10-12 22:46:52 UTC (rev 17262)
@@ -0,0 +1,192 @@
+! sum_kernels_globe
+!
+! this program can be used for event kernel summation,
+! where it sums up transverse isotropic kernel files:
+!
+! - proc***_reg1_bulk_c_kernel.bin
+! - proc***_reg1_bulk_betav_kernel.bin
+! - proc***_reg1_bulk_betah_kernel.bin
+! - proc***_reg1_eta_kernel.bin
+!
+! input file: kernels_run.globe
+! lists all event kernel directories which should be summed together
+!
+! input directory: INPUT_KERNELS/
+! contains links to all event kernel directories (listed in "kernel_list.globe")
+!
+! output directory: OUTPUT_SUM/
+! the resulting kernel files will be stored in this directory
+
+program sum_kernels_globe
+
+ implicit none
+ include 'mpif.h'
+ include 'constants_globe.h'
+ include 'precision_globe.h'
+ include 'values_from_mesher_globe.h'
+
+! ======================================================
+ ! USER PARAMETERS
+
+ ! by default, this algorithm uses transverse isotropic (bulk,bulk_betav,bulk_betah,eta) kernels to sum up
+ ! if you prefer using isotropic kernels, set flags below accordingly
+
+ ! if you prefer using isotropic kernels (bulk,bulk_beta,rho) kernels, set this flag to true
+ logical, parameter :: USE_ISO_KERNELS = .false.
+
+ ! if you prefer isotropic (alpha,beta,rho) kernels, set this flag to true
+ logical, parameter :: USE_ALPHA_BETA_RHO = .false.
+
+
+! ======================================================
+
+ ! crust mantle region 1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
+ kernel_crust_mantle,total_kernel
+
+ character(len=150) :: kernel_file_list, kernel_list(1000), sline, k_file, kernel_name
+ integer :: iker, nker, myrank, sizeprocs, ier
+ integer :: i, j, k,ispec, iglob, ishell, n, it, j1, ib, npts_sem, ios
+
+ ! ============ program starts here =====================
+ ! initialize the MPI communicator and start the NPROCTOT MPI processes
+ call MPI_INIT(ier)
+ call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
+ call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
+
+ if(myrank==0) write(*,*) 'reading kernel list: '
+
+ ! default values
+ kernel_file_list='kernels_run.globe'
+
+ ! reads in event list
+ nker=0
+ open(unit = 20, file = trim(kernel_file_list), status = 'old',iostat = ios)
+ if (ios /= 0) then
+ print *,'Error opening ',trim(kernel_file_list),myrank
+ stop 1
+ endif
+ do while (1 == 1)
+ read(20,'(a)',iostat=ios) sline
+ if (ios /= 0) exit
+ nker = nker+1
+ kernel_list(nker) = sline
+ enddo
+ close(20)
+ if( myrank == 0 ) then
+ write(*,*) ' ',nker,' events'
+ write(*,*)
+ endif
+
+ ! sums up kernels
+ if( USE_ISO_KERNELS ) then
+
+ ! isotropic kernels
+ kernel_name = 'reg1_bulk_c_kernel'
+ call sum_kernel(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_bulk_beta_kernel'
+ call sum_kernel(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_rho_kernel'
+ call sum_kernel(kernel_name,kernel_list,nker,myrank)
+
+ else if( USE_ALPHA_BETA_RHO ) then
+
+ ! isotropic kernels
+
+ kernel_name = 'reg1_alpha_kernel'
+ call sum_kernel(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_beta_kernel'
+ call sum_kernel(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_rho_kernel'
+ call sum_kernel(kernel_name,kernel_list,nker,myrank)
+
+ else
+
+ ! transverse isotropic kernels
+
+ kernel_name = 'reg1_bulk_c_kernel'
+ call sum_kernel(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_bulk_betav_kernel'
+ call sum_kernel(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_bulk_betah_kernel'
+ call sum_kernel(kernel_name,kernel_list,nker,myrank)
+
+ kernel_name = 'reg1_eta_kernel'
+ call sum_kernel(kernel_name,kernel_list,nker,myrank)
+
+ endif
+
+ if(myrank==0) write(*,*) 'done writing all kernels'
+
+ ! stop all the MPI processes, and exit
+ call MPI_FINALIZE(ier)
+
+end program sum_kernels_globe
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+subroutine sum_kernel(kernel_name,kernel_list,nker,myrank)
+
+ implicit none
+
+ include 'constants_globe.h'
+ include 'values_from_mesher_globe.h'
+
+ character(len=150) :: kernel_name,kernel_list(1000)
+ integer :: nker,myrank
+
+ ! local parameters
+ character(len=150) :: k_file
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
+ kernel_crust_mantle,total_kernel
+ integer :: iker,ios
+
+ total_kernel = 0._CUSTOM_REAL
+ do iker = 1, nker
+ if(myrank==0) then
+ write(*,*) 'reading in event kernel for: ',trim(kernel_name)
+ write(*,*) ' ',iker, ' out of ', nker
+ endif
+
+ write(k_file,'(a,i6.6,a)') 'INPUT_KERNELS/'//trim(kernel_list(iker)) &
+ //'/proc',myrank,'_'//trim(kernel_name)//'.bin'
+
+ open(12,file=trim(k_file),status='old',form='unformatted',action='read',iostat=ios)
+ if( ios /= 0 ) then
+ write(*,*) ' kernel not found:',trim(k_file)
+ cycle
+ endif
+
+ read(12) kernel_crust_mantle
+ close(12)
+
+ total_kernel(:,:,:,1:NSPEC_CRUST_MANTLE_ADJOINT) = &
+ total_kernel(:,:,:,1:NSPEC_CRUST_MANTLE_ADJOINT) &
+ + kernel_crust_mantle(:,:,:,1:NSPEC_CRUST_MANTLE_ADJOINT)
+ enddo
+
+ ! stores summed kernels
+ if(myrank==0) write(*,*) 'writing out summed kernel for: ',trim(kernel_name)
+
+ write(k_file,'(a,i6.6,a)') 'OUTPUT_SUM/proc',myrank,'_'//trim(kernel_name)//'.bin'
+
+ open(12,file=trim(k_file),form='unformatted',status='unknown',action='write',iostat=ios)
+ if( ios /= 0 ) then
+ write(*,*) 'ERROR kernel not written:',trim(k_file)
+ stop 'error kernel write'
+ endif
+ write(12) total_kernel
+ close(12)
+
+ if(myrank==0) write(*,*)
+
+end subroutine sum_kernel
+
More information about the CIG-COMMITS
mailing list