VUAMT 学习笔记

1. 放一个晶体塑性代码

C #=======================================================================
C		  							               _______
C                 \    /  |    |   |\  /|    /\       |
C                  \  /   |    |   | \/ |   /__\      |
C                   \/    |____|   |    |  /    \     |
C
C     This subroutine is the MAIN MATERIAL SUBROUTINE: 
C     It updates: stresses, internal variables (SDV array), and energies 
C =======================================================================
      subroutine vumat (
C Read only -
     1     nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     2     stepTime, totalTime, dt, cmname, coordMp, charLength,
     3     props, density, strainInc, relSpinInc,
     4     tempOld, stretchOld, defgradOld, fieldOld,
     5     stressOld, stateOld, enerInternOld, enerInelasOld,
     6     tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
     7     stressNew, stateNew, enerInternNew, enerInelasNew )
      include 'vaba_param.inc'
      include 'DataInitialize.inc'
      include 'ElasticStarter.inc'
C
C #=====================================================================
C		 #   LOADING MATERIAL DATA FROM INPUT FILE  #
C		(independent of Element), (nprops = 7*nslip+27)
C         Note: Not all data in props is loaded here. 
C          Some will be loaded when needed
C #=====================================================================
	  nslip        = int(props(1))											   
      e            = props(2)
      xnu          = props(3)
      yield        = props(4)
	  g_source     = props(10)
	  g_immob      = props(11)
	  g_minter     = props(12)
	  g_recov      = props(13)
	  Cp           = props(14)	
	  tempR        = props(15)												
	  H_K	       = props(16)  	                                        
	  xi           = props(17)												
	  Chi          = props(18)												
	  EM_factor    = props(20)		
	  b_v(1:nslip) = props(21:20+nslip)
C 
C Derived material constants	  
	  g        = e/(two*(one + xnu))								   
	  twomu    = two*g												   
	  bulk     = e/three/( one - two * xnu )										   
C
C Note H/K is assumed constant, so temp must stay at tempR
	  thermal_coef = H_K/tempR     
C		
C Derived Pointers/Flags
	  ibeginN  = 20 + nslip           
	  ibeginS  = nslip*3 + ibeginN     
	  		    
C #=====================================================================
C							         # WORKING WITH ELEMENT "k" #
C #=====================================================================
	   do k = 1,nblock
C Load Internal variables independent of slip system
          Psi21   = stateOld( k,1)*pi/180.d0
          Psi32   = stateOld( k,2)*pi/180.d0
          Psi13   = stateOld( k,3)*pi/180.d0
          gamma   = stateOld( k,4)
          temp    = stateOld( k,6) + tempR
C Load Internal variables dependent on slip system
          do j = 1, nslip													
			 tau(j)       = stateOld( k,8  + (j-1)*10 )				    
			 gdot(j)      = stateOld( k,9  + (j-1)*10 )				
			 den_im(j)    = stateOld( k,10 + (j-1)*10 )				
			 den_m(j)     = stateOld( k,11 + (j-1)*10 )				
			 slip_n(j,1)  = stateOld( k,12 + (j-1)*10 )	 			
			 slip_n(j,2)  = stateOld( k,13 + (j-1)*10 )				
			 slip_n(j,3)  = stateOld( k,14 + (j-1)*10 )				
			 slip_s(j,1)  = stateOld( k,15 + (j-1)*10 )				
			 slip_s(j,2)  = stateOld( k,16 + (j-1)*10 )				
			 slip_s(j,3)  = stateOld( k,17 + (j-1)*10 )				
          end do
C   
C Compute derived slip-specific mechanical constants	  
		  do j = 1, nslip
				 recip_m(j)  = one / props(6)
				 ref_gdot(j) = props(7)																		 				 
		  end do
C														
C    -------------------------------------------------------------------
C	    # IF APPLICABLE: REPLACE ZEROED VARIABLES WITH INITIAL DATA #              	   
C    -------------------------------------------------------------------

          if (dt .eq. totalTime) then
              do i = 1, nslip																	   
				   den_im(i) = props(8)
				   den_m (i) = props(9)
				   do j =  1, 3
					  slip_n(i,j) = props( ibeginN + 3*(i-1) + j )
					  slip_s(i,j) = props( ibeginS + 3*(i-1) + j )
				   enddo
              enddo					  
          endif	
