# now $ID aware

# LF as line separation and some tidy up

code/constitutive_dislotwin.f90

@@ -1,4 +1,4 @@
-!* $Id: constitutive_dislotwin.f90 412 2009-09-18 15:37:14Z MPIE\c.kords $
+!* $Id$
 !************************************
 !*      Module: CONSTITUTIVE        *
 !************************************
@@ -122,7 +122,7 @@ character(len=1024) line
 
 !write(6,*)
 !write(6,'(a20,a20,a12)') '<<<+-  constitutive_',constitutive_dislotwin_label,' init  -+>>>'
-!write(6,*) '$Id: constitutive_dislotwin.f90 412 2009-09-18 15:37:14Z MPIE\c.kords $'
+!write(6,*) '$Id$'
 !write(6,*)
 
 maxNinstance = count(phase_constitution == constitutive_dislotwin_label)
@@ -160,7 +160,6 @@ constitutive_dislotwin_slipSystemLattice = 0.0_pReal
 constitutive_dislotwin_twinSystemLattice = 0.0_pReal
 constitutive_dislotwin_totalNslip        = 0_pInt
 constitutive_dislotwin_totalNtwin        = 0_pInt
-
 allocate(constitutive_dislotwin_CoverA(maxNinstance))
 allocate(constitutive_dislotwin_C11(maxNinstance))
 allocate(constitutive_dislotwin_C12(maxNinstance))
@@ -203,7 +202,6 @@ constitutive_dislotwin_Cthresholdtwin      = 0.0_pReal
 constitutive_dislotwin_relevantRho         = 0.0_pReal
 constitutive_dislotwin_Cslip_66            = 0.0_pReal
 constitutive_dislotwin_Cslip_3333          = 0.0_pReal
-
 allocate(constitutive_dislotwin_rhoEdge0(lattice_maxNslipFamily,maxNinstance))
 allocate(constitutive_dislotwin_rhoEdgeDip0(lattice_maxNslipFamily,maxNinstance)) 
 allocate(constitutive_dislotwin_burgersPerSlipFamily(lattice_maxNslipFamily,maxNinstance))
@@ -222,7 +220,6 @@ constitutive_dislotwin_v0PerSlipFamily          = 0.0_pReal
 constitutive_dislotwin_Ndot0PerTwinFamily       = 0.0_pReal
 constitutive_dislotwin_twinsizePerTwinFamily    = 0.0_pReal
 constitutive_dislotwin_CLambdaSlipPerSlipFamily = 0.0_pReal
-
 allocate(constitutive_dislotwin_interactionSlipSlip(lattice_maxNinteraction,maxNinstance)) 
 allocate(constitutive_dislotwin_interactionSlipTwin(lattice_maxNinteraction,maxNinstance)) 
 allocate(constitutive_dislotwin_interactionTwinSlip(lattice_maxNinteraction,maxNinstance)) 
@@ -232,6 +229,8 @@ constitutive_dislotwin_interactionSlipTwin = 0.0_pReal
 constitutive_dislotwin_interactionTwinSlip = 0.0_pReal
 constitutive_dislotwin_interactionTwinTwin = 0.0_pReal
 
+
+
 !* Readout data from material.config file
 rewind(file)
 line = ''
@@ -272,27 +271,38 @@ do                                                       ! read thru sections of
        case ('c44')
               constitutive_dislotwin_C44(i) = IO_floatValue(line,positions,2)
        case ('nslip')
-              forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_Nslip(j,i) = IO_intValue(line,positions,1+j)
+              forall (j = 1:lattice_maxNslipFamily) &
+                constitutive_dislotwin_Nslip(j,i) = IO_intValue(line,positions,1+j)
        case ('ntwin')
-              forall (j = 1:lattice_maxNtwinFamily) constitutive_dislotwin_Ntwin(j,i) = IO_intValue(line,positions,1+j)
+              forall (j = 1:lattice_maxNtwinFamily) &
+                constitutive_dislotwin_Ntwin(j,i) = IO_intValue(line,positions,1+j)
        case ('rhoedge0')
-              forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_rhoEdge0(j,i) = IO_floatValue(line,positions,1+j)
+              forall (j = 1:lattice_maxNslipFamily) &
+                constitutive_dislotwin_rhoEdge0(j,i) = IO_floatValue(line,positions,1+j)
        case ('rhoedgedip0')
