[cig-commits] commit: move traction to elastic module; optionally smooth fault patches individually; add option to remove postseismic contribution; add stress output for initial condition.

Mercurial hg at geodynamics.org
Wed Oct 19 15:18:57 PDT 2011


changeset:   30:fc4919855571
user:        Sylvain Barbot <sbarbot at caltech.edu>
date:        Wed Oct 05 19:18:00 2011 -0700
files:       elastic3d.f90 input.f90 relax.f90 types.f90
description:
move traction to elastic module; optionally smooth fault patches individually; add option to remove postseismic contribution; add stress output for initial condition.


diff -r 731e14401888 -r fc4919855571 elastic3d.f90
--- a/elastic3d.f90	Tue Oct 04 01:24:47 2011 -0700
+++ b/elastic3d.f90	Wed Oct 05 19:18:00 2011 -0700
@@ -2064,6 +2064,79 @@ CONTAINS
     END DO
 
   END SUBROUTINE mogisource
+
+  !---------------------------------------------------------------------
+  !> subroutine Traction 
+  !! assigns the traction vector at the surface.
+  !!
+  !! \author sylvain barbot (07-19-07) - original form
+  !---------------------------------------------------------------------
+  SUBROUTINE traction(mu,e,sx1,sx2,dx1,dx2,t,Dt,t3,rate)
+    TYPE(EVENT_STRUC), INTENT(IN) :: e
+    INTEGER, INTENT(IN) :: sx1,sx2
+    REAL*8, INTENT(IN) :: mu,dx1,dx2,t,Dt
+#ifdef ALIGN_DATA
+    REAL*4, DIMENSION(sx1+2,sx2), INTENT(INOUT) :: t3
+#else
+    REAL*4, DIMENSION(sx1,sx2), INTENT(INOUT) :: t3
+#endif
+    LOGICAL, INTENT(IN), OPTIONAL :: rate
+
+    INTEGER :: i,i1,i2,i3,israte
+    REAL*8 :: period,phi,amp,L,W,Lp,Wp,x1,x2,x3,x,y,beta
+
+    REAL*8, PARAMETER :: pi=3.141592653589793115997963468544185161_8
+
+    IF (PRESENT(rate)) THEN
+       israte=rate
+    ELSE
+       israte=.FALSE.
+    END IF
+
+    ! loop over traction sources
+    DO i=1,e%nl
+
+       x=e%l(i)%x
+       y=e%l(i)%y
+
+       L=e%l(i)%length
+       W=e%l(i)%width
+
+       beta=e%l(i)%beta
+
+       ! effective tapered dimensions
+       Lp=L*(1._8+2._8*beta)/2._8
+       Wp=W*(1._8+2._8*beta)/2._8
+
+       i3=1
+       DO i2=1,sx2
+          DO i1=1,sx1
+             CALL shiftedcoordinates(i1,i2,i3,sx1,sx2,1, &
+                                     dx1,dx2,1.d8,x1,x2,x3)
+
+             IF ((ABS(x1-x).GT.MAX(Lp,Wp)).OR.(ABS(x2-y).GT.MAX(Lp,Wp))) CYCLE
+
+             amp=omega((x1-x)/L,beta)* &
+                 omega((x2-y)/W,beta)* &
+                 mu*e%l(i)%slip
+
+             IF (israte) THEN
+                ! surface tractions rate
+                period=e%l(i)%period
+                phi=e%l(i)%phase
+
+                t3(i1,i2)=t3(i1,i2)-amp*(sin(2*pi*(t+Dt)/period+phi)-sin(2*pi*t/period+phi))
+             ELSE
+                IF (e%l(i)%period .LE. 0) THEN
+                   ! surface tractions
+                   t3(i1,i2)=t3(i1,i2)-amp
+                END IF
+             END IF
+          END DO
+       END DO
+    END DO
+
+  END SUBROUTINE traction
 
   !---------------------------------------------------------------------
   !! function MomentDensityShear