C					   -------------------------------------------------
C				        #  COMPUTE TOTAL SPIN and RATE OF DEFORMATION  #
C                       #      COMPUTE ROTATION FROM F = RU            # 
C                       #      (Rnew is the Transpose of R)            #
C					   -------------------------------------------------

C Compute Spin, Dij and Rotation tensors
         call getSpinDij(dt,defgradNew(k,:),defgradOld(k,:),Spin,Dij)
	     call getRot(defgradNew(k,:),stretchNew(k,:),Rnew)       	
C				   		    	   	
C Deviatoric Dij tensor
		 traceDij = Dij(1) + Dij(2) + Dij(3) 
		 Dij_dev(1) = Dij(1) - traceDij/3.d0
		 Dij_dev(2) = Dij(2) - traceDij/3.d0
		 Dij_dev(3) = Dij(3) - traceDij/3.d0
		 Dij_dev(4) = Dij(4)
		 Dij_dev(5) = Dij(5)
		 Dij_dev(6) = Dij(6)
C
C					   -------------------------------------------------
C						#               COMPUTE Pij & Wij             #
C					   -------------------------------------------------				
C SYMMETRIC SCHMID TENSOR (EXPRESSED W.R.T GLOBAL AXES)
C This is for the sake of Dij tensor compuations.
         do j = 1, nslip
		    p(j,1) = 0.5d0*(slip_s(j,1)*slip_n(j,1) +
     1                      slip_s(j,1)*slip_n(j,1))
            p(j,2) = 0.5d0*(slip_s(j,2)*slip_n(j,2) + 
     1                      slip_s(j,2)*slip_n(j,2))							
			p(j,3) = 0.5d0*(slip_s(j,3)*slip_n(j,3) + 
     1                      slip_s(j,3)*slip_n(j,3))	
			p(j,4) = 0.5d0*(slip_s(j,1)*slip_n(j,2) + 
     1                      slip_s(j,2)*slip_n(j,1))			
			p(j,5) = 0.5d0*(slip_s(j,2)*slip_n(j,3) + 
     1                      slip_s(j,3)*slip_n(j,2))	 	 
			p(j,6) = 0.5d0*(slip_s(j,3)*slip_n(j,1) + 
     1                      slip_s(j,1)*slip_n(j,3))
         end do
C		 		 
C ANTI-SYMMETRIC SCHMID TENSOR (EXPRESSED W.R.T GLOBAL AXES)
C This is for the sake of spin tensor compuations, which are 
C also computed and expressed w.r.t. Initial (global) Axes
         do j = 1, nslip
            w12(j) = 0.5d0 * (slip_s(j,1)*slip_n(j,2) -					
     1						  slip_s(j,2)*slip_n(j,1))
			w23(j) = 0.5d0 * (slip_s(j,2)*slip_n(j,3) -					
     1						  slip_s(j,3)*slip_n(j,2))
			w31(j) = 0.5d0 * (slip_s(j,3)*slip_n(j,1) -					
     1						  slip_s(j,1)*slip_n(j,3))
         end do
C
C					   -------------------------------------------------
C					     #                    PijDij                  #
C		                        For the Computation of Tau_dot			
C					   -------------------------------------------------
		 do j = 1, nslip
			PijDijDev(j)=   p(j,1)*Dij_dev(1) +      p(j,2)*Dij_dev(2)
     1			           +p(j,3)*Dij_dev(3) + 2.d0*p(j,4)*Dij_dev(4)		
     2				  +2.d0*p(j,5)*Dij_dev(5) + 2.d0*p(j,6)*Dij_dev(6)            
		 end do
C
C					   -------------------------------------------------
C					    #  UPDATE Resolved Shear stress: TAU (alpha)  #
C					   -------------------------------------------------
C
C Thermal Multiplier
		 thermal_factor = (tempR/temp)**xi								
		 tauR(1:maxSS) = 0.d0
C          	
C Compute Tau Reference (crystal strength)	 		 
		 do j = 1, nslip
		    do kk = 1,nslip
			   alphaInt = abs( p(j,1)*p(kk,1) +    p(j,2)*p(kk,2)
     1			         +     p(j,3)*p(kk,3) + 2.*p(j,4)*p(kk,4)		
     2				     +  2.*p(j,5)*p(kk,5) + 2.*p(j,6)*p(kk,6) )/0.5d0
               tauR(j) = tauR(j) + alphaInt*g*b_v(kk)*sqrt(den_im(kk))
			enddo				
         enddo
		 tauR(1:nslip) = (tauR(1:nslip) +yield/EM_factor)*thermal_factor
