223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

whitespace not needed

This commit is contained in:
Martin Diehl 2019-01-07 07:24:02 +01:00
parent 705d55a3a5
commit 4037085f6c
4 changed files with 74 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
src/plastic_disloUCLA.f90
View File

@ -227,7 +227,7 @@ subroutine plastic_disloUCLA_init()
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), &
structure(1:3))
prm%rho0 = config%getFloats('rhoedge0', requiredShape=shape(prm%Nslip))
prm%rho0 = config%getFloats('rhoedge0', requiredShape=shape(prm%Nslip))
prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredShape=shape(prm%Nslip))
prm%v0 = config%getFloats('v0', requiredShape=shape(prm%Nslip))
prm%burgers = config%getFloats('slipburgers', requiredShape=shape(prm%Nslip))
@ -268,7 +268,7 @@ subroutine plastic_disloUCLA_init()
prm%clambda = math_expand(prm%clambda, prm%Nslip)
prm%atomicVolume = math_expand(prm%atomicVolume, prm%Nslip)
prm%minDipDistance = math_expand(prm%minDipDistance, prm%Nslip)
prm%tau0 = prm%tau_peierls + prm%SolidSolutionStrength
! sanity checks
@ -280,7 +280,7 @@ subroutine plastic_disloUCLA_init()
if (any(prm%tau_peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls'
if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' d0'
if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' qsd'
!if (plastic_disloUCLA_CAtomicVolume(instance) <= 0.0_pReal) &
! call IO_error(211_pInt,el=instance,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOUCLA_label//')')
else slipActive
@ -338,7 +338,7 @@ subroutine plastic_disloUCLA_init()
plasticState(p)%sizePostResults = sum(plastic_disloUCLA_sizePostResult(:,phase_plasticityInstance(p)))
allocate(prm%forestProjectionEdge(prm%totalNslip,prm%totalNslip),source = 0.0_pReal)
i = 0_pInt
mySlipFamilies: do f = 1_pInt,size(prm%Nslip,1)
index_myFamily = sum(prm%Nslip(1:f-1_pInt))
@ -416,7 +416,7 @@ subroutine plastic_disloUCLA_dependentState(instance,of)
dst%mfp(:,of) = prm%grainSize/(1.0_pReal+prm%grainSize*dst%dislocationSpacing(:,of)/prm%Clambda)
dst%dislocationSpacing(:,of) = dst%mfp(:,of) ! ToDo: Hack to recover wrong behavior for the moment
end associate
@ -450,7 +450,7 @@ pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,inst
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
call kinetics(Mp,Temperature,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
@ -461,7 +461,7 @@ pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,inst
+ dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo slipSystems
end associate
end subroutine plastic_disloUCLA_LpAndItsTangent
@ -600,12 +600,12 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe
end function plastic_disloUCLA_postResults
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the
! resolved stresss
!> @details Derivatives and resolved stress are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
!> @details Derivatives and resolved stress are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics(Mp,Temperature,instance,of, &
gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg,tau_slip_pos1,tau_slip_neg1)
@ -642,13 +642,13 @@ pure subroutine kinetics(Mp,Temperature,instance,of, &
integer(pInt) :: j
associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
do j = 1_pInt, prm%totalNslip
tau_slip_pos(j) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j))
tau_slip_neg(j) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j))
enddo
if (present(tau_slip_pos1)) tau_slip_pos1 = tau_slip_pos
if (present(tau_slip_neg1)) tau_slip_neg1 = tau_slip_neg
@ -735,7 +735,7 @@ pure subroutine kinetics(Mp,Temperature,instance,of, &
) &
/(2.0_pReal*prm%burgers**2.0_pReal*tau_slip_neg &
+ prm%omega * prm%B* effectiveLength**2.0_pReal* needsGoodName )**2.0_pReal
dgdot_dtauslip_neg = DotGamma0 * dvel_slip * 0.5_pReal
else where significantNegativeTau2
dgdot_dtauslip_neg = 0.0_pReal

54
src/plastic_isotropic.f90
View File

@ -44,7 +44,7 @@ module plastic_isotropic
aTolShear
integer(pInt) :: &
of_debug = 0_pInt
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
integer(kind(undefined_ID)), allocatable, dimension(:) :: &
outputID
logical :: &
dilatation
@ -119,7 +119,7 @@ subroutine plastic_isotropic_init()
p, i, &
NipcMyPhase, &
sizeState, sizeDotState
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
integer(kind(undefined_ID)) :: &
@ -154,7 +154,7 @@ subroutine plastic_isotropic_init()
dot => dotState(phase_plasticityInstance(p)), &
stt => state(phase_plasticityInstance(p)), &
config => config_phase(p))
#ifdef DEBUG
if (p==material_phase(debug_g,debug_i,debug_e)) then
prm%of_debug = phasememberAt(debug_g,debug_i,debug_e)
@ -175,22 +175,22 @@ subroutine plastic_isotropic_init()
prm%a = config%getFloat('a')
prm%aTolFlowStress = config%getFloat('atol_flowstress',defaultVal=1.0_pReal)
prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal)
prm%dilatation = config%keyExists('/dilatation/')
!--------------------------------------------------------------------------------------------------
! sanity checks
extmsg = ''
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//'aTolShear '
if (prm%tau0 < 0.0_pReal) extmsg = trim(extmsg)//'tau0 '
if (prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0 '
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//'n '
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//'aTolShear '
if (prm%tau0 < 0.0_pReal) extmsg = trim(extmsg)//'tau0 '
if (prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0 '
if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//'n '
if (prm%tausat <= prm%tau0) extmsg = trim(extmsg)//'tausat '
if (prm%a <= 0.0_pReal) extmsg = trim(extmsg)//'a '
if (prm%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//'m '
if (prm%a <= 0.0_pReal) extmsg = trim(extmsg)//'a '
if (prm%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//'m '
if (prm%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//'atol_flowstress '
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//'atol_shear '
if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//'atol_shear '
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') &
@ -242,7 +242,7 @@ subroutine plastic_isotropic_init()
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(2:2,1:NipcMyPhase)
plasticState(p)%accumulatedSlip => plasticState(p)%state (2:2,1:NipcMyPhase)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
@ -289,15 +289,15 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
k, l, m, n
associate(prm => param(instance), stt => state(instance))
Mp_dev = math_deviatoric33(Mp)
squarenorm_Mp_dev = math_mul33xx33(Mp_dev,Mp_dev)
norm_Mp_dev = sqrt(squarenorm_Mp_dev)
norm_Mp_dev = sqrt(squarenorm_Mp_dev)
if (norm_Mp_dev > 0.0_pReal) then
gamma_dot = prm%gdot0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%fTaylor*stt%flowstress(of))) **prm%n
Lp = Mp_dev/norm_Mp_dev * gamma_dot/prm%fTaylor
Lp = Mp_dev/norm_Mp_dev * gamma_dot/prm%fTaylor
#ifdef DEBUG
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
.and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
@ -318,7 +318,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
end if
end associate
end subroutine plastic_isotropic_LpAndItsTangent
@ -338,7 +338,7 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of)
Li !< inleastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLi_dTstar !< derivative of Li with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Tstar !< Mandel stress ToDo: Mi?
integer(pInt), intent(in) :: &
@ -355,10 +355,10 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of)
k, l, m, n
associate(prm => param(instance), stt => state(instance))
Tstar_sph = math_spherical33(Tstar)
squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph,Tstar_sph)
norm_Tstar_sph = sqrt(squarenorm_Tstar_sph)
norm_Tstar_sph = sqrt(squarenorm_Tstar_sph)
if (prm%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! no stress or J2 plastitiy --> Li and its derivative are zero
gamma_dot = prm%gdot0 * (sqrt(1.5_pReal) * norm_Tstar_sph /(prm%fTaylor*stt%flowstress(of))) **prm%n
@ -404,15 +404,15 @@ subroutine plastic_isotropic_dotState(Mp,instance,of)
norm_Mp !< norm of the (deviatoric) Mandel stress
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
if (prm%dilatation) then
norm_Mp = sqrt(math_mul33xx33(Mp,Mp))
else
norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp)))
endif
gamma_dot = prm%gdot0 * (sqrt(1.5_pReal) * norm_Mp /(prm%fTaylor*stt%flowstress(of))) **prm%n
if (abs(gamma_dot) > 1e-12_pReal) then
if (dEq0(prm%tausat_SinhFitA)) then
saturation = prm%tausat
@ -431,7 +431,7 @@ subroutine plastic_isotropic_dotState(Mp,instance,of)
dot%flowstress (of) = hardening * gamma_dot
dot%accumulatedShear(of) = gamma_dot
end associate
end subroutine plastic_isotropic_dotState
@ -461,13 +461,13 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults)
o,c
associate(prm => param(instance), stt => state(instance))
if (prm%dilatation) then
norm_Mp = sqrt(math_mul33xx33(Mp,Mp))
else
norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp)))
endif
c = 0_pInt
outputsLoop: do o = 1_pInt,size(prm%outputID)
@ -483,7 +483,7 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults)
end select
enddo outputsLoop
end associate
end function plastic_isotropic_postResults