-              forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_rhoEdgeDip0(j,i) = IO_floatValue(line,positions,1+j)
+              forall (j = 1:lattice_maxNslipFamily) &
+                constitutive_dislotwin_rhoEdgeDip0(j,i) = IO_floatValue(line,positions,1+j)
        case ('slipburgers')
-              forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_burgersPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
+              forall (j = 1:lattice_maxNslipFamily) &
+                constitutive_dislotwin_burgersPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
        case ('twinburgers')
-              forall (j = 1:lattice_maxNtwinFamily) constitutive_dislotwin_burgersPerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
+              forall (j = 1:lattice_maxNtwinFamily) &
+                constitutive_dislotwin_burgersPerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
        case ('qedge')
-              forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_QedgePerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
+              forall (j = 1:lattice_maxNslipFamily) &
+                constitutive_dislotwin_QedgePerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
        case ('v0')
-              forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_v0PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
+              forall (j = 1:lattice_maxNslipFamily) &
+                constitutive_dislotwin_v0PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
        case ('ndot0')
-              forall (j = 1:lattice_maxNtwinFamily) constitutive_dislotwin_Ndot0PerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
+              forall (j = 1:lattice_maxNtwinFamily) &
+                constitutive_dislotwin_Ndot0PerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
        case ('twinsize')
-              forall (j = 1:lattice_maxNtwinFamily) constitutive_dislotwin_twinsizePerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
+              forall (j = 1:lattice_maxNtwinFamily) &
+                constitutive_dislotwin_twinsizePerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
        case ('clambdaslip')
-              forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_CLambdaSlipPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
+              forall (j = 1:lattice_maxNslipFamily) &
+                constitutive_dislotwin_CLambdaSlipPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
        case ('grainsize')
               constitutive_dislotwin_GrainSize(i) = IO_floatValue(line,positions,2)
        case ('maxtwinfraction')
@@ -375,6 +385,7 @@ enddo
 maxTotalNslip = maxval(constitutive_dislotwin_totalNslip)
 maxTotalNtwin = maxval(constitutive_dislotwin_totalNtwin)      
 
+
 allocate(constitutive_dislotwin_burgersPerSlipSystem(maxTotalNslip, maxNinstance))
 allocate(constitutive_dislotwin_burgersPerTwinSystem(maxTotalNtwin, maxNinstance))
 allocate(constitutive_dislotwin_QedgePerSlipSystem(maxTotalNslip, maxNinstance))