diff -r 731e14401888 -r fc4919855571 input.f90
--- a/input.f90	Tue Oct 04 01:24:47 2011 -0700
+++ b/input.f90	Wed Oct 05 19:18:00 2011 -0700
@@ -57,7 +57,7 @@ CONTAINS
 !$  INTEGER :: omp_get_num_procs,omp_get_max_threads
     REAL*8 :: dummy,dum1,dum2
     REAL*8 :: minlength,minwidth
-    TYPE(OPTION_S) :: opts(8)
+    TYPE(OPTION_S) :: opts(9)
 
     INTEGER :: k,iostatus,i,e
 
@@ -70,13 +70,14 @@ CONTAINS
 
     ! parse the command line for options
     opts(1)=OPTION_S("no-proj-output",.FALSE.,CHAR(20))
-    opts(2)=OPTION_S("no-txt-output",.FALSE.,CHAR(21))
-    opts(3)=OPTION_S("no-vtk-output",.FALSE.,CHAR(22))
-    opts(4)=OPTION_S("no-grd-output",.FALSE.,CHAR(23))
-    opts(5)=OPTION_S("no-xyz-output",.FALSE.,CHAR(24))
-    opts(5)=OPTION_S("no-stress-output",.FALSE.,CHAR(25))
-    opts(6)=OPTION_S("dry-run",.FALSE.,CHAR(26))
-    opts(7)=OPTION_S("help",.FALSE.,'h')
+    opts(2)=OPTION_S("no-relax-output",.FALSE.,CHAR(21))
+    opts(3)=OPTION_S("no-txt-output",.FALSE.,CHAR(22))
+    opts(4)=OPTION_S("no-vtk-output",.FALSE.,CHAR(23))
+    opts(5)=OPTION_S("no-grd-output",.FALSE.,CHAR(24))
+    opts(6)=OPTION_S("no-xyz-output",.FALSE.,CHAR(25))
+    opts(7)=OPTION_S("no-stress-output",.FALSE.,CHAR(26))
+    opts(8)=OPTION_S("dry-run",.FALSE.,CHAR(27))
+    opts(9)=OPTION_S("help",.FALSE.,'h')
 
     DO
        ch=getopt("h",opts)
@@ -87,21 +88,24 @@ CONTAINS
           ! option no-proj-output
           in%isoutputproj=.FALSE.
        CASE(CHAR(21))
+          ! option no-relax-output
+          in%isoutputrelax=.FALSE.
+       CASE(CHAR(22))
           ! option no-txt-output
           in%isoutputtxt=.FALSE.
-       CASE(CHAR(22))
+       CASE(CHAR(23))
           ! option no-vtk-output
           in%isoutputvtk=.FALSE.
-       CASE(CHAR(23))
+       CASE(CHAR(24))
           ! option no-grd-output
           in%isoutputgrd=.FALSE.
-       CASE(CHAR(24))
+       CASE(CHAR(25))
           ! option no-xyz-output
           in%isoutputxyz=.FALSE.
-       CASE(CHAR(25))
+       CASE(CHAR(26))
           ! option stress output
           in%isoutputstress=.FALSE.
-       CASE(CHAR(26))
+       CASE(CHAR(27))
           ! option dry-run
           in%isdryrun=.TRUE.
        CASE('h')
@@ -120,21 +124,24 @@ CONTAINS
 
     IF (in%ishelp) THEN
        PRINT '("usage:")'
-       PRINT '("relax [-h] [--dry-run] [--help] [--no-grd-output] [--no-stress-output]")' 
-       PRINT '("      [--no-txt-output] [--no-vtk-output] [--no-xyz-output]")'
+       PRINT '("relax [-h] [--dry-run] [--help] [--no-grd-output] [--no-proj-output]")' 
+       PRINT '("      [--no-relax-output] [--no-stress-output] [--no-txt-output]")'
+       PRINT '("      [--no-vtk-output] [--no-xyz-output]")'
        PRINT '("")'
        PRINT '("options:")'
        PRINT '("   -h                 prints this message and aborts calculation")'
        PRINT '("   --dry-run          abort calculation, only output geometry")'
        PRINT '("   --help             prints this message and aborts calculation")'
        PRINT '("   --no-grd-output    cancel output in GMT grd binary format")'