C		 		 
C Define Integration Parameters to pass to ODEINT (initial value ODE problem)  
		 tauStart(1:nslip) = tau(1:nslip) 
	 	 NVAR = nslip	 
		 Time = totalTime - dt
		 eps  = 0.01d0
		 hmin = dt/24.d0 
C		 
C PACK Parameters into Array PAR to send to tauDot and tauResidual
C These are all required input to solve for the resolved shear stress
		 PAR (1)                      = dble(nslip)
		 PAR (        2 :    NVAR+1 ) = recip_m(1:nslip)
		 PAR (   NVAR+2 :  2*NVAR+1 ) = ref_gdot(1:nslip)
		 PAR ( 2*NVAR+2 :  3*NVAR+1 ) = PijDijDev(1:nslip)
		 PAR ( 3*NVAR+2 :  4*NVAR+1 ) = p(1:nslip,1)
		 PAR ( 4*NVAR+2 :  5*NVAR+1 ) = p(1:nslip,2)
		 PAR ( 5*NVAR+2 :  6*NVAR+1 ) = p(1:nslip,3)
		 PAR ( 6*NVAR+2 :  7*NVAR+1 ) = p(1:nslip,4)
		 PAR ( 7*NVAR+2 :  8*NVAR+1 ) = p(1:nslip,5)
		 PAR ( 8*NVAR+2 :  9*NVAR+1 ) = p(1:nslip,6)
		 PAR ( 9*NVAR+2 : 10*NVAR+1 ) = tau(1:nslip)
		 PAR (10*NVAR+2 : 11*NVAR+1 ) = tauR(1:nslip)
		 PAR (11*NVAR+2)              = twomu
		 PAR (size(PAR,1))            = 0.d0 !This takes the value of time step increment 'h'		 
C					 
C Integrate tauDot using an ODE solver subroutine       
         call odeint(tauDot, tauResidual, tauStart, NVAR,
     1               Time, Time + dt, eps, dt, hmin, PAR)
C Update Tau 
         tau(1:nslip) = tauStart(1:nslip)
C		 		
C					   -------------------------------------------------
C					   #  COMPUTE CRYSTALLOGRAPHIC SHEAR STRAIN RATE   # 
C					   -------------------------------------------------							    		 
         do j = 1, nslip
C Power-law 
			gdot(j)   =  ref_gdot(j)*(tau(j)/tauR(j))*                   
     1                  (abs(tau(j)/tauR(j)))**(recip_m(j)-1.)          
		 end do 
C					   -------------------------------------------------
C					   #  COMPUTE PLASTIC RATES OF DEFORMATION & SPIN  # 
C					   -------------------------------------------------				
C Initialize Values
		 D11_p    = 0.d0					
		 D22_p    = 0.d0
		 D33_p    = 0.d0
		 D12_p    = 0.d0
		 D23_p    = 0.d0
		 D31_p    = 0.d0
		 spin_p12 = 0.d0
		 spin_p23 = 0.d0
		 spin_p31 = 0.d0
	     spin_21  = Spin(1)
	     spin_32  = Spin(2)
	     spin_13  = Spin(3)
	     spin_12  =-Spin(1)
	     spin_23  =-Spin(2)
	     spin_31  =-Spin(3)		
C
         do j = 1 , nslip
C Plastic Deformation Rate,
			if (abs(gdot(j)).gt.abs(ref_gdot(j))) then  
				D11_p	 = D11_p	+ p(j,1)*gdot(j)						
				D22_p	 = D22_p	+ p(j,2)*gdot(j)
				D33_p    = D33_p	+ p(j,3)*gdot(j)
				D12_p	 = D12_p	+ p(j,4)*gdot(j)
				D23_p	 = D23_p	+ p(j,5)*gdot(j)
				D31_p	 = D31_p	+ p(j,6)*gdot(j)
				spin_p12 = spin_p12 + w12(j)*gdot(j)						
				spin_p23 = spin_p23 + w23(j)*gdot(j)						
				spin_p31 = spin_p31 + w31(j)*gdot(j)
			endif						
         end do