@@ -441,6 +452,7 @@ do i = 1,maxNinstance
       select case(constitutive_dislotwin_output(o,i))
         case('edge_density', &
              'dipole_density', &
+
              'shear_rate_slip', &
              'mfp_slip', &
              'resolved_stress_slip', &
@@ -532,6 +544,7 @@ do i = 1,maxNinstance
          constitutive_dislotwin_interactionSlipSlip(lattice_interactionSlipSlip(constitutive_dislotwin_slipSystemLattice(s1,i), &
                                                                                 constitutive_dislotwin_slipSystemLattice(s2,i), &
                                                                                 myStructure),i)
+
    enddo; enddo
    
    do s1 = 1,constitutive_dislotwin_totalNslip(i)
@@ -577,6 +590,7 @@ function constitutive_dislotwin_stateInit(myInstance)
 !* initial microstructural state                                     *
 !*********************************************************************
 use prec,    only: pReal,pInt
+
 use math,    only: pi
 use lattice, only: lattice_maxNslipFamily,lattice_maxNtwinFamily
 implicit none
@@ -586,7 +600,12 @@ integer(pInt) :: myInstance
 real(pReal), dimension(constitutive_dislotwin_sizeState(myInstance))  :: constitutive_dislotwin_stateInit
 !* Local variables
 integer(pInt) s0,s1,s,t,f,ns,nt
-real(pReal), dimension(constitutive_dislotwin_totalNslip(myInstance)) :: rhoEdge0,rhoEdgeDip0,invLambdaSlip0,MeanFreePathSlip0,tauSlipThreshold0              
+real(pReal), dimension(constitutive_dislotwin_totalNslip(myInstance)) :: rhoEdge0, &
+                                                                         rhoEdgeDip0, &
+                                                                         invLambdaSlip0, &
+                                                                         MeanFreePathSlip0, &
+                                                                         tauSlipThreshold0
+
 real(pReal), dimension(constitutive_dislotwin_totalNtwin(myInstance)) :: MeanFreePathTwin0,TwinVolume0
 
 ns = constitutive_dislotwin_totalNslip(myInstance)
@@ -594,6 +613,7 @@ nt = constitutive_dislotwin_totalNtwin(myInstance)
 constitutive_dislotwin_stateInit = 0.0_pReal
 
 !* Initialize basic slip state variables
+
 s1 = 0_pInt
 do f = 1,lattice_maxNslipFamily
    s0 = s1 + 1_pInt
@@ -603,6 +623,7 @@ do f = 1,lattice_maxNslipFamily
       rhoEdgeDip0(s) = constitutive_dislotwin_rhoEdgeDip0(f,myInstance)
    enddo 
 enddo
+
 constitutive_dislotwin_stateInit(1:ns)      = rhoEdge0
 constitutive_dislotwin_stateInit(ns+1:2*ns) = rhoEdgeDip0
 
@@ -615,6 +636,7 @@ constitutive_dislotwin_stateInit(2*ns+nt+1:3*ns+nt) = invLambdaSlip0
 forall (s = 1:ns) &
 MeanFreePathSlip0(s) = &
 constitutive_dislotwin_GrainSize(myInstance)/(1.0_pReal+invLambdaSlip0(s)*constitutive_dislotwin_GrainSize(myInstance)) 
+
 constitutive_dislotwin_stateInit(4*ns+2*nt+1:5*ns+2*nt) = MeanFreePathSlip0
 
 forall (s = 1:ns) &
@@ -623,6 +645,7 @@ constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSys
 sqrt(dot_product(rhoEdge0,constitutive_dislotwin_interactionMatrixSlipSlip(1:ns,s,myInstance)))
 constitutive_dislotwin_stateInit(5*ns+3*nt+1:6*ns+3*nt) = tauSlipThreshold0
 
+
 !* Initialize dependent twin microstructural variables
 forall (t = 1:nt) &
 MeanFreePathTwin0(t) = constitutive_dislotwin_GrainSize(myInstance)
@@ -633,6 +656,7 @@ TwinVolume0(t) = &
 (pi/6.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance)*MeanFreePathTwin0(t)**(2.0_pReal)
 constitutive_dislotwin_stateInit(6*ns+4*nt+1:6*ns+5*nt) = TwinVolume0
 
+
 !write(6,*) '#STATEINIT#'
 !write(6,*)
 !write(6,'(a,/,4(3(f30.20,x)/))') 'RhoEdge',rhoEdge0
@@ -691,12 +715,14 @@ ns = constitutive_dislotwin_totalNslip(myInstance)
 nt = constitutive_dislotwin_totalNtwin(myInstance)
 
 !* Total twin volume fraction
+
 sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0
 
 !* Homogenized elasticity matrix
 constitutive_dislotwin_homogenizedC = (1.0_pReal-sumf)*constitutive_dislotwin_Cslip_66(:,:,myInstance)
 do i=1,nt
    constitutive_dislotwin_homogenizedC = &
+
    constitutive_dislotwin_homogenizedC + state(g,ip,el)%p(2*ns+i)*constitutive_dislotwin_Ctwin_66(:,:,i,myInstance)
 enddo 
 
@@ -736,6 +762,7 @@ myStructure = constitutive_dislotwin_structure(myInstance)
 ns = constitutive_dislotwin_totalNslip(myInstance)
 nt = constitutive_dislotwin_totalNtwin(myInstance)
 !* State: 1           :  ns         rho_edge
+
 !* State: ns+1        :  2*ns       rho_dipole
 !* State: 2*ns+1      :  2*ns+nt    f
 !* State: 2*ns+nt+1   :  3*ns+nt    1/lambda_slip
@@ -747,10 +774,14 @@ nt = constitutive_dislotwin_totalNtwin(myInstance)
 !* State: 6*ns+3*nt+1 :  6*ns+4*nt  threshold_stress_twin
 !* State: 6*ns+4*nt+1 :  6*ns+5*nt  twin volume
 
+
 !* Total twin volume fraction
+
 sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0
 
+
 !* Stacking fault energy
+
 sfe = 0.0002_pReal*Temperature-0.0396_pReal
 
 !* rescaled twin volume fraction for topology
@@ -761,19 +792,22 @@ state(g,ip,el)%p(2*ns+t)/constitutive_dislotwin_twinsizePerTwinSystem(t,myInstan
 !* 1/mean free distance between 2 forest dislocations seen by a moving dislocation
 forall (s = 1:ns) &
   state(g,ip,el)%p(2*ns+nt+s) = &
-sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1:2*ns)),constitutive_dislotwin_forestProjectionEdge(1:ns,s,myInstance)))/ &
+    sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1:2*ns)),&
+                     constitutive_dislotwin_forestProjectionEdge(1:ns,s,myInstance)))/ &
     constitutive_dislotwin_CLambdaSlipPerSlipSystem(s,myInstance) 
 
 !* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
 !$OMP CRITICAL (evilmatmul)
 state(g,ip,el)%p((3*ns+nt+1):(4*ns+nt)) = 0.0_pReal