+       PRINT '("   --no-proj-output   cancel output in geographic projection")'
+       PRINT '("   --no-relax-output  cancel output of the postseismic contribution")'
        PRINT '("   --no-stress-output cancel output of stress tensor in any format")'
        PRINT '("   --no-txt-output    cancel output in text format")'
        PRINT '("   --no-vtk-output    cancel output in Paraview VTK format")'
        PRINT '("   --no-xyz-output    cancel output in GMT xyz format")'
        PRINT '("")'
        PRINT '("description:")'
-       PRINT '("   Evaluates the deformation due to fault slip, surfacial loading")'
+       PRINT '("   Evaluates the deformation due to fault slip, surface loading")'
        PRINT '("   or inflation and the time series of postseismic relaxation")'
        PRINT '("   that follows due to fault creep or viscoelastic flow.")'
        RETURN
@@ -777,6 +784,7 @@ CONTAINS
                in%inter%s(k)%x,in%inter%s(k)%y,in%inter%s(k)%z, &
                in%inter%s(k)%length,in%inter%s(k)%width, &
                in%inter%s(k)%strike,in%inter%s(k)%dip,in%inter%s(k)%rake
+
           ! copy the input format for display
           in%inter%sc(k)=in%inter%s(k)
              
@@ -933,19 +941,44 @@ CONTAINS
           PRINT 2000
           DO k=1,in%events(e)%ns
              CALL getdata(iunit,dataline)
-             READ (dataline,*) i,in%events(e)%s(k)%slip, &
+             READ (dataline,*,IOSTAT=iostatus) i,in%events(e)%s(k)%slip, &
                   in%events(e)%s(k)%x,in%events(e)%s(k)%y,in%events(e)%s(k)%z, &
                   in%events(e)%s(k)%length,in%events(e)%s(k)%width, &
-                  in%events(e)%s(k)%strike,in%events(e)%s(k)%dip,in%events(e)%s(k)%rake
+                  in%events(e)%s(k)%strike,in%events(e)%s(k)%dip,in%events(e)%s(k)%rake, &
+                  in%events(e)%s(k)%beta
+
+             SELECT CASE(iostatus)
+             CASE (1:)
+                WRITE_DEBUG_INFO
+                WRITE (0,'("invalid shear source definition at line")')
+                WRITE (0,'(a)') dataline
+                STOP 1
+             CASE (0)
+                IF (in%events(e)%s(k)%beta.GT.0.5d8) STOP "invalid smoothing parameter (beta)."
+             CASE (:-1)
+                ! use default value for smoothing
+                in%events(e)%s(k)%beta=in%beta
+             END SELECT
+
              ! copy the input format for display
              in%events(e)%sc(k)=in%events(e)%s(k)
              
-             PRINT '(I3.3,4ES9.2E1,2ES8.2E1,f7.1,f6.1,f7.1)',i, &
-                  in%events(e)%sc(k)%slip,&
-                  in%events(e)%sc(k)%x,in%events(e)%sc(k)%y,in%events(e)%sc(k)%z, &
-                  in%events(e)%sc(k)%length,in%events(e)%sc(k)%width, &
-                  in%events(e)%sc(k)%strike,in%events(e)%sc(k)%dip, &
-                  in%events(e)%sc(k)%rake
+             IF (iostatus.NE.0) THEN
+                PRINT '(I3.3,4ES9.2E1,2ES8.2E1,f7.1,f6.1,f7.1)',i, &
+                     in%events(e)%sc(k)%slip,&
+                     in%events(e)%sc(k)%x,in%events(e)%sc(k)%y,in%events(e)%sc(k)%z, &
+                     in%events(e)%sc(k)%length,in%events(e)%sc(k)%width, &
+                     in%events(e)%sc(k)%strike,in%events(e)%sc(k)%dip, &
+                     in%events(e)%sc(k)%rake
+             ELSE
+                ! print the smoothing value for this patch
+                PRINT '(I3.3,4ES9.2E1,2ES8.2E1,f7.1,f6.1,f7.1,f6.1)',i, &
+                     in%events(e)%sc(k)%slip,&
+                     in%events(e)%sc(k)%x,in%events(e)%sc(k)%y,in%events(e)%sc(k)%z, &
+                     in%events(e)%sc(k)%length,in%events(e)%sc(k)%width, &
+                     in%events(e)%sc(k)%strike,in%events(e)%sc(k)%dip, &
+                     in%events(e)%sc(k)%rake,in%events(e)%sc(k)%beta
+             END IF
              
              IF (i .ne. k) THEN
                 WRITE_DEBUG_INFO
