[cig-commits] r20901 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Wed Oct 24 06:48:12 PDT 2012
Author: surendra
Date: 2012-10-24 06:48:12 -0700 (Wed, 24 Oct 2012)
New Revision: 20901
Added:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90
Log:
Added common fault solver
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90 (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90 2012-10-24 13:48:12 UTC (rev 20901)
@@ -0,0 +1,610 @@
+! Base module for kinematic and dynamic fault solvers
+!
+! Written by: Percy Galvez, Surendra Somala and Jean-Paul Ampuero
+
+module fault_solver_common
+
+ implicit none
+ private
+
+ ! outputs on selected fault nodes at every time step:
+ ! slip, slip velocity, fault stresses
+ type dataT_type
+ integer :: npoin=0,npoin_local=0
+ integer, dimension(:), pointer :: iglob=>null() ! on-fault global index of output nodes
+ integer, dimension(:,:), pointer :: iglob_all=>null()
+ integer, dimension(:), pointer :: islice=>null(),glob_indx=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: dist=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: dist_all=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: d1=>null(),v1=>null(),t1=>null(), &
+ d2=>null(),v2=>null(),t2=>null(), &
+ t3=>null(),theta=>null()
+ character(len=70), dimension(:), pointer :: name=>null()
+ end type dataT_type
+
+
+ ! outputs(dyn) /inputs (kind) at selected times for all fault nodes:
+ ! strength, state, slip, slip velocity, fault stresses, rupture time, process zone time
+ ! rupture time = first time when slip velocity = threshold V_RUPT (defined below)
+ ! process zone time = first time when slip = Dc
+ type dataXZ_type
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: stg=>null(), sta=>null(), d1=>null(), d2=>null(), v1=>null(), v2=>null(), &
+ t1=>null(), t2=>null(), t3=>null(), tRUP=>null(), tPZ=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoord=>null(), ycoord=>null(), zcoord=>null()
+ integer :: npoin=0
+ end type dataXZ_type
+
+
+ type fault_type
+ private
+ integer :: nspec=0, nglob=0
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T0=>null(),T=>null(),V=>null(),D=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: coord=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:,:), pointer :: R=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: MU=>null(),B=>null(),invM1=>null(),invM2=>null(),Z=>null()
+ real(kind=CUSTOM_REAL) :: dt
+ integer, dimension(:), pointer :: ibulk1=>null(), ibulk2=>null()
+ type(dataT_type) :: dataT
+ type(dataXZ_type) :: dataXZ,dataXZ_all
+ integer, dimension(:), pointer :: poin_offset=>null(),npoin_perproc=>null()
+ end type fault_type
+
+
+!Number of time steps defined by the user : NTOUT
+ integer, save :: NTOUT,NSNAP
+
+
+ logical, save :: PARALLEL_FAULT = .false.
+
+ public :: initialize_fault,PARALLEL_FAULT
+
+contains
+
+!---------------------------------------------------------------------
+
+subroutine initialize_fault (bc,IIN_BIN,Minv,dt)
+
+ use specfem_par
+
+ type(bc_dynflt_type), intent(inout) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
+ integer, intent(in) :: IIN_BIN
+ real(kind=CUSTOM_REAL), intent(in) :: dt
+
+ real(kind=CUSTOM_REAL) :: tmp_vec(3,NGLOB_AB)
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: jacobian2Dw
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: normal
+ integer, dimension(:,:), allocatable :: ibool1
+ real(kind=CUSTOM_REAL) :: norm
+ integer :: n1,n2,n3
+ integer :: ij,k,e
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: nx,ny,nz
+
+ read(IIN_BIN) bc%nspec,bc%nglob
+ if (.NOT.PARALLEL_FAULT .and. bc%nspec==0) return
+ if (bc%nspec>0) then
+
+ allocate(bc%ibulk1(bc%nglob))
+ allocate(bc%ibulk2(bc%nglob))
+ allocate(bc%R(3,3,bc%nglob))
+ allocate(bc%coord(3,(bc%nglob)))
+ allocate(bc%invM1(bc%nglob))
+ allocate(bc%invM2(bc%nglob))
+ allocate(bc%B(bc%nglob))
+ allocate(bc%Z(bc%nglob))
+
+ allocate(ibool1(NGLLSQUARE,bc%nspec))
+ allocate(normal(NDIM,NGLLSQUARE,bc%nspec))
+ allocate(jacobian2Dw(NGLLSQUARE,bc%nspec))
+
+ read(IIN_BIN) ibool1
+ read(IIN_BIN) jacobian2Dw
+ read(IIN_BIN) normal
+ read(IIN_BIN) bc%ibulk1
+ read(IIN_BIN) bc%ibulk2
+ read(IIN_BIN) bc%coord(1,:)
+ read(IIN_BIN) bc%coord(2,:)
+ read(IIN_BIN) bc%coord(3,:)
+
+ bc%B = 0e0_CUSTOM_REAL
+ allocate(nx(bc%nglob),ny(bc%nglob),nz(bc%nglob))
+ nx = 0e0_CUSTOM_REAL
+ ny = 0e0_CUSTOM_REAL
+ nz = 0e0_CUSTOM_REAL
+ do e=1,bc%nspec
+ do ij = 1,NGLLSQUARE
+ k = ibool1(ij,e)
+ nx(k) = nx(k) + normal(1,ij,e)
+ ny(k) = ny(k) + normal(2,ij,e)
+ nz(k) = nz(k) + normal(3,ij,e)
+ bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
+ enddo
+ enddo
+ endif
+
+ if (PARALLEL_FAULT) then
+
+ tmp_vec = 0._CUSTOM_REAL
+ if (bc%nspec>0) tmp_vec(1,bc%ibulk1) = bc%B
+ ! assembles with other MPI processes
+ call assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,tmp_vec, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ my_neighbours_ext_mesh)
+ if (bc%nspec>0) bc%B = tmp_vec(1,bc%ibulk1)
+
+ tmp_vec = 0._CUSTOM_REAL
+ if (bc%nspec>0) then
+ tmp_vec(1,bc%ibulk1) = nx
+ tmp_vec(2,bc%ibulk1) = ny
+ tmp_vec(3,bc%ibulk1) = nz
+ endif
+ ! assembles with other MPI processes
+ call assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,tmp_vec, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ my_neighbours_ext_mesh)
+ if (bc%nspec>0) then
+ nx = tmp_vec(1,bc%ibulk1)
+ ny = tmp_vec(2,bc%ibulk1)
+ nz = tmp_vec(3,bc%ibulk1)
+ endif
+
+ endif
+
+ if (bc%nspec>0) then
+ do k=1,bc%nglob
+ norm = sqrt( nx(k)*nx(k) + ny(k)*ny(k) + nz(k)*nz(k) )
+ nx(k) = nx(k) / norm
+ ny(k) = ny(k) / norm
+ nz(k) = nz(k) / norm
+ enddo
+
+ call compute_R(bc%R,bc%nglob,nx,ny,nz)
+
+ ! Needed in dA_Free = -K2*d2/M2 + K1*d1/M1
+ bc%invM1 = Minv(bc%ibulk1)
+ bc%invM2 = Minv(bc%ibulk2)
+
+ ! Fault impedance, Z in : Trac=T_Stick-Z*dV
+ ! Z = 1/( B1/M1 + B2/M2 ) / (0.5*dt)
+ ! T_stick = Z*Vfree traction as if the fault was stuck (no displ discontinuity)
+ ! NOTE: same Bi on both sides, see note above
+ bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
+ ! WARNING: In non-split nodes at fault edges M is assembled across the fault.
+ ! hence invM1+invM2=2/(M1+M2) instead of 1/M1+1/M2
+ ! In a symmetric mesh (M1=M2) Z will be twice its intended value
+
+end subroutine initialize
+
+!---------------------------------------------------------------------
+subroutine compute_R(R,nglob,nx,ny,nz)
+
+ integer :: nglob
+ real(kind=CUSTOM_REAL), intent(out) :: R(3,3,nglob)
+ real(kind=CUSTOM_REAL), dimension(nglob), intent(in) :: nx,ny,nz
+
+ real(kind=CUSTOM_REAL), dimension(nglob) :: sx,sy,sz,dx,dy,dz,norm
+
+ ! Percy , defining fault directions (in concordance with SCEC conventions) .
+ ! fault coordinates (s,d,n) = (1,2,3)
+ ! s = strike , d = dip , n = n.
+ ! 1 = strike , 2 = dip , 3 = n.
+ norm = sqrt(nx*nx+ny*ny)
+ sx = ny/norm
+ sy = -nx/norm
+ sz = 0.e0_CUSTOM_REAL
+
+ norm = sqrt(sy*sy*nz*nz+sx*sx*nz*nz+(sy*nx-ny*sx)*(nx*sy-ny*sx))
+ dx = -sy*nz/norm
+ dy = sx*nz/norm
+ dz = (sy*nx-ny*sx)/norm
+ !Percy, dz is always dipwards = -1/norm , because (nx*sy-ny*sx)= - 1
+
+ R(1,1,:)=sx
+ R(1,2,:)=sy
+ R(1,3,:)=sz
+ R(2,1,:)=dx
+ R(2,2,:)=dy
+ R(2,3,:)=dz
+ R(3,1,:)=nx
+ R(3,2,:)=ny
+ R(3,3,:)=nz
+
+end subroutine compute_R
+
+!===============================================================
+function get_jump (bc,v) result(dv)
+
+ type(bc_dynflt_type), intent(in) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: v(:,:)
+ real(kind=CUSTOM_REAL) :: dv(3,bc%nglob)
+
+ ! diference between side 2 and side 1 of fault nodes. dv
+ dv(1,:) = v(1,bc%ibulk2)-v(1,bc%ibulk1)
+ dv(2,:) = v(2,bc%ibulk2)-v(2,bc%ibulk1)
+ dv(3,:) = v(3,bc%ibulk2)-v(3,bc%ibulk1)
+
+end function get_jump
+
+!---------------------------------------------------------------------
+function get_weighted_jump (bc,f) result(da)
+
+ type(bc_dynflt_type), intent(in) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: f(:,:)
+
+ real(kind=CUSTOM_REAL) :: da(3,bc%nglob)
+
+ ! diference between side 2 and side 1 of fault nodes. M-1 * F
+ da(1,:) = bc%invM2*f(1,bc%ibulk2)-bc%invM1*f(1,bc%ibulk1)
+ da(2,:) = bc%invM2*f(2,bc%ibulk2)-bc%invM1*f(2,bc%ibulk1)
+ da(3,:) = bc%invM2*f(3,bc%ibulk2)-bc%invM1*f(3,bc%ibulk1)
+
+ ! NOTE: In non-split nodes at fault edges M and f are assembled across the fault.
+ ! Hence, f1=f2, invM1=invM2=1/(M1+M2) instead of invMi=1/Mi, and da=0.
+
+end function get_weighted_jump
+
+!----------------------------------------------------------------------
+function rotate(bc,v,fb) result(vr)
+
+ type(bc_dynflt_type), intent(in) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: v(3,bc%nglob)
+ integer, intent(in) :: fb
+ real(kind=CUSTOM_REAL) :: vr(3,bc%nglob)
+
+ ! Percy, tangential direction Vt, equation 7 of Pablo's notes in agreement with SPECFEM3D
+
+ ! forward rotation
+ if (fb==1) then
+ vr(1,:) = v(1,:)*bc%R(1,1,:)+v(2,:)*bc%R(1,2,:)+v(3,:)*bc%R(1,3,:) ! vs
+ vr(2,:) = v(1,:)*bc%R(2,1,:)+v(2,:)*bc%R(2,2,:)+v(3,:)*bc%R(2,3,:) ! vd
+ vr(3,:) = v(1,:)*bc%R(3,1,:)+v(2,:)*bc%R(3,2,:)+v(3,:)*bc%R(3,3,:) ! vn
+
+ ! backward rotation
+ else
+ vr(1,:) = v(1,:)*bc%R(1,1,:)+v(2,:)*bc%R(2,1,:)+v(3,:)*bc%R(3,1,:) !vx
+ vr(2,:) = v(1,:)*bc%R(1,2,:)+v(2,:)*bc%R(2,2,:)+v(3,:)*bc%R(3,2,:) !vy
+ vr(3,:) = v(1,:)*bc%R(1,3,:)+v(2,:)*bc%R(2,3,:)+v(3,:)*bc%R(3,3,:) !vz
+
+ endif
+
+end function rotate
+
+!---------------------------------------------------------------------
+! adds a value to a fault parameter inside an area with prescribed shape
+subroutine init_2d_distribution(a,coord,iin,n)
+!JPA refactor: background value shuld be an argument
+
+ real(kind=CUSTOM_REAL), intent(inout) :: a(:)
+ real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
+ integer, intent(in) :: iin,n
+
+ real(kind=CUSTOM_REAL) :: b(size(a))
+ character(len=20) :: shape
+ real(kind=CUSTOM_REAL) :: val,valh, xc, yc, zc, r, l, lx,ly,lz
+ real(kind=CUSTOM_REAL) :: r1(size(a))
+ integer :: i
+
+ NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
+
+ if (n==0) return
+
+ do i=1,n
+ shape = ''
+ xc = 0e0_CUSTOM_REAL
+ yc = 0e0_CUSTOM_REAL
+ zc = 0e0_CUSTOM_REAL
+ r = 0e0_CUSTOM_REAL
+ l = 0e0_CUSTOM_REAL
+ lx = 0e0_CUSTOM_REAL
+ ly = 0e0_CUSTOM_REAL
+ lz = 0e0_CUSTOM_REAL
+ valh = 0e0_CUSTOM_REAL
+
+ read(iin,DIST2D)
+ select case(shape)
+ case ('circle')
+ b = heaviside( r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 ) )
+ case ('circle-exp')
+ r1 = sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 )
+ where(r1<r)
+ b =exp(r1**2/(r1**2 - r**2) )
+ elsewhere
+ b =0._CUSTOM_REAL
+ endwhere
+ case ('ellipse')
+ b = heaviside( 1e0_CUSTOM_REAL - sqrt( (coord(1,:)-xc)**2/lx**2 + (coord(2,:)-yc)**2/ly**2 + (coord(3,:)-zc)**2/lz**2 ) )
+ case ('square')
+ b = heaviside((l/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((l/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((l/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('cylinder')
+ b = heaviside(r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2)) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('rectangle')
+ b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('rectangle-taper')
+ b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ (val + ( coord(3,:) - zc + lz/2._CUSTOM_REAL ) * ((valh-val)/lz))
+ case default
+ stop 'bc_dynflt_3d::init_2d_distribution:: unknown shape'
+ end select
+
+ where (b /= 0e0_CUSTOM_REAL) a = b
+ enddo
+
+end subroutine init_2d_distribution
+
+!---------------------------------------------------------------------
+elemental function heaviside(x)
+
+ real(kind=CUSTOM_REAL), intent(in) :: x
+ real(kind=CUSTOM_REAL) :: heaviside
+
+ if (x>=0e0_CUSTOM_REAL) then
+ heaviside = 1e0_CUSTOM_REAL
+ else
+ heaviside = 0e0_CUSTOM_REAL
+ endif
+
+end function heaviside
+
+
+
+!===============================================================
+! OUTPUTS
+subroutine init_dataT(DataT,coord,nglob,NT,iflt)
+
+ use specfem_par, only : NPROC,myrank
+ ! NT = total number of time steps
+
+ integer, intent(in) :: nglob,NT,iflt
+ real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob)
+ type (dataT_type), intent(out) :: DataT
+
+ real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
+ integer :: i, iglob , IIN, ier, jflt, np, k
+ character(len=70) :: tmpname
+
+ integer :: ipoin, ipoin_local, iproc
+
+ ! 1. read fault output coordinates from user file,
+ ! 2. define iglob: the fault global index of the node nearest to user
+ ! requested coordinate
+
+ IIN = 251
+ open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
+ read(IIN,*) np
+ DataT%npoin =0
+ do i=1,np
+ read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
+ if (jflt==iflt) DataT%npoin = DataT%npoin +1
+ enddo
+ close(IIN)
+
+ if (DataT%npoin == 0) return
+
+ allocate(DataT%iglob(DataT%npoin))
+ allocate(DataT%name(DataT%npoin))
+ allocate(DataT%dist(DataT%npoin)) !Surendra : for parallel fault
+
+ open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
+ if( ier /= 0 ) stop 'error opening FAULT_STATIONS file'
+ read(IIN,*) np
+ k = 0
+ do i=1,np
+ read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
+ if (jflt/=iflt) cycle
+ k = k+1
+ DataT%name(k) = tmpname
+ !search nearest node
+ distkeep = huge(distkeep)
+
+ do iglob=1,nglob
+ dist = sqrt((coord(1,iglob)-xtarget)**2 &
+ + (coord(2,iglob)-ytarget)**2 &
+ + (coord(3,iglob)-ztarget)**2)
+ if (dist < distkeep) then
+ distkeep = dist
+ DataT%iglob(k) = iglob
+ endif
+ enddo
+ DataT%dist(k) = distkeep !Surendra : for parallel fault
+ enddo
+
+ !Surendra : for parallel fault
+ if (PARALLEL_FAULT) then
+ allocate(DataT%islice(DataT%npoin))
+ allocate(DataT%iglob_all(DataT%npoin,0:NPROC-1))
+ allocate(DataT%dist_all(DataT%npoin,0:NPROC-1))
+ call gather_all_i(DataT%iglob,DataT%npoin,DataT%iglob_all,DataT%npoin,NPROC)
+ call gather_all_cr(DataT%dist,DataT%npoin,DataT%dist_all,DataT%npoin,NPROC)
+ if(myrank==0) then
+ do ipoin = 1,DataT%npoin
+ iproc = minloc(DataT%dist_all(ipoin,:), 1) - 1
+ DataT%islice(ipoin) = iproc
+ DataT%iglob(ipoin) = DataT%iglob_all(ipoin,iproc)
+ enddo
+ endif
+
+ call bcast_all_i(DataT%islice,DataT%npoin)
+ call bcast_all_i(DataT%iglob,DataT%npoin)
+
+ DataT%npoin_local = 0
+ do ipoin = 1,DataT%npoin
+ if(myrank == DataT%islice(ipoin)) DataT%npoin_local = DataT%npoin_local + 1
+ enddo
+ allocate(DataT%glob_indx(DataT%npoin_local))
+ do ipoin = 1,DataT%npoin
+ if(myrank == DataT%islice(ipoin)) then
+ ipoin_local = ipoin_local + 1
+ DataT%glob_indx(ipoin_local) = ipoin
+ endif
+ enddo
+ else
+ DataT%npoin_local = DataT%npoin
+ endif !Parallel_fault
+
+
+ ! 3. allocate arrays and set to zero
+ allocate(DataT%d1(NT,DataT%npoin))
+ allocate(DataT%v1(NT,DataT%npoin))
+ allocate(DataT%t1(NT,DataT%npoin))
+ allocate(DataT%d2(NT,DataT%npoin))
+ allocate(DataT%v2(NT,DataT%npoin))
+ allocate(DataT%t2(NT,DataT%npoin))
+ allocate(DataT%t3(NT,DataT%npoin))
+ allocate(DataT%theta(NT,DataT%npoin))
+ DataT%d1 = 0e0_CUSTOM_REAL
+ DataT%v1 = 0e0_CUSTOM_REAL
+ DataT%t1 = 0e0_CUSTOM_REAL
+ DataT%d2 = 0e0_CUSTOM_REAL
+ DataT%v2 = 0e0_CUSTOM_REAL
+ DataT%t2 = 0e0_CUSTOM_REAL
+ DataT%t3 = 0e0_CUSTOM_REAL
+ DataT%theta = 0e0_CUSTOM_REAL
+
+ close(IIN)
+
+end subroutine init_dataT
+
+!---------------------------------------------------------------
+subroutine store_dataT(dataT,d,v,t,theta,itime)
+
+ type(dataT_type), intent(inout) :: dataT
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: theta
+ integer, intent(in) :: itime
+
+ integer :: i,k,ipoin
+
+ do ipoin=1,dataT%npoin_local
+ if (PARALLEL_FAULT) then
+ i = DataT%glob_indx(ipoin)
+ else
+ i = ipoin
+ endif
+ k = dataT%iglob(i)
+ dataT%d1(itime,i) = d(1,k)
+ dataT%d2(itime,i) = d(2,k)
+ dataT%v1(itime,i) = v(1,k)
+ dataT%v2(itime,i) = v(2,k)
+ dataT%t1(itime,i) = t(1,k)
+ dataT%t2(itime,i) = t(2,k)
+ dataT%t3(itime,i) = t(3,k)
+ dataT%theta(itime,i) = theta(k)
+ enddo
+
+end subroutine store_dataT
+
+!-----------------------------------------------------------------
+subroutine write_dataT_all(nt,statelaw)
+
+ integer, intent(in) :: nt
+ integer, intent(in) :: statelaw
+
+ integer :: i
+
+ if (.not.allocated(faults)) return
+ do i = 1,size(faults)
+ call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt,statelaw)
+ enddo
+
+end subroutine write_dataT_all
+
+!------------------------------------------------------------------------
+subroutine SCEC_write_dataT(dataT,DT,NT,statelaw)
+
+ type(dataT_type), intent(in) :: dataT
+ real(kind=CUSTOM_REAL), intent(in) :: DT
+ integer, intent(in) :: NT
+ integer, intent(in) :: statelaw
+
+ integer :: i,k,IOUT,ipoin
+ character :: NTchar*5
+ integer :: today(3), now(3)
+
+ call idate(today) ! today(1)=day, (2)=month, (3)=year
+ call itime(now) ! now(1)=hour, (2)=minute, (3)=second
+
+ IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
+
+ write(NTchar,1) NT
+ NTchar = adjustl(NTchar)
+
+1 format(I5)
+ do ipoin=1,dataT%npoin_local
+ if (PARALLEL_FAULT) then
+ i = DataT%glob_indx(ipoin)
+ else
+ i = ipoin
+ endif
+
+ open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
+ write(IOUT,*) "# problem=TPV104"
+ write(IOUT,*) "# author=Surendra Nadh Somala"
+ write(IOUT,1000) today(2), today(1), today(3), now
+ write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
+ write(IOUT,*) "# code_version=1.1"
+ write(IOUT,*) "# element_size=100 m (*5 GLL nodes)"
+ write(IOUT,*) "# time_step=",DT
+ write(IOUT,*) "# location=",trim(dataT%name(i))
+ write(IOUT,*) "# Column #1 = Time (s)"
+ write(IOUT,*) "# Column #2 = horizontal right-lateral slip (m)"
+ write(IOUT,*) "# Column #3 = horizontal right-lateral slip rate (m/s)"
+ write(IOUT,*) "# Column #4 = horizontal right-lateral shear stress (MPa)"
+ write(IOUT,*) "# Column #5 = vertical up-dip slip (m)"
+ write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
+ write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
+ write(IOUT,*) "# Column #8 = normal stress (MPa)"
+ if (RATE_AND_STATE) write(IOUT,*) "# Column #9 = log10 of state variable (log-seconds)"
+ write(IOUT,*) "#"
+ write(IOUT,*) "# The line below lists the names of the data fields:"
+ if (.not. RATE_AND_STATE) then
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(8(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
+ -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
+ dataT%t3(k,i)/1.0e6_CUSTOM_REAL
+ enddo
+ else
+ if(statelaw == 1) then
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress log-theta"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(9(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
+ -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
+ dataT%t3(k,i)/1.0e6_CUSTOM_REAL, log10(dataT%theta(k,i))
+ enddo
+ else
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress psi"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(9(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
+ -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
+ dataT%t3(k,i)/1.0e6_CUSTOM_REAL, dataT%theta(k,i)
+ enddo
+ endif
+ endif
+ close(IOUT)
+
+
+1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
+
+
+ enddo
+
+end subroutine SCEC_write_dataT
+
+!-------------------------------------------------------------------------------------------------
+
+end module fault_solver_common
More information about the CIG-COMMITS
mailing list