-if (nt > 0_pInt) state(g,ip,el)%p((3*ns+nt+1):(4*ns+nt)) = &
+if (nt > 0_pInt) &
+  state(g,ip,el)%p((3*ns+nt+1):(4*ns+nt)) = &
     matmul(constitutive_dislotwin_interactionMatrixSlipTwin(1:ns,1:nt,myInstance),fOverStacksize(1:nt))/(1.0_pReal-sumf)
 !$OMP END CRITICAL (evilmatmul)
 
 !* 1/mean free distance between 2 twin stacks from different systems seen by a growing twin
 !$OMP CRITICAL (evilmatmul)
-if (nt > 0_pInt) state(g,ip,el)%p((4*ns+nt+1):(4*ns+2*nt)) = &
+if (nt > 0_pInt) &
+  state(g,ip,el)%p((4*ns+nt+1):(4*ns+2*nt)) = &
     matmul(constitutive_dislotwin_interactionMatrixTwinTwin(1:nt,1:nt,myInstance),fOverStacksize(1:nt))/(1.0_pReal-sumf)
 !$OMP END CRITICAL (evilmatmul)
 
@@ -802,6 +836,7 @@ state(g,ip,el)%p(5*ns+3*nt+s) = &
     constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(s,myInstance)*&
     sqrt(dot_product(state(g,ip,el)%p(1:ns),constitutive_dislotwin_interactionMatrixSlipSlip(1:ns,s,myInstance)))
 
+
 !* threshold stress for growing twin
 forall (t = 1:nt) &
   state(g,ip,el)%p(6*ns+3*nt+t) = &
@@ -815,6 +850,7 @@ forall (t = 1:nt) &
   state(g,ip,el)%p(6*ns+4*nt+t) = &
     (pi/6.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance)*state(g,ip,el)%p(5*ns+2*nt+t)**(2.0_pReal)
 
+
 !if ((ip==1).and.(el==1)) then
 !   write(6,*) '#MICROSTRUCTURE#'
 ! write(6,*)
@@ -823,6 +859,7 @@ state(g,ip,el)%p(6*ns+4*nt+t) = &
 ! write(6,'(a,/,4(3(f10.4,x)/))') 'Fraction',state(g,ip,el)%p(2*ns+1:2*ns+nt)
 !endif
 
+
 return
 end subroutine
 
@@ -872,6 +909,7 @@ ns = constitutive_dislotwin_totalNslip(myInstance)
 nt = constitutive_dislotwin_totalNtwin(myInstance)
 
 !* Total twin volume fraction
+
 sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0
 
 Lp = 0.0_pReal
@@ -888,9 +926,9 @@ do f = 1,lattice_maxNslipFamily                                 ! loop over all
       j = j+1_pInt
 
       !* Calculation of Lp
-
       !* Resolved shear stress on slip system
       tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) 
+
       !* Stress ratios
       StressRatio_p = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance)
       StressRatio_pminus1 = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal)
@@ -898,7 +936,8 @@ do f = 1,lattice_maxNslipFamily                                 ! loop over all
       BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature)
       !* Initial shear rates
       DotGamma0 = &
-	  state(g,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*constitutive_dislotwin_v0PerSlipSystem(f,myInstance)	        
+        state(g,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*&
+        constitutive_dislotwin_v0PerSlipSystem(f,myInstance)
 
       !* Shear rates due to slip
       gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau_slip(j))
@@ -931,8 +970,8 @@ do f = 1,lattice_maxNtwinFamily                                 ! loop over all
       j = j+1_pInt
 
       !* Calculation of Lp
-
       !* Resolved shear stress on twin system
+
       tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure))        
       !* Stress ratios
       StressRatio_r = (state(g,ip,el)%p(6*ns+3*nt+j)/tau_twin(j))**constitutive_dislotwin_r(myInstance)