@@ -1100,20 +1133,48 @@ CONTAINS
           PRINT 2000
           PRINT '(a)',"t3 in units of force/surface/rigidity, positive down"
           PRINT '(a)',"T>0 for t3 sin(2pi/T+phi), T<=0 for t3 H(t)"
-          PRINT '(a)',"no.       xs       ys       t3        T      phi"
+          PRINT '(a)',"no.       xs       ys   length    width       t3        T      phi"
           PRINT 2000
           DO k=1,in%events(e)%nl
              CALL getdata(iunit,dataline)
-             READ  (dataline,*) i, &
-                  in%events(e)%l(k)%x,in%events(e)%l(k)%y,in%events(e)%l(k)%slip, &
-                  in%events(e)%l(k)%period,in%events(e)%l(k)%phase
+             READ  (dataline,*,IOSTAT=iostatus) i, &
+                  in%events(e)%l(k)%x,in%events(e)%l(k)%y, &
+                  in%events(e)%l(k)%length,in%events(e)%l(k)%width, &
+                  in%events(e)%l(k)%slip, &
+                  in%events(e)%l(k)%period,in%events(e)%l(k)%phase, &
+                  in%events(e)%l(k)%beta
+             
+             SELECT CASE(iostatus)
+             CASE (1:)
+                WRITE_DEBUG_INFO
+                WRITE (0,'("invalid surface load definition at line")')
+                WRITE (0,'(a)') dataline
+                STOP 1
+             CASE (0)
+                IF (in%events(e)%l(k)%beta.GT.0.5d8) STOP "invalid smoothing parameter beta."
+             CASE (:-1)
+                ! use default value for smoothing
+                in%events(e)%l(k)%beta=in%beta
+             END SELECT
+
              ! copy the input format for display
              in%events(e)%lc(k)=in%events(e)%l(k)
-             
-             PRINT '(I3.3,6ES9.2E1)',k, &
-                  in%events(e)%lc(k)%x,in%events(e)%lc(k)%y,in%events(e)%lc(k)%slip, &
-                  in%events(e)%lc(k)%period,in%events(e)%lc(k)%phase
-             
+
+             IF (iostatus.EQ.0) THEN
+                PRINT '(I3.3,9ES9.2E1)',k, &
+                     in%events(e)%lc(k)%x,in%events(e)%lc(k)%y, &
+                     in%events(e)%lc(k)%length,in%events(e)%lc(k)%width, &
+                     in%events(e)%lc(k)%slip, &
+                     in%events(e)%lc(k)%period,in%events(e)%lc(k)%phase, &
+                     in%events(e)%lc(k)%beta
+             ELSE
+                PRINT '(I3.3,8ES9.2E1)',k, &
+                     in%events(e)%lc(k)%x,in%events(e)%lc(k)%y, &
+                     in%events(e)%lc(k)%length,in%events(e)%lc(k)%width, &
+                     in%events(e)%lc(k)%slip, &
+                     in%events(e)%lc(k)%period,in%events(e)%lc(k)%phase
+             END IF
+
              IF (i .NE. k) THEN
                 PRINT *, "error in input file: source index misfit"
                 STOP 1
