c*********************************************************************** c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- c*********************************************************************** !--------------------------------------------------------------------------------------------------- ! BOP ! ! ROUTINE: DERIVS ! ! INTERFACE: SUBROUTINE DERIVS ( & y, ! supplied - array of stresses and velocities on each cell & dydt) ! returned - array of derivatives of stresses and velocities ! ! USES: USE Green_sshr ! Make data in module "Green_sshr" available ! ! RETURN VALUE: c implicit none ! Not yet converted to using "implicit none" and explicitly typing all variables c :: ! Description implicit real*8 (a-h,o-z) real*4 sa,sab parameter (NF_MAX=15003, NE2_MAX=30006) common /deriv/ tauvp(NF_MAX),ca(NF_MAX),cab(NF_MAX),dx(NF_MAX), . dz(NF_MAX),tauss(NF_MAX),sa(NF_MAX),sab(NF_MAX), . e10,xldis,vpl,xmu,ne,ne2,MMP common /pat/ xc(NF_MAX),zc(NF_MAX),vlim,npat,nb,ibc(20), . iec(NF_MAX),iecinv(NF_MAX) dimension y(1),dydt(1) ! ! NECESSARY PRE-CONDITIONS: ! none ! ! SIDE EFFECTS: ! none ! ! POST_CONDITIONS UPON EXIT ! none ! ! DESCRIPTION: ! Called by MAIN, by RKQC and by RK4 ! Calls depend on value of MMP: ! If MMP = 1 (multipole version) it calls buildtree, sumtree, and forget tree ! If MMP = 0 (not multipole version) it calls: GREENS_FUNCTION ! ! calc derivatives (rt sides) of ode's t isn't in ode's ! shear stress - sshr(src,obs) ! dimension y,dydt by 2*nf in main; eqns for non-bc cells only ! vpl,vel,u,tauss always >0 ! calc dtau/dt from disloc stresses, then dv/dt from fault law ! ! uses radiation damping (quasi-dynamic) approximation to full elastodynamics ! uses slip law version of evolution equation for rate and state friction ! ! EOP !--------------------------------------------------------------------------------------------------- rigid = 3.0d4 ! rigidity in MPa --- NK vs = 3.0d3 ! shear-wave speed in m/s --- NK rigid = rigid * 10.0d0 ! MPa ---> bar NK vs = vs * 1.0d3 * 60.0d0 * 60.0d0 * 24.0d0 * 365.25d0 ! mm/year NK if (MMP .eq. 1) then CALL buildtree(ne, xc(1), zc(1), dx(1), dz(1), y(ne+1)) ! Multipole version end if do 10 n=1,ne ! obs pt: y(n) is tau, y(n+ne) is vel strs = tauvp(n) ! only bc cells slipping k = n + ne ! radiation damping denom = 2.0d0 * vs * ca( n ) + rigid * y( k ) ! radiation damping coefv = 2.0d0 * vs * y( k ) /denom ! radiation damping coeft1 = 2.0d0 * vs * ca( n ) / denom ! radiation damping coeft2 = rigid * y( k ) / denom ! radiation damping if (MMP .eq. 0) then ! no Multipole version do 4 i=1,ne ! this is for source point k=i+ne go to (1,2) igf_test(i,n) ! test to see if this Green's Function has been calculated 1 igf_test(i,n)=2 ! set flag saying have caclulated this one CALL GREENS_FUNCTION(i,n) 2 strs = strs + sshr(i,n)* y(k) ! dtau/dt - does total vel, not diff from vpl 4 continue k=n+ne ! this needed again, since k redefined in inner loop end if if (MMP .eq. 1) then CALL sumtree(xc(n), zc(n), strs) ! Multipole version strs = 2.*xmu*strs ! this done in GREENS_FUNCTION for noMP version; not done in tree.c end if arg = vpl/y(k)+e10 tauss(n) = -cab(n)* log(arg) trm = (y(n)-tauss(n))* y(k)/xldis dydt( n ) = coeft1 * strs - coeft2 * trm ! dtau/dt radiation damping dydt( k ) = coefv * ( strs + trm ) ! dvel/dt radiation damping 10 continue if (MMP .eq. 1) then CALL forgettree ! Multipole version end if return end c*********************************************************************** c---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- c***********************************************************************