@@ -942,10 +981,10 @@ do f = 1,lattice_maxNtwinFamily                                 ! loop over all
         gdot_twin(j) = &
           (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*lattice_shearTwin(index_myFamily+i,myStructure)*&
           state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r) 
-       
         dgdot_dtautwin(j) = ((gdot_twin(j)*constitutive_dislotwin_r(myInstance))/tau_twin(j))*StressRatio_r
       endif
 
+
       !* Plastic velocity gradient for mechanical twinning	   					   
       Lp = Lp + gdot_twin(j)*lattice_Stwin(:,:,index_myFamily+i,myStructure)
 
@@ -960,6 +999,7 @@ enddo
 
 dLp_dTstar = math_Plain3333to99(dLp_dTstar3333)
 
+
 !if ((ip==1).and.(el==1)) then
 !   write(6,*) '#MICROSTRUCTURE#'
 ! write(6,*)
@@ -967,6 +1007,7 @@ dLp_dTstar = math_Plain3333to99(dLp_dTstar3333)
 ! write(6,'(a,/,4(3(f10.4,x)/))') 'rhoEdgeDip',state(g,ip,el)%p(13:24)
 !endif
 
+
 return
 end subroutine
 
@@ -985,6 +1026,7 @@ function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,g,ip,el)
 !*  - constitutive_dotState : evolution of state variable            *
 !*********************************************************************
 use prec,     only: pReal,pInt,p_vec
+
 use math,     only: pi
 use mesh,     only: mesh_NcpElems,mesh_maxNips
 use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
@@ -996,14 +1038,17 @@ implicit none
 integer(pInt), intent(in) :: g,ip,el
 real(pReal), intent(in) :: Temperature
 real(pReal), dimension(6), intent(in) :: Tstar_v
+
 type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
 real(pReal), dimension(constitutive_dislotwin_sizeDotState(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
 constitutive_dislotwin_dotState
 !* Local variables
 integer(pInt) MyInstance,MyStructure,ns,nt,f,i,j,k,index_myFamily
-real(pReal) sumf,StressRatio_p,StressRatio_pminus1,BoltzmannRatio,DotGamma0,EdgeDipMinDistance,AtomicVolume,VacancyDiffusion,StressRatio_r
+real(pReal) sumf,StressRatio_p,StressRatio_pminus1,BoltzmannRatio,DotGamma0,&
+            EdgeDipMinDistance,AtomicVolume,VacancyDiffusion,StressRatio_r
 real(pReal), dimension(constitutive_dislotwin_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
 gdot_slip,tau_slip,DotRhoMultiplication,EdgeDipDistance,DotRhoEdgeEdgeAnnihilation,DotRhoEdgeDipAnnihilation,&
+
 ClimbVelocity,DotRhoEdgeDipClimb,DotRhoDipFormation
 real(pReal), dimension(constitutive_dislotwin_totalNtwin(phase_constitutionInstance(material_phase(g,ip,el)))) :: gdot_twin,tau_twin
    
@@ -1014,6 +1059,7 @@ ns = constitutive_dislotwin_totalNslip(myInstance)
 nt = constitutive_dislotwin_totalNtwin(myInstance)
 
 !* Total twin volume fraction
+
 sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0
 
 constitutive_dislotwin_dotState = 0.0_pReal
@@ -1026,6 +1072,7 @@ do f = 1,lattice_maxNslipFamily                                 ! loop over all
    do i = 1,constitutive_dislotwin_Nslip(f,myInstance)          ! process each (active) slip system in family
       j = j+1_pInt
 
+
       !* Resolved shear stress on slip system
       tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) 
       !* Stress ratios
@@ -1035,15 +1082,19 @@ do f = 1,lattice_maxNslipFamily                                 ! loop over all
       BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature)
       !* Initial shear rates
       DotGamma0 = &
-	  state(g,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*constitutive_dislotwin_v0PerSlipSystem(f,myInstance)	        
+        state(g,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*&
+        constitutive_dislotwin_v0PerSlipSystem(f,myInstance)            
        
       !* Shear rates due to slip
       gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau_slip(j))
 
+
       !* Multiplication
-      DotRhoMultiplication(j) = abs(gdot_slip(j))/(constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*state(g,ip,el)%p(4*ns+2*nt+j))
+      DotRhoMultiplication(j) = abs(gdot_slip(j))/&
+                                (constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*state(g,ip,el)%p(4*ns+2*nt+j))
 
       !* Dipole formation
+
       if (tau_slip(j) == 0.0_pReal) then
          DotRhoDipFormation(j) = 0.0_pReal
       else 
@@ -1070,6 +1121,7 @@ do f = 1,lattice_maxNslipFamily                                 ! loop over all
       !* Dislocation dipole climb
       AtomicVolume = &
         constitutive_dislotwin_CAtomicVolume(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)**(3.0_pReal)
+
       VacancyDiffusion = &
         constitutive_dislotwin_D0(myInstance)*exp(-constitutive_dislotwin_Qsd(myInstance)/(kB*Temperature))
       ClimbVelocity(j) = &
@@ -1082,6 +1134,7 @@ do f = 1,lattice_maxNslipFamily                                 ! loop over all
       constitutive_dislotwin_dotState(j) = &
         DotRhoMultiplication(j)-DotRhoDipFormation(j)-DotRhoEdgeEdgeAnnihilation(j)
 
+
       !* Edge dislocation dipole density rate of change
       constitutive_dislotwin_dotState(ns+j) = &
         DotRhoDipFormation(j)-DotRhoEdgeDipAnnihilation(j)-DotRhoEdgeDipClimb(j)
@@ -1108,9 +1161,11 @@ do f = 1,lattice_maxNtwinFamily                                 ! loop over all
           (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*&
           state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r) 
       endif
+
    enddo
 enddo
 
+
 !if ((ip==1).and.(el==1)) then
 !    write(6,*) '#DOTSTATE#'
 !   write(6,*)
@@ -1125,6 +1180,8 @@ enddo
 !   write(6,'(a,/,4(3(f30.20,x)/))') 'DipClimb',DotRhoEdgeDipClimb 
 !endif
 
+
+
 return
 end function
 
@@ -1206,11 +1263,9 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
      case ('edge_density')
        constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns)
        c = c + ns
-
      case ('dipole_density')
        constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p(ns+1:2*ns)
        c = c + ns
-
      case ('shear_rate_slip')
        j = 0_pInt
        do f = 1,lattice_maxNslipFamily                                 ! loop over all slip families
@@ -1227,18 +1282,17 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
              BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature)
              !* Initial shear rates
              DotGamma0 = &