diff -r 731e14401888 -r fc4919855571 relax.f90
--- a/relax.f90	Tue Oct 04 01:24:47 2011 -0700
+++ b/relax.f90	Wed Oct 05 19:18:00 2011 -0700
@@ -243,11 +243,15 @@ PROGRAM relax
   ! allocate memory
   ALLOCATE (v1(in%sx1+2,in%sx2,in%sx3),v2(in%sx1+2,in%sx2,in%sx3),v3(in%sx1+2,in%sx2,in%sx3), &
             u1(in%sx1+2,in%sx2,in%sx3/2),u2(in%sx1+2,in%sx2,in%sx3/2),u3(in%sx1+2,in%sx2,in%sx3/2), &
-            inter1(in%sx1+2,in%sx2,2),inter2(in%sx1+2,in%sx2,2),inter3(in%sx1+2,in%sx2,2), &
             tau(in%sx1,in%sx2,in%sx3/2),sig(in%sx1,in%sx2,in%sx3/2),gamma(in%sx1+2,in%sx2,in%sx3/2), &
             t1(in%sx1+2,in%sx2),t2(in%sx1+2,in%sx2),t3(in%sx1+2,in%sx2), &
             STAT=iostatus)
   IF (iostatus>0) STOP "could not allocate memory"
+
+  IF (in%isoutputrelax) THEN
+     ALLOCATE(inter1(in%sx1+2,in%sx2,2),inter2(in%sx1+2,in%sx2,2),inter3(in%sx1+2,in%sx2,2),STAT=iostatus)
+     IF (iostatus>0) STOP "could not allocate memory for postseismic displacement"
+  END IF
 
   v1=0;v2=0;v3=0;u1=0;u2=0;u3=0;gamma=0;t1=0;t2=0;t3=0
   CALL tensorfieldadd(tau,tau,in%sx1,in%sx2,in%sx3/2,c1=0._4,c2=0._4)
@@ -298,7 +302,7 @@ PROGRAM relax
 
   ! sources
   CALL dislocations(in%events(e),in%lambda,in%mu,in%beta,in%sx1,in%sx2,in%sx3, &
-       in%dx1,in%dx2,in%dx3,v1,v2,v3,t1,t2,t3,tau)
+                    in%dx1,in%dx2,in%dx3,v1,v2,v3,t1,t2,t3,tau)
   CALL traction(in%mu,in%events(e),in%sx1,in%sx2,in%dx1,in%dx2,t,0.d0,t3)
   
   PRINT '("coseismic event ",I3.3)', e
@@ -353,8 +357,10 @@ PROGRAM relax
 #ifdef GRD
   IF (in%isoutputgrd) THEN
      CALL exportgrd(u1,u2,u3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,in%x0,in%y0,in%wdir,0)
-     CALL exportgrd(inter1,inter2,inter3,in%sx1,in%sx2,in%sx3/2, &
-          in%dx1,in%dx2,in%dx3,0._8,in%x0,in%y0,in%wdir,0,convention=2)
+     IF (in%isoutputrelax) THEN
+        CALL exportgrd(inter1,inter2,inter3,in%sx1,in%sx2,in%sx3/2, &
+                       in%dx1,in%dx2,in%dx3,0._8,in%x0,in%y0,in%wdir,0,convention=2)
+     END IF
   END IF
 #endif
 #ifdef PROJ
@@ -605,9 +611,11 @@ PROGRAM relax
      CALL tensorfieldadd(tau,moment,in%sx1,in%sx2,in%sx3/2,c2=REAL(Dt))
      
      ! keep track of the viscoelastic contribution alone
-     CALL sliceadd(inter1(:,:,1),v1,in%sx1+2,in%sx2,in%sx3,int(in%oz/in%dx3)+1,c2=REAL(Dt))
-     CALL sliceadd(inter2(:,:,1),v2,in%sx1+2,in%sx2,in%sx3,int(in%oz/in%dx3)+1,c2=REAL(Dt))
-     CALL sliceadd(inter3(:,:,1),v3,in%sx1+2,in%sx2,in%sx3,int(in%oz/in%dx3)+1,c2=REAL(Dt))
+     IF (in%isoutputrelax) THEN
+        CALL sliceadd(inter1(:,:,1),v1,in%sx1+2,in%sx2,in%sx3,int(in%oz/in%dx3)+1,c2=REAL(Dt))
+        CALL sliceadd(inter2(:,:,1),v2,in%sx1+2,in%sx2,in%sx3,int(in%oz/in%dx3)+1,c2=REAL(Dt))
+        CALL sliceadd(inter3(:,:,1),v3,in%sx1+2,in%sx2,in%sx3,int(in%oz/in%dx3)+1,c2=REAL(Dt))
+     END IF
 
      ! time increment
      t=t+Dt