C
C					   -------------------------------------------------		 
C                       #                COMPUTE ELASTIC SPIN         #	 
C					   -------------------------------------------------
		 spin_p21 = -spin_p12				
		 spin_p13 = -spin_p31						 
		 spin_p32 = -spin_p23				
		 spin_e21 =  spin_21 - spin_p21
		 spin_e12 = -spin_e21
		 spin_e32 =  spin_32 - spin_p32
		 spin_e23 = -spin_e32
		 spin_e13 =  spin_13 - spin_p13
		 spin_e31 = -spin_e13
C
C
C					   -------------------------------------------------
C					    # UPDATE/NORMALIZE SLIP NORMALS AND DIRECTIONS #
C                            NOTE: STILL REFERENCED TO GLOBAL AXES
C					   -------------------------------------------------		 
		 do j = 1 , nslip
C
C n_dot and s_dot, 
		   an_dot1 = spin_e12*slip_n(j,2) + spin_e13*slip_n(j,3)	  
		   an_dot2 = spin_e21*slip_n(j,1) + spin_e23*slip_n(j,3)
		   an_dot3 = spin_e31*slip_n(j,1) + spin_e32*slip_n(j,2)
		   as_dot1 = spin_e12*slip_s(j,2) + spin_e13*slip_s(j,3)
		   as_dot2 = spin_e21*slip_s(j,1) + spin_e23*slip_s(j,3)
		   as_dot3 = spin_e31*slip_s(j,1) + spin_e32*slip_s(j,2)
C
C n and s updated
		   slip_n(j,1) = slip_n(j,1) + an_dot1*dt					  
		   slip_n(j,2) = slip_n(j,2) + an_dot2*dt
		   slip_n(j,3) = slip_n(j,3) + an_dot3*dt
		   slip_s(j,1) = slip_s(j,1) + as_dot1*dt
		   slip_s(j,2) = slip_s(j,2) + as_dot2*dt
		   slip_s(j,3) = slip_s(j,3) + as_dot3*dt	
		 end do
C         
C Normalize slip vectors to unit magnitude		 		 
		 call unit_vector( 3,nslip,slip_n,slip_s)		        		 
C
C			   ---------------------------------------------------------
C              # EXPRESS THE ROTATED STRESS AT BEGINNING OF INCREMENT  #
C              # AS CAUCHY STRESS:: S = R*S(corot)*R^T ,for R = FU^-1  #
C              # NOTE: RNEW is the Transpose of that in F=R*U          #
C			   ---------------------------------------------------------				
C Cauchy Stress
         Rnew = Transpose(Rnew)
         call ROTATE(stressNew(k,:),Rnew,stressOld(k,1),stressOld(k,2),
     1               stressOld(k,3),stressOld(k,4),stressOld(k,5),
     2               stressOld(k,6) )

         sigGLB1 = stressNew(k,1)
         sigGLB2 = stressNew(k,2)
         sigGLB3 = stressNew(k,3)
         sigGLB4 = stressNew(k,4)
         sigGLB5 = stressNew(k,5)
         sigGLB6 = stressNew(k,6)
         Rnew = Transpose(Rnew)
          
C COMPUTE THE DEVIATORIC PARTS
		press =  ( sigGLB1 + sigGLB2 + sigGLB3 )/3.d0
		sigGLB1 = sigGLB1  - press
		sigGLB2 = sigGLB2  - press
		sigGLB3 = sigGLB3  - press
C
C					   -------------------------------------------------
C				       #            UPDATE CAUCHY STRESSES             # 
C					   -------------------------------------------------
C COMPUTE THE NON-CONSTITUTIVE PART OF STRESS RATE: "SIGMA.W* - W*.SIGMA"		  
		  o11    = - 2.d0*( spin_e12*sigGLB4 + spin_e13*sigGLB6 )
		  o22    = - 2.d0*( spin_e21*sigGLB4 + spin_e23*sigGLB5 )	
		  o33    = - 2.d0*( spin_e31*sigGLB6 + spin_e32*sigGLB5 )				 	
		  o12    =   spin_e12*sigGLB1 - spin_e12*sigGLB2 -
     1               spin_e13*sigGLB5 - spin_e23*sigGLB6								
		  o23    =   spin_e23*sigGLB2 - spin_e23*sigGLB3 -
     1               spin_e21*sigGLB6 - spin_e31*sigGLB4					
		  o31    =   spin_e13*sigGLB1 - spin_e13*sigGLB3 -
     1               spin_e32*sigGLB4 - spin_e12*sigGLB5