-	         state(g,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*constitutive_dislotwin_v0PerSlipSystem(f,myInstance)	        
+               state(g,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)* &
+               constitutive_dislotwin_v0PerSlipSystem(f,myInstance)
        
              !* Shear rates due to slip
              constitutive_dislotwin_postResults(c+j) = &
                DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau)
        enddo ; enddo
        c = c + ns
-
      case ('mfp_slip')
        constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p((4*ns+2*nt+1):(5*ns+2*nt))
        c = c + ns
-
      case ('resolved_stress_slip')
        j = 0_pInt
        do f = 1,lattice_maxNslipFamily                                 ! loop over all slip families
@@ -1248,15 +1302,12 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
              constitutive_dislotwin_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
        enddo; enddo
        c = c + ns
-   
      case ('threshold_stress_slip')
        constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p((5*ns+3*nt+1):(6*ns+3*nt))
        c = c + ns
-
      case ('twin_fraction')
        constitutive_dislotwin_postResults(c+1:c+nt) = state(g,ip,el)%p((2*ns+1):(2*ns+nt))
        c = c + nt
-
      case ('shear_rate_twin')
        if (nt > 0_pInt) then 
          j = 0_pInt
@@ -1276,14 +1327,13 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
                  (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*&
                  state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r)
              endif
+
          enddo ; enddo
        endif
        c = c + nt
-
      case ('mfp_twin')
        constitutive_dislotwin_postResults(c+1:c+nt) = state(g,ip,el)%p((5*ns+2*nt+1):(5*ns+3*nt))
        c = c + nt
-
      case ('resolved_stress_twin')
        if (nt > 0_pInt) then
          j = 0_pInt
@@ -1295,11 +1345,9 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
          enddo; enddo
        endif
        c = c + nt
-
      case ('threshold_stress_twin')
        constitutive_dislotwin_postResults(c+1:c+nt) = state(g,ip,el)%p((6*ns+3*nt+1):(6*ns+4*nt))
        c = c + nt
-
    end select
 enddo