@@ -665,22 +673,28 @@ PROGRAM relax
 #ifdef XYZ
         IF (in%isoutputxyz) THEN
            CALL exportxyz(u1,u2,u3,in%sx1,in%sx2,in%sx3/2,in%oz,in%dx1,in%dx2,in%dx3,i,in%wdir)
-           !CALL exportxyz(inter1,inter2,inter3,in%sx1,in%sx2,in%sx3/2,0.0_8,in%dx1,in%dx2,in%dx3,i,in%wdir)
+           IF (in%isoutputrelax) THEN
+              !CALL exportxyz(inter1,inter2,inter3,in%sx1,in%sx2,in%sx3/2,0.0_8,in%dx1,in%dx2,in%dx3,i,in%wdir)
+           END IF
         END IF
 #endif
         CALL exporteigenstrain(gamma,in%nop,in%op,in%x0,in%y0,in%dx1,in%dx2,in%dx3,in%sx1,in%sx2,in%sx3/2,in%wdir,oi)
 #ifdef GRD
         IF (in%isoutputgrd) THEN
-           CALL exportgrd(inter1,inter2,inter3,in%sx1,in%sx2,in%sx3/2, &
-                          in%dx1,in%dx2,in%dx3,0._8,in%x0,in%y0,in%wdir,oi,convention=2)
+           IF (in%isoutputrelax) THEN
+              CALL exportgrd(inter1,inter2,inter3,in%sx1,in%sx2,in%sx3/2, &
+                             in%dx1,in%dx2,in%dx3,0._8,in%x0,in%y0,in%wdir,oi,convention=2)
+           END IF
            CALL exportgrd(u1,u2,u3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,in%x0,in%y0,in%wdir,oi)
         END IF
 #endif
 #ifdef PROJ
         IF (in%isoutputproj) THEN
-           CALL exportproj(inter1,inter2,inter3,in%sx1,in%sx2,in%sx3/2, &
-                           in%dx1,in%dx2,in%dx3,in%oz,in%x0,in%y0, &
-                           in%lon0,in%lat0,in%zone,in%umult,in%wdir,oi,convention=2)
+           IF (in%isoutputrelax) THEN
+              CALL exportproj(inter1,inter2,inter3,in%sx1,in%sx2,in%sx3/2, &
+                              in%dx1,in%dx2,in%dx3,in%oz,in%x0,in%y0, &
+                              in%lon0,in%lat0,in%zone,in%umult,in%wdir,oi,convention=2)
+           END IF
            CALL exportproj(u1,u2,u3,in%sx1,in%sx2,in%sx3/2,in%dx1,in%dx2,in%dx3,in%oz,in%x0,in%y0, &
                            in%lon0,in%lat0,in%zone,in%umult,in%wdir,oi)
         END IF
