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

included newest hypela2 from marc 2005r3 

the whole thing should compile now, however some things are still missing
This commit is contained in:
Franz Roters 2007-03-28 08:12:48 +00:00
parent 9241c7de91
commit 08c5390f95
2 changed files with 140 additions and 149 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
trunk/CPFEM.f90
View File

@ -105,18 +105,18 @@
CPFEM_jacobi_all = 0.0_pReal CPFEM_jacobi_all = 0.0_pReal
! !
! *** User defined results !!! MISSING incorporate consti_Nresults *** ! *** User defined results !!! MISSING incorporate consti_Nresults ***
allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,texture_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
CPFEM_results = 0.0_pReal CPFEM_results = 0.0_pReal
! !
! *** Second Piola-Kirchoff stress tensor at (t=t0) and (t=t1) *** ! *** Second Piola-Kirchoff stress tensor at (t=t0) and (t=t1) ***
allocate(CPFEM_sigma_old(6,texture_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(CPFEM_sigma_old(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(CPFEM_sigma_new(6,texture_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(CPFEM_sigma_new(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
CPFEM_sigma_old = 0.0_pReal CPFEM_sigma_old = 0.0_pReal
CPFEM_sigma_new = 0.0_pReal CPFEM_sigma_new = 0.0_pReal
! !
! *** Plastic deformation gradient at (t=t0) and (t=t1) *** ! *** Plastic deformation gradient at (t=t0) and (t=t1) ***
allocate(CPFEM_Fp_old(3,3,texture_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(CPFEM_Fp_new(3,3,texture_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
CPFEM_Fp_old = 0.0_pReal CPFEM_Fp_old = 0.0_pReal
CPFEM_Fp_new = 0.0_pReal CPFEM_Fp_new = 0.0_pReal
! !
@ -187,7 +187,7 @@
! !
! *** Initialization of the matrices for t=t0 *** ! *** Initialization of the matrices for t=t0 ***
! data from constitutive? ! data from constitutive?
vf = constitutive_volfrac(iori,CPFEM_in,cp_en) !ÄÄÄ vf = constitutive_matVolFrac(iori,CPFEM_in,cp_en)*constitutive_texVolFrac(iori,CPFEM_in,cp_en) !ÄÄÄ
! *** Calculation of the solution at t=t1 *** ! *** Calculation of the solution at t=t1 ***
! QUESTION use the mod() as flag parameter in the call ?? ! QUESTION use the mod() as flag parameter in the call ??
@ -282,10 +282,10 @@
real(pReal) cs(6), cd(6,6), CPFEM_dt real(pReal) cs(6), cd(6,6), CPFEM_dt
integer(pInt) cp_en ,CPFEM_in, iori, ising, icut, iconv, isjaco integer(pInt) cp_en ,CPFEM_in, iori, ising, icut, iconv, isjaco
! *** Local variables *** ! *** Local variables ***
real(pReal) Fp_old(3,3), Fp_new(3,3), state_old(constitutive_Nstatevars) real(pReal) Fp_old(3,3), Fp_new(3,3), state_old(constitutive_Nstatevars(iori, CPFEM_in, cp_en))
real(pReal) state_new(constitutive_Nstatevars), Tstar_v(6), CPFEM_ffn(3,3), CPFEM_ffn1(3,3) real(pReal) state_new(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), Tstar_v(6), CPFEM_ffn(3,3), CPFEM_ffn1(3,3)
real(pReal) Tstar_v_h(6), state_new_h(constitutive_Nstatevars), phi1, PHI, phi2, dt_i, delta_Fg(3,3), Fg_i(3,3) real(pReal) Tstar_v_h(6), state_new_h(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), phi1, PHI, phi2, dt_i
real(pReal) state_new_i(constitutive_Nstatevars), time real(pReal) delta_Fg(3,3), Fg_i(3,3), state_new_i(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), time
integer(pInt) jcut integer(pInt) jcut
! *** Numerical parameters *** ! *** Numerical parameters ***
integer(pInt), parameter :: ncut=7_pInt integer(pInt), parameter :: ncut=7_pInt
@ -404,14 +404,14 @@
! !
! *** Definition of variables *** ! *** Definition of variables ***
! *** Subroutine parameters *** ! *** Subroutine parameters ***
integer(pInt) cp_en, CPFEM_in, iori, ising, icut, iconv, isjaco
real(pReal) cs(6), dcs_de(6,6), dt, phi1, PHI, phi2, Fg_old(3,3), Fg_new(3,3) real(pReal) cs(6), dcs_de(6,6), dt, phi1, PHI, phi2, Fg_old(3,3), Fg_new(3,3)
real(pReal) Fp_old(3,3), Fp_new(3,3), state_old(constitutive_Nstatevars) real(pReal) Fp_old(3,3), Fp_new(3,3), state_old(constitutive_Nstatevars(iori, CPFEM_in, cp_en))
real(pReal) state_new(constitutive_Nstatevars), Tstar_v(6) real(pReal) state_new(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), Tstar_v(6)
integer(pInt) cp_en, CPFEM_in, iori, ising, icut, iconv, isjaco
! *** Local variables *** ! *** Local variables ***
integer(pInt) ic integer(pInt) ic
real(pReal) Fe(3,3), R(3,3), U(3,3), dev(6), dF(3,3), Fg2(3,3), sgm2(6) real(pReal) Fe(3,3), R(3,3), U(3,3), dev(6), dF(3,3), Fg2(3,3), sgm2(6)
real(pReal) state2(constitutive_Nstatevars), Fp2(3,3), cs1(6) real(pReal) state2(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), Fp2(3,3), cs1(6)
! *** Numerical parameters *** ! *** Numerical parameters ***
real(pReal), parameter :: pert_ct=1.0e-5_pReal real(pReal), parameter :: pert_ct=1.0e-5_pReal
! *** Error treatment *** ! *** Error treatment ***
@ -483,19 +483,19 @@
!*** NEWTON-RAPHSON Calculation *** !*** NEWTON-RAPHSON Calculation ***
!*********************************************************************** !***********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use constitutive, only: constitutive_Nstatevars use constitutive, only: constitutive_Nstatevars, constitutive_HomogenizedC, constitutive_dotState
use math use math
implicit none implicit none
! *** Definition of variables *** ! *** Definition of variables ***
! *** Subroutine parameters *** ! *** Subroutine parameters ***
real(pReal) dt,Fg_old(3,3),Fg_new(3,3),Fp_old(3,3),Fp_new(3,3), Fe(3,3)
real(pReal) state_old(constitutive_Nstatevars), state_new(constitutive_Nstatevars)
real(pReal) Tstar_v(6), cs(6)
integer(pInt) cp_en, CPFEM_in, iori, iconv, ising integer(pInt) cp_en, CPFEM_in, iori, iconv, ising
real(pReal) dt,Fg_old(3,3),Fg_new(3,3),Fp_old(3,3),Fp_new(3,3), Fe(3,3)
real(pReal) state_old(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), state_new(constitutive_Nstatevars(iori, CPFEM_in, cp_en))
real(pReal) Tstar_v(6), cs(6)
! *** Local variables *** ! *** Local variables ***
real(pReal) crite, tol_in, tol_out, invFp_old(3,3), det, A(3,3), C_66(6,6), Lp(3,3), dLp(3,3,6) real(pReal) crite, tol_in, tol_out, invFp_old(3,3), det, A(3,3), C_66(6,6), Lp(3,3), dLp(3,3,6)
real(pReal) I3tLp(3,3), help(3,3), help1(6), Tstar0_v(6), R1(6), norm1, tdLp(3,3) real(pReal) I3tLp(3,3), help(3,3), help1(6), Tstar0_v(6), R1(6), norm1, tdLp(3,3)
real(pReal) dstate(constitutive_Nstatevars), R2(6), norm2, invFp_new(3,3), Estar(3,3) real(pReal) dstate(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), R2(6), norm2, invFp_new(3,3), Estar(3,3)
real(pReal) Estar_v(6), Jacobi(6,6), invJacobi(6,6), dTstar_v(6) real(pReal) Estar_v(6), Jacobi(6,6), invJacobi(6,6), dTstar_v(6)
integer(pInt) iouter, iinner , dummy, err, i, j, k integer(pInt) iouter, iinner , dummy, err, i, j, k
! *** Numerical parameters *** ! *** Numerical parameters ***
@ -505,11 +505,11 @@
real(pReal), parameter :: tol_inner = 1.0e-3_pReal real(pReal), parameter :: tol_inner = 1.0e-3_pReal
real(pReal), parameter :: eta = 13.7_pReal real(pReal), parameter :: eta = 13.7_pReal
crite=eta*constitutive_s0_slip/constitutive_n_slip !ÄÄÄ ! crite=eta*constitutive_s0_slip/constitutive_n_slip !ÄÄÄ
! !
! *** Tolerances *** ! *** Tolerances ***
tol_in = tol_inner*s0_slip !ÄÄÄ ! tol_in = tol_inner*s0_slip !ÄÄÄ
tol_out = tol_outer*s0_slip !ÄÄÄ ! tol_out = tol_outer*s0_slip !ÄÄÄ
! !
! *** Error treatment *** ! *** Error treatment ***
iconv = 0 iconv = 0
@ -527,7 +527,7 @@
! *** Calculation of A and T*0 (see Kalidindi) *** ! *** Calculation of A and T*0 (see Kalidindi) ***
A = matmul(Fg_new,invFp_old) A = matmul(Fg_new,invFp_old)
A = matmul(transpose(A), A) A = matmul(transpose(A), A)
C_66=constitutive_homogenizedC(iori, CPFEM_in, cp_en) !ÄÄÄ C_66=constitutive_HomogenizedC(iori, CPFEM_in, cp_en) !ÄÄÄ
! !
! *** Second level of iterative procedure: Resistences *** ! *** Second level of iterative procedure: Resistences ***
do iouter=1,nouter do iouter=1,nouter
@ -554,7 +554,7 @@
enddo enddo
enddo enddo
enddo enddo
Jacobi= matmul(C_66, help1) + math_identity(6) ! Jacobi= matmul(C_66, help1) + math_identity(6)
call math_invert6x6(Jacobi, invJacobi, dummy, err) !ÄÄÄ call math_invert6x6(Jacobi, invJacobi, dummy, err) !ÄÄÄ
if (err==1_pInt) then if (err==1_pInt) then
forall (i=1:6) Jacobi(i,i)=1.05d0*maxval(Jacobi(i,:)) ! regularization forall (i=1:6) Jacobi(i,i)=1.05d0*maxval(Jacobi(i,:)) ! regularization

241
trunk/mpie_cpfem_marc.f90
View File

@ -4,7 +4,7 @@
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
! MPI fuer Eisenforschung, Duesseldorf ! MPI fuer Eisenforschung, Duesseldorf
! !
! last modified: 26.03.2007 ! last modified: 28.03.2007
!******************************************************************** !********************************************************************
! Usage: ! Usage:
! - choose material as hypela2 ! - choose material as hypela2
@ -32,102 +32,91 @@
include "mesh.f90" include "mesh.f90"
include "constitutive.f90" include "constitutive.f90"
include "CPFEM.f90" include "CPFEM.f90"
! !
subroutine hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kc,mats,ndi,nshear,& subroutine hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,&
disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,& nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,&
frotn1,strechn1,eigvn1,ncrd1,itel,ndeg1,ndm,& frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,&
nnode,jtype,lclass,ifr,ifu) nnode,jtype,lclass,ifr,ifu)
!******************************************************************** !********************************************************************
! This is the Marc material routine ! This is the Marc material routine
!******************************************************************** !********************************************************************
! !
! ************* user subroutine for defining material behavior ************** ! ************* user subroutine for defining material behavior **************
! !
! !
!CAUTION: Due to calculation of the Deformation gradients, Stretch Tensors and ! CAUTION : Due to calculation of the Deformation gradients, Stretch Tensors and
! Rotation tensors at previous and current states, the analysis can be ! Rotation tensors at previous and current states, the analysis can be
! computationally expensive. Please use the user subroutine -> hypela ! computationally expensive. Please use the user subroutine -> hypela
! if these kinematic quantities are not needed in the constitutive model ! if these kinematic quantities are not needed in the constitutive model
! !
! !
! IMPORTANT NOTES : ! IMPORTANT NOTES :
! !
! ! (1) F,R,U are only available for continuum and membrane elements (not for
! (1) You must include the -> process,1,1,1, card in the parameter section ! shells and beams).
! of MARC input deck. !
! ! (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(=
! (2) For total Lagrangian formulation use the -> 'large disp' card in the ! total Lagrange with large disp) in the parameter section of input deck.
! parameter section of MARC input deck. ! For updated Lagrangian formulation use the -> 'Plasticity,3' card(=
! For updated Lagrangian formulation use the -> 'large disp' and 'update' ! update+finite+large disp+constant d) in the parameter section of
! cards in the parameter section of MARC input deck. However for any ! input deck.
! large strain calculation (whether elasticity or inelasticity) must entail !
! the use of 'finite' parameter card also. !
! ! d stress strain law to be formed
! (3) For Plasticity, the 2nd or 3rd cards in 'geometry' option in the model ! g change in stress due to temperature effects
! definition sections must be flagged for correct behavior in incompressible ! e total elastic strain
! deformation. ! de increment of strain
! ! s stress - should be updated by user
! (4) The kinematic quantities are calculated for the following continuum ! t state variables (comes in at t=n, must be updated
! elements (both lower and higher order) : ! to have state variables at t=n+1)
! plane stress, plane strain, generalized plane strain, axisymmetric, ! dt increment of state variables
! axisymmetric with twist and brick elements. ! ngens size of stress - strain law
! ! n element number
! ! nn integration point number
! ! kcus(1) layer number
! ! kcus(2) internal layer number
! d stress strain law to be formed ! matus(1) user material identification number
! g change in stress due to temperature effects ! matus(2) internal material identification number
! e total strain ! ndi number of direct components
! de increment of strain ! nshear number of shear components
! s stress - should be updated by user ! disp incremental displacements
! t state variables (comes in at t=n, must be updated ! dispt displacements at t=n (at assembly, lovl=4) and
! to have state variables at t=n+1) ! displacements at t=n+1 (at stress recovery, lovl=6)
! dt increment of state variables ! coord coordinates
! ngens size of stress - strain law ! ncrd number of coordinates
! n element number ! ndeg number of degrees of freedom
! nn integration point number ! itel dimension of F and R, either 2 or 3
! kc layer number ! nnode number of nodes per element
! mats material identification number ! jtype element type
! ndi number of direct components ! lclass element class
! nshear number of shear components ! ifr set to 1 if R has been calculated
! disp incremental displacements ! ifu set to 1 if strech has been calculated
! dispt displacements at t=n (at assembly, lovl=4) and !
! displacements at t=n+1 (at stress recovery, lovl=6) ! at t=n :
! coord coordinates !
! ncrd number of coordinates ! ffn deformation gradient
! ndeg number of degrees of freedom ! frotn rotation tensor
! itel dimension of F and R, either 2 or 3 ! strechn square of principal stretch ratios, lambda(i)
! nnode number of nodes per element ! eigvn(i,j) i principal direction components for j eigenvalues
! jtype element type !
! lclass element class ! at t=n+1 :
! ifr set to 1 if R has been calculated !
! ifu set to 1 if strech has been calculated ! ffn1 deformation gradient
! ! frotn1 rotation tensor
! at t=n : ! strechn1 square of principal stretch ratios, lambda(i)
! ! eigvn1(i,j) i principal direction components for j eigenvalues
! ffn deformation gradient !
! frotn rotation tensor ! The following operation obtains U (stretch tensor) at t=n+1 :
! strechn square of principal stretch ratios, lambda(i) !
! eigvn(i,j) i principal direction components for j eigenvalues ! call scla(un1,0.d0,itel,itel,1)
! ! do 3 k=1,3
! at t=n+1 : ! do 2 i=1,3
! ! do 1 j=1,3
! ffn1 deformation gradient ! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k)
! frotn1 rotation tensor !1 continue
! strechn1 square of principal stretch ratios, lambda(i) !2 continue
! eigvn1(i,j) i principal direction components for j eigenvalues !3 continue
! !
! The following operation obtains U (stretch tensor) at t=n+1 :
!
! call scla(un1,0.d0,itel,itel,1)
! do 3 k=1,3
! do 2 i=1,3
! do 1 j=1,3
! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k)
!1 continue
!2 continue
!3 continue
!
! !
use prec, only: pReal,pInt use prec, only: pReal,pInt
use CPFEM, only : CPFEM_stress_all, CPFEM_jaco_old use CPFEM, only : CPFEM_stress_all, CPFEM_jaco_old
@ -162,8 +151,8 @@
real(pReal) mpie_timefactor, mpie_stress(ngens) real(pReal) mpie_timefactor, mpie_stress(ngens)
real(pReal) mpie_jacobi(ngens,ngens) real(pReal) mpie_jacobi(ngens,ngens)
! !
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,ngens),s(ngens),n(2),coord(ncrd,*),disp(ndeg,*),dispt(ndeg,*),ffn(itel,*),& dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),&
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*) frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2)
! !
! call general material routine only in increment 0 and for lovl==6 (stress recovery) ! call general material routine only in increment 0 and for lovl==6 (stress recovery)
@ -182,13 +171,13 @@
! mpie_in intergration point number ! mpie_in intergration point number
! mpie_dimension dimension of stress/strain vector ! mpie_dimension dimension of stress/strain vector
!******************************************************************** !********************************************************************
cp_en=mesh_mapFEtoCPelement(n(1)) cp_en=mesh_FEasCP('elem', n(1))
if ((lovl==6).or.(inc==0)) then if ((lovl==6).or.(inc==0)) then
call cpfem_general(ffn, ffn1, inc, incsub, ncycle, timinc, cp_en, nn) call cpfem_general(ffn, ffn1, inc, incsub, ncycle, timinc, cp_en, nn)
endif endif
! return stress and jacobi ! return stress and jacobi
s=CPFEM_stress_all(1:ngens, nn, cp_en) s(1:ngens)=CPFEM_stress_all(1:ngens, nn, cp_en)
d=CPFEM_jaco_old(1:ngens,1:ngens, nn, cp_en) d(1:ngens,1:ngens)=CPFEM_jaco_old(1:ngens,1:ngens, nn, cp_en)
return return
end end
@ -217,6 +206,8 @@
!******************************************************************** !********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use CPFEM, only: CPFEM_results, CPFEM_Nresults use CPFEM, only: CPFEM_results, CPFEM_Nresults
use constitutive, only: constitutive_maxNresults
use mesh, only: mesh_FEasCP
implicit none implicit none
! !
real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*) real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*)
@ -224,23 +215,23 @@
integer(pInt) m, nn, layer, ndi, nshear, jpltcd integer(pInt) m, nn, layer, ndi, nshear, jpltcd
! !
! assign result variable ! assign result variable
v=CPFEM_result(mod(jpltcd, CPFEM_Nresults+constitutive_Nresults),& v=CPFEM_results(mod(jpltcd, CPFEM_Nresults+constitutive_maxNresults),&
int(jpltcd/(CPFEM_Nresults+constitutive_Nresults)),& int(jpltcd/(CPFEM_Nresults+constitutive_maxNresults)),&
nn, mesh_mapFEtoCPelement(m)) nn, mesh_FEasCP('elem', m))
return return
end end
! !
! !
subroutine utimestep(timestep,timestepold,icall,time,timeloadcase) ! subroutine utimestep(timestep,timestepold,icall,time,timeloadcase)
!******************************************************************** !********************************************************************
! This routine modifies the addaptive time step of Marc ! This routine modifies the addaptive time step of Marc
!******************************************************************** !********************************************************************
use prec, only: pReal,pInt ! use prec, only: pReal,pInt
use CPFEM, only : CPFEM_timefactor_max ! use CPFEM, only : CPFEM_timefactor_max
implicit none ! implicit none
! !
real(pReal) timestep, timestepold, time,timeloadcase ! real(pReal) timestep, timestepold, time,timeloadcase
integer(pInt) icall ! integer(pInt) icall
! !
! user subroutine for modifying the time step in auto step ! user subroutine for modifying the time step in auto step
! !
@ -261,22 +252,22 @@
! !
! user coding ! user coding
! reduce timestep during increment in case mpie_timefactor is too large ! reduce timestep during increment in case mpie_timefactor is too large
if(icall==2_pInt) then ! if(icall==2_pInt) then
if(mpie_timefactor_max>1.25_pReal) then ! if(mpie_timefactor_max>1.25_pReal) then
timestep=min(timestep,timestepold*0.8_pReal) ! timestep=min(timestep,timestepold*0.8_pReal)
end if ! end if
return ! return
! modify timestep at beginning of new increment ! modify timestep at beginning of new increment
else if(icall==3_pInt) then ! else if(icall==3_pInt) then
if(mpie_timefactor_max<=0.8_pReal) then ! if(mpie_timefactor_max<=0.8_pReal) then
timestep=min(timestep,timestepold*1.25_pReal) ! timestep=min(timestep,timestepold*1.25_pReal)
else if (mpie_timefactor_max<=1.0_pReal) then ! else if (mpie_timefactor_max<=1.0_pReal) then
timestep=min(timestep,timestepold/mpie_timefactor_max) ! timestep=min(timestep,timestepold/mpie_timefactor_max)
else if (mpie_timefactor_max<=1.25_pReal) then ! else if (mpie_timefactor_max<=1.25_pReal) then
timestep=min(timestep,timestepold*1.01_pReal) ! timestep=min(timestep,timestepold*1.01_pReal)
else ! else
timestep=min(timestep,timestepold*0.8_pReal) ! timestep=min(timestep,timestepold*0.8_pReal)
end if ! end if
end if ! end if
return ! return
end ! end