36
src/plastic_kinematichardening.f90
View File

@ -40,7 +40,7 @@ module plastic_kinehardening
crss0, & !< initial critical shear stress for slip (input parameter, per family)
theta0, & !< initial hardening rate of forward stress for each slip
theta1, & !< asymptotic hardening rate of forward stress for each slip >
theta0_b, & !< initial hardening rate of back stress for each slip >
theta0_b, & !< initial hardening rate of back stress for each slip >
theta1_b, & !< asymptotic hardening rate of back stress for each slip >
tau1, &
tau1_b, &
@ -226,7 +226,7 @@ subroutine plastic_kinehardening_init
prm%theta1 = config%getFloats('theta1', requiredShape=shape(prm%Nslip))
prm%theta0_b = config%getFloats('theta0_b', requiredShape=shape(prm%Nslip))
prm%theta1_b = config%getFloats('theta1_b', requiredShape=shape(prm%Nslip))
prm%gdot0 = config%getFloat('gdot0')
prm%n_slip = config%getFloat('n_slip')
@ -266,7 +266,7 @@ subroutine plastic_kinehardening_init
outputID = undefined_ID
outputSize = prm%totalNslip
select case(outputs(i))
case ('resistance')
outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0_pInt)
case ('accumulatedshear')
@ -316,7 +316,7 @@ subroutine plastic_kinehardening_init
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1_pInt
endIndex = endIndex + prm%totalNslip
endIndex = endIndex + prm%totalNslip
stt%crss_back => plasticState(p)%state (startIndex:endIndex,:)
dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
@ -335,7 +335,7 @@ subroutine plastic_kinehardening_init
endIndex = endIndex + prm%totalNslip
stt%sense => plasticState(p)%state (startIndex :endIndex ,:)
dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1_pInt
endIndex = endIndex + prm%totalNslip
stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:)
@ -345,7 +345,7 @@ subroutine plastic_kinehardening_init
endIndex = endIndex + prm%totalNslip
stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
end associate
@ -391,7 +391,7 @@ pure subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
+ dgdot_dtau_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+ dgdot_dtau_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo
end associate
end subroutine plastic_kinehardening_LpAndItsTangent
@ -424,10 +424,10 @@ subroutine plastic_kinehardening_deltaState(Mp,instance,of)
sense
associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
sense = merge(state(instance)%sense(:,of), & ! keep existing...
sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined
sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined
dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction
#ifdef DEBUG
@ -450,7 +450,7 @@ subroutine plastic_kinehardening_deltaState(Mp,instance,of)
dlt%chi0 (:,of) = 0.0_pReal
dlt%gamma0(:,of) = 0.0_pReal
end where
end associate
end subroutine plastic_kinehardening_deltaState
@ -474,13 +474,13 @@ subroutine plastic_kinehardening_dotState(Mp,instance,of)
gdot_pos,gdot_neg
real(pReal) :: &
sumGamma
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
dot%accshear(:,of) = abs(gdot_pos+gdot_neg)
sumGamma = sum(stt%accshear(:,of))
sumGamma = sum(stt%accshear(:,of))
do i = 1_pInt, prm%totalNslip
dot%crss(i,of) = dot_product(prm%interaction_SlipSlip(i,:),dot%accshear(:,of)) &
* ( prm%theta1(i) &
@ -493,8 +493,8 @@ subroutine plastic_kinehardening_dotState(Mp,instance,of)
(prm%theta0_b - prm%theta1_b &
+ prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))&
) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) &
)
)
end associate
end subroutine plastic_kinehardening_dotState
@ -551,9 +551,9 @@ function plastic_kinehardening_postResults(Mp,instance,of) result(postResults)
enddo
end select
c = c + prm%totalNslip
enddo outputsLoop
end associate
@ -595,7 +595,7 @@ pure subroutine kinetics(Mp,instance,of, &
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0_pInt
do i = 1_pInt, prm%totalNslip

28
src/plastic_phenopowerlaw.f90
View File

@ -221,7 +221,7 @@ subroutine plastic_phenopowerlaw_init
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), &
defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))])
prm%gdot0_slip = config%getFloat('gdot0_slip')
prm%n_slip = config%getFloat('n_slip')
prm%a_slip = config%getFloat('a_slip')
@ -392,7 +392,7 @@ end subroutine plastic_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @details asummes that deformation by dislocation glide affects twinned and untwinned volume
!> @details asummes that deformation by dislocation glide affects twinned and untwinned volume
! equally (Taylor assumption). Twinning happens only in untwinned volume
!--------------------------------------------------------------------------------------------------
pure subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
@ -419,9 +419,9 @@ pure subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance))
call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
slipSystems: do i = 1_pInt, prm%totalNslip
Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i)
@ -438,7 +438,7 @@ pure subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i)
enddo twinSystems
end associate
end subroutine plastic_phenopowerlaw_LpAndItsTangent
@ -578,8 +578,8 @@ end function plastic_phenopowerlaw_postResults
!--------------------------------------------------------------------------------------------------
!> @brief Shear rates on slip systems and their derivatives with respect to resolved stress
!> @details Derivatives are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
!> @details Derivatives are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(Mp,instance,of, &
@ -610,7 +610,7 @@ pure subroutine kinetics_slip(Mp,instance,of, &
logical :: nonSchmidActive
associate(prm => param(instance), stt => state(instance))
nonSchmidActive = size(prm%nonSchmidCoeff) > 0_pInt
do i = 1_pInt, prm%totalNslip
@ -656,7 +656,7 @@ end subroutine kinetics_slip
!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress.
! twinning is assumed to take place only in untwinned volume.
!> @details Derivates are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(Mp,instance,of,&
@ -672,7 +672,7 @@ pure subroutine kinetics_twin(Mp,instance,of,&
integer(pInt), intent(in) :: &
instance, &
of
real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: &
gdot_twin
real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: &
@ -681,17 +681,17 @@ pure subroutine kinetics_twin(Mp,instance,of,&
real(pReal), dimension(param(instance)%totalNtwin) :: &
tau_twin
integer(pInt) :: i
associate(prm => param(instance), stt => state(instance))
do i = 1_pInt, prm%totalNtwin
tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
enddo
where(tau_twin > 0.0_pReal)
gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction
* prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin
else where
else where
gdot_twin = 0.0_pReal
end where
@ -702,7 +702,7 @@ pure subroutine kinetics_twin(Mp,instance,of,&
dgdot_dtau_twin = 0.0_pReal
end where
endif
end associate
end subroutine kinetics_twin