@@ -803,54 +817,6 @@ 1200 FORMAT ("--------------------------
 
 CONTAINS
 
-  !---------------------------------------------------------------------
-  !> subroutine Traction 
-  !! assigns the traction vector at the surface.
-  !!
-  !! \author sylvain barbot (07-19-07) - original form
-  !---------------------------------------------------------------------
-  SUBROUTINE traction(mu,e,sx1,sx2,dx1,dx2,t,Dt,t3,rate)
-    TYPE(EVENT_STRUC), INTENT(IN) :: e
-    INTEGER, INTENT(IN) :: sx1,sx2
-    REAL*8, INTENT(IN) :: mu,dx1,dx2,t,Dt
-#ifdef ALIGN_DATA
-    REAL*4, DIMENSION(sx1+2,sx2), INTENT(INOUT) :: t3
-#else
-    REAL*4, DIMENSION(sx1,sx2), INTENT(INOUT) :: t3
-#endif
-    LOGICAL, INTENT(IN), OPTIONAL :: rate
-
-    INTEGER :: i,i1,i2,i3,israte
-    REAL*8 :: period,phi
-    REAL*8, PARAMETER :: pi=3.141592653589793115997963468544185161_8
-
-    IF (PRESENT(rate)) THEN
-       israte=rate
-    ELSE
-       israte=.FALSE.
-    END IF
-
-    DO i=1,e%nl
-       CALL shiftedindex(e%l(i)%x,e%l(i)%y,0._8,sx1,sx2,1,dx1,dx2,1._8,i1,i2,i3)
-
-       IF (e%l(i)%period .LE. 0) THEN
-          ! surface tractions
-          t3(i1,i2)=t3(i1,i2)-mu*e%l(i)%slip
-
-       ELSE
-          ! velocity for t>0 only
-          IF (israte) THEN
-             ! surface tractions rate
-             period=e%l(i)%period
-             phi=e%l(i)%phase
-
-             t3(i1,i2)=t3(i1,i2)-mu*e%l(i)%slip*(sin(2*pi*(t+Dt)/period+phi)-sin(2*pi*t/period+phi))
-          END IF
-       END IF
-    END DO
-             
-  END SUBROUTINE traction
-
   !--------------------------------------------------------------------
   !> subroutine dislocations
   !! assigns equivalent body forces or moment density to simulate
@@ -889,7 +855,7 @@ CONTAINS
                event%s(i)%x,event%s(i)%y,event%s(i)%z, &
                event%s(i)%width,event%s(i)%length, &
                event%s(i)%strike,event%s(i)%dip,event%s(i)%rake, &
-               beta,sx1,sx2,sx3,dx1,dx2,dx3,v1,v2,v3,t1,t2,t3)
+               event%s(i)%beta,sx1,sx2,sx3,dx1,dx2,dx3,v1,v2,v3,t1,t2,t3)
        END DO
     ELSE
        ! forcing term in moment density
@@ -898,7 +864,7 @@ CONTAINS
                event%s(i)%x,event%s(i)%y,event%s(i)%z, &
                event%s(i)%width,event%s(i)%length, &
                event%s(i)%strike,event%s(i)%dip,event%s(i)%rake, &
-               beta,sx1,sx2,sx3/2,dx1,dx2,dx3,eigenstress)
+               event%s(i)%beta,sx1,sx2,sx3/2,dx1,dx2,dx3,eigenstress)
        END DO
     END IF
 
@@ -908,7 +874,7 @@ CONTAINS
             event%s(i)%x,event%s(i)%y,event%s(i)%z, &
             event%s(i)%width,event%s(i)%length, &
             event%s(i)%strike,event%s(i)%dip,event%s(i)%rake, &
-            beta,sx1,sx2,sx3/2,dx1,dx2,dx3,tau)
+            event%s(i)%beta,sx1,sx2,sx3/2,dx1,dx2,dx3,tau)
     END DO
     
     ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
diff -r 731e14401888 -r fc4919855571 types.f90
--- a/types.f90	Tue Oct 04 01:24:47 2011 -0700
+++ b/types.f90	Wed Oct 05 19:18:00 2011 -0700
@@ -23,7 +23,7 @@ MODULE types
 
   TYPE SOURCE_STRUCT
      SEQUENCE
-     REAL*8 :: slip,x,y,z,width,length,strike,dip,rake,period,phase
+     REAL*8 :: slip,x,y,z,width,length,strike,dip,rake,period,phase,beta
   END TYPE SOURCE_STRUCT
 
   TYPE PLANE_STRUCT
@@ -187,6 +187,7 @@ MODULE types
 
      ! overrides output to formats
      LOGICAL :: isoutputproj=.TRUE.
+     LOGICAL :: isoutputrelax=.TRUE.
      LOGICAL :: isoutputtxt=.TRUE.
      LOGICAL :: isoutputvtk=.TRUE.
      LOGICAL :: isoutputgrd=.TRUE.



More information about the CIG-COMMITS mailing list