C
C UPDATE THE CAUCHY STRESSES, BUT EXPRESSED W.R.T. GLOBAL AXES	
		  sigGLB1 = sigGLB1 + dt*o11 + twomu*dt*(Dij_dev(1)-D11_p)
		  sigGLB2 = sigGLB2 + dt*o22 + twomu*dt*(Dij_dev(2)-D22_p)	 
		  sigGLB3 = sigGLB3 + dt*o33 + twomu*dt*(Dij_dev(3)-D33_p)
		  sigGLB4 = sigGLB4 + dt*o12 + twomu*dt*(Dij_dev(4)-D12_p)
		  sigGLB5 = sigGLB5 + dt*o23 + twomu*dt*(Dij_dev(5)-D23_p)
		  sigGLB6 = sigGLB6 + dt*o31 + twomu*dt*(Dij_dev(6)-D31_p)		
C		 
C					   -------------------------------------------------
C						      # UPDATE TEMPERATURE AND PLASTIC WORK INC# 
C					   -------------------------------------------------
C (Deviatoric) StressPower,	
		 DijSij    =        sigGLB1*D11_p  +       sigGLB2*D22_p +    
     1                      sigGLB3*D33_p  +  2.d0*sigGLB4*D12_p +
     2                 2.d0*sigGLB5*D23_p  +  2.d0*sigGLB6*D31_p	 
C Update plastic work from stress power
         plasticWorkInc = dt*abs(DijSij) 								    
		 eta            = density(k)*Cp
C Adiabatic Temp Update, 
	     tempInc        = plasticWorkInc/eta* Chi           		    
C Update normal stress components with volumetric part		 
		 press  = press  + bulk*dt*(traceDij)  		  	
		 sigGLB1 = sigGLB1 + press
		 sigGLB2 = sigGLB2 + press
		 sigGLB3 = sigGLB3 + press


C					   -------------------------------------------------
C				       #  STORE UPDATED COROTATIONAL STRESS :: R^T*S*R  #
C					   -------------------------------------------------
         call ROTATE(stressNew(k,:),Rnew,sigGLB1,sigGLB2,sigGLB3,
     1                                   sigGLB4,sigGLB5,sigGLB6)
C
C					   -------------------------------------------------
C					       # UPDATE INTERNAL VARIABLES: RHO_M & RHO_IM # 
C					   -------------------------------------------------						 
C Define Integration Parameters to pass to ODEINT (initial value ODE problem)  
		 eps  = 0.05d0
		 hmin = dt/24.d0 
		 Time = totalTime - dt
		 NVAR = 2
C
C Reset PAR
		 PAR( 1: size(PAR,1) )  = 0.d0		 		 
C PACK Parameters into Array PAR to send to rhoDot and rhoResidual
		 PAR (       1        ) = dble(nslip)
		 PAR (       2        ) = g_source
		 PAR (       3        ) = g_immob
		 PAR (       4        ) = g_minter
		 PAR (       5        ) = g_recov
		 PAR (       6        ) = thermal_coef
		 PAR (      11   :   nslip+10 ) = gdot(1:nslip)
		 PAR (  nslip+11 : 2*nslip+10 ) = b_v(1:nslip)
		 PAR (   size(PAR,1)  ) = 0.d0							!This later takes on a value of 'dt' (actually 'h')
C
C Call rho_dot integrator 'nslip' times
		 do j = 1, nslip
             rhoStart(1) = den_im(j)							!Rho_im Starting Value 
             rhoStart(2) = den_m(j)								!Rho_ m Starting Value
		     PAR (     2*nslip+11      ) = dble(j)				!assigns jslip
 			 PAR (2*nslip+12:2*nslip+13) = rhoStart(1:2)
