[cig-commits] r19186 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db/src: . devel
percygalvez at geodynamics.org
percygalvez at geodynamics.org
Fri Nov 11 04:03:29 PST 2011
Author: percygalvez
Date: 2011-11-11 04:03:29 -0800 (Fri, 11 Nov 2011)
New Revision: 19186
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/devel/fault_solver.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Faut solver with fault initial parameters output
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/devel/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/devel/fault_solver.f90 2011-11-11 03:12:40 UTC (rev 19185)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/devel/fault_solver.f90 2011-11-11 12:03:29 UTC (rev 19186)
@@ -127,9 +127,9 @@
dummy_idfault = 0
- open(unit=IIN_PAR,file='DATA/FAULT/Par_file_faults.in',status='old',iostat=ier)
+ open(unit=IIN_PAR,file='DATA/FAULT/Par_file_faults',status='old',iostat=ier)
if( ier /= 0 ) then
- write(6,*) 'File Par_file_faults.in not found: assume no faults'
+ write(6,*) 'File Par_file_faults not found: assume no faults'
close(IIN_PAR)
return
endif
@@ -329,6 +329,10 @@
call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
call init_dataXZ(bc%dataXZ,bc,bc%nglob)
+!---- Initial Fault Parameters -------------
+ call write_init_fault_par(bc,iflt)
+
+
end subroutine init_one_fault
!---------------------------------------------------------------------
@@ -399,14 +403,20 @@
read(iin,DIST2D)
select case(shape)
case ('circle')
- b = heaviside( r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 ) )
+ b = heaviside(r-sqrt((coord(1,:)-xc)**2 + &
+ (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2)) * val
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 ) )
+ b = heaviside(1e0_CUSTOM_REAL-sqrt((coord(1,:)-xc)**2/lx**2 + &
+ (coord(2,:)-yc)**2/ly**2 + (coord(3,:)-zc)**2/lz**2)) * val
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 ('cilinder')
+ 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) * &
@@ -421,7 +431,7 @@
stop 'bc_dynflt_3d::init_2d_distribution:: unknown shape'
end select
- where (b /= 0) a = b
+ where (b /= 0e0_CUSTOM_REAL) a = b
enddo
end subroutine init_2d_distribution
@@ -524,8 +534,6 @@
! Solve for shear stress
tnorm = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
-! TO DO :
-! To avoid division by cero (temporal)
tnorm = max(tnorm,1e0_CUSTOM_REAL)
t1 = T(1,:)/tnorm
t2 = T(2,:)/tnorm
@@ -689,7 +697,7 @@
! requested coordinate
IIN = 251
- open(IIN,file='DATA/FAULT/FAULT_STATIONS.in',status='old',action='read',iostat=ier)
+ open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
read(IIN,*) np
DataT%npoin =0
do i=1,np
@@ -703,7 +711,7 @@
allocate(DataT%iglob(DataT%npoin))
allocate(DataT%name(DataT%npoin))
- open(IIN,file='DATA/FAULT/FAULT_STATIONS.in',status='old',action='read',iostat=ier)
+ 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
@@ -975,5 +983,39 @@
end subroutine write_dataXZ
+!---------------------------------------------------------------
+ subroutine write_init_fault_par(bc,iflt)
+
+ type(bc_dynflt_type), intent(inout) :: bc
+ integer, intent(in) :: iflt
+
+ character(len=70) :: filename
+
+
+ write(filename,"('OUTPUT_FILES/Fault_init_par','_F',I0,'.bin')") iflt
+! open(unit=IOUT, file= trim(filename), status='replace', form='formatted',action='write')
+! NOTE : It had to be adopted formatted output to avoid conflicts readings with different
+! compilers.
+
+! write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
+
+ open(unit=IOUT, file= trim(filename), status='replace', form='unformatted',action='write')
+
+ write(IOUT) bc%coord(1,:)
+ write(IOUT) bc%coord(2,:)
+ write(IOUT) bc%coord(3,:)
+ write(IOUT) bc%swf%mus
+ write(IOUT) bc%swf%mud
+ write(IOUT) bc%swf%Dc
+ write(IOUT) bc%T0(1,:)
+ write(IOUT) bc%T0(2,:)
+ write(IOUT) bc%T0(3,:)
+
+ close(IOUT)
+
+ end subroutine write_init_fault_par
+
+!---
+
end module fault_solver
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2011-11-11 03:12:40 UTC (rev 19185)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2011-11-11 12:03:29 UTC (rev 19186)
@@ -329,6 +329,10 @@
call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
call init_dataXZ(bc%dataXZ,bc,bc%nglob)
+!---- Initial Fault Parameters -------------
+ call write_init_fault_par(bc,iflt)
+
+
end subroutine init_one_fault
!---------------------------------------------------------------------
@@ -399,14 +403,20 @@
read(iin,DIST2D)
select case(shape)
case ('circle')
- b = heaviside( r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 ) )
+ b = heaviside(r-sqrt((coord(1,:)-xc)**2 + &
+ (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2)) * val
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 ) )
+ b = heaviside(1e0_CUSTOM_REAL-sqrt((coord(1,:)-xc)**2/lx**2 + &
+ (coord(2,:)-yc)**2/ly**2 + (coord(3,:)-zc)**2/lz**2)) * val
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 ('cilinder')
+ 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) * &
@@ -973,5 +983,39 @@
end subroutine write_dataXZ
+!---------------------------------------------------------------
+ subroutine write_init_fault_par(bc,iflt)
+
+ type(bc_dynflt_type), intent(inout) :: bc
+ integer, intent(in) :: iflt
+
+ character(len=70) :: filename
+
+
+ write(filename,"('OUTPUT_FILES/Fault_init_par','_F',I0,'.bin')") iflt
+! open(unit=IOUT, file= trim(filename), status='replace', form='formatted',action='write')
+! NOTE : It had to be adopted formatted output to avoid conflicts readings with different
+! compilers.
+
+! write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
+
+ open(unit=IOUT, file= trim(filename), status='replace', form='unformatted',action='write')
+
+ write(IOUT) bc%coord(1,:)
+ write(IOUT) bc%coord(2,:)
+ write(IOUT) bc%coord(3,:)
+ write(IOUT) bc%swf%mus
+ write(IOUT) bc%swf%mud
+ write(IOUT) bc%swf%Dc
+ write(IOUT) bc%T0(1,:)
+ write(IOUT) bc%T0(2,:)
+ write(IOUT) bc%T0(3,:)
+
+ close(IOUT)
+
+ end subroutine write_init_fault_par
+
+!---
+
end module fault_solver
More information about the CIG-COMMITS
mailing list