C
C Call Nested Subroutines for Adaptive RK5 Integration or Back Euler
             call odeint(rhoDot,rhoResidual,rhoStart, NVAR,  
     1	                 Time, Time + dt, eps, dt, hmin, PAR )  
             den_im(j) = rhoStart(1) 							!Rho_im Updated Value 
             den_m(j)  = rhoStart(2) 							!Rho_ m Updated Value
		 end do
C
C Reset PAR
		 PAR( 1: size(PAR,1) ) = 0.d0
C		 
C					   -------------------------------------------------
C					             # UPDATE ALL OTHER INTERNAL VARIABLES # 
C					   -------------------------------------------------						 
C Increment in shear slip:: Using an 'effective' measure; 
C dgamma = dt*sqrt([D_p]:[D_p])
		 dgamma = dt*sqrt(( D11_p**2.d0 +     D22_p**2.d0 +     D33_p**2.d0+
     1                 2.d0*D12_p**2.d0 +2.d0*D23_p**2.d0 +2.d0*D31_p**2.d0 ))
C
C Internal variables not dependent on slip system	
		 Psi21 = Psi21 + spin_e21*dt
		 Psi32 = Psi32 + spin_e32*dt
		 Psi13 = Psi13 + spin_e13*dt										
		 stateNew( k,1) = Psi21*(180.d0/pi)				!Angle 21 through which Slip system rotated,
		 stateNew( k,2) = Psi32*(180.d0/pi)				!Angle 32 through which Slip system rotated,
		 stateNew( k,3) = Psi13*(180.d0/pi)				!Angle 13 through which Slip system rotated,
         stateNew( k,4) = stateOld( k,4) + dgamma		!New Shear Slip,
		 stateNew( k,5) = sum(tauR(1:nslip))/nslip		!New Reference Tau
         stateNew( k,6) = stateOld( k,6) + tempInc		!New Temperature
C
C Internal variables dependent on slip system
		 do j = 1, nslip									!Loop over slip-system alpha:
 			 stateNew( k,8  + (j-1)*10 ) = tau(j) 			      !shear stress     
 			 stateNew( k,9  + (j-1)*10 ) = gdot(j)          	  !shear strain-rate
			 stateNew( k,10 + (j-1)*10 ) = den_im(j)          	  !rho immobile
			 stateNew( k,11 + (j-1)*10 ) = den_m(j)         	  !rho mobile
			 stateNew( k,12 + (j-1)*10 ) = slip_n(j,1)       	  !slip plane normal (x)
			 stateNew( k,13 + (j-1)*10 ) = slip_n(j,2)       	  !slip plane normal (y)
			 stateNew( k,14 + (j-1)*10 ) = slip_n(j,3)       	  !slip plane normal (z)
			 stateNew( k,15 + (j-1)*10 ) = slip_s(j,1)       	  !slip direction    (x)
			 stateNew( k,16 + (j-1)*10 ) = slip_s(j,2)       	  !slip direction    (y) 
			 stateNew( k,17 + (j-1)*10 ) = slip_s(j,3)       	  !slip direction    (z)
		 end do												!Endloop
C		 
C					   -------------------------------------------------
C					                             # UPDATE THE ENERGIES # 
C					   -------------------------------------------------				
C Update the dissipated inelastic specific energy -
	     enerInelasNew( k) = enerInelasOld( k) + 
     1	                     plasticWorkInc / density( k)		            
C Update the specific internal energy -
         stressPower = half * (
     1      ( stressOld( k,1)+stressNew( k,1) ) * (Dij(1) - D11_p)*dt +
     2      ( stressOld( k,2)+stressNew( k,2) ) * (Dij(2) - D22_p)*dt +
     3      ( stressOld( k,3)+stressNew( k,3) ) * (Dij(3) - D33_p)*dt ) +
     4      ( stressOld( k,4)+stressNew( k,4) ) * (Dij(4) - D12_p)*dt +     
     5      ( stressOld( k,5)+stressNew( k,5) ) * (Dij(5) - D23_p)*dt +    
     6      ( stressOld( k,6)+stressNew( k,6) ) * (Dij(6) - D31_p)*dt     
C
		 enerInternNew( k) = enerInternOld( k) 
     1                      + stressPower / density( k)	   
	  enddo	  
      return
	  end
C #=====================================================================
      include 'TensorSubs.inc'  
      include 'UpdateRhos.inc'
      include 'UpdateTaus.inc'
      include 'ODEINT.inc'
C ======================================================================