OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps187.F File Reference
#include "implicit_f.inc"
#include "com01_c.inc"
#include "param_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps187 (nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, yld, pla, dpla, etse, ipm, mat, israte, facyldi, epsp)

Function/Subroutine Documentation

◆ sigeps187()

subroutine sigeps187 ( integer nel,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
time,
timestep,
uparam,
rho0,
rho,
volume,
eint,
epspxx,
epspyy,
epspzz,
epspxy,
epspyz,
epspzx,
depsxx,
depsyy,
depszz,
depsxy,
depsyz,
depszx,
epsxx,
epsyy,
epszz,
epsxy,
epsyz,
epszx,
sigoxx,
sigoyy,
sigozz,
sigoxy,
sigoyz,
sigozx,
signxx,
signyy,
signzz,
signxy,
signyz,
signzx,
sigvxx,
sigvyy,
sigvzz,
sigvxy,
sigvyz,
sigvzx,
soundsp,
viscmax,
uvar,
off,
yld,
pla,
dpla,
etse,
integer, dimension(npropmi,*) ipm,
integer, dimension(nel) mat,
integer, dimension(nel) israte,
facyldi,
epsp )

Definition at line 31 of file sigeps187.F.

44C=======================================================================
45C BARLAT YLD2000 material model
46C=======================================================================
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51C-----------------------------------------------
52C G L O B A L P A R A M E T E R S
53C-----------------------------------------------
54#include "com01_c.inc"
55#include "param_c.inc"
56C-----------------------------------------------
57C C O M M O N
58C-----------------------------------------------
59C---------+---------+---+---+--------------------------------------------
60C VAR | SIZE |TYP| RW| DEFINITION
61C---------+---------+---+---+--------------------------------------------
62C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
63C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
64C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
65C---------+---------+---+---+--------------------------------------------
66C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
67C IFUNC | NFUNC | I | R | FUNCTION INDEX
68C NPF | * | I | R | FUNCTION ARRAY
69C TF | * | F | R | FUNCTION ARRAY
70C---------+---------+---+---+--------------------------------------------
71C TIME | 1 | F | R | CURRENT TIME
72C TIMESTEP| 1 | F | R | CURRENT TIME STEP
73C UVAR | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
74C RHO0 | NEL | F | R | INITIAL DENSITY
75C RHO | NEL | F | R | DENSITY
76C VOLUME | NEL | F | R | VOLUME
77C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
78C EPSPXX | NEL | F | R | STRAIN RATE XX
79C EPSPYY | NEL | F | R | STRAIN RATE YY
80C ... | | | |
81C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
82C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
83C ... | | | |
84C EPSXX | NEL | F | R | STRAIN XX
85C EPSYY | NEL | F | R | STRAIN YY
86C ... | | | |
87C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
88C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
89C ... | | | |
90C---------+---------+---+---+--------------------------------------------
91C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
92C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
93C ... | | | |
94C SIGVXX | NEL | F | W | VISCOUS STRESS XX
95C SIGVYY | NEL | F | W | VISCOUS STRESS YY
96C ... | | | |
97C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
98C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
99C---------+---------+---+---+--------------------------------------------
100C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
101C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
102C---------+---------+---+---+--------------------------------------------
103C-----------------------------------------------
104C I N P U T A r g u m e n t s
105C-----------------------------------------------
106C
107 INTEGER :: NEL, NUPARAM, NUVAR
108 INTEGER :: IPM(NPROPMI,*),ISRATE(NEL),MAT(NEL)
109 my_real
110 . time,timestep,uparam(nuparam),
111 . rho(nel),rho0(nel),volume(nel),eint(nel),
112 . epspxx(nel),epspyy(nel),epspzz(nel),
113 . epspxy(nel),epspyz(nel),epspzx(nel),
114 . depsxx(nel),depsyy(nel),depszz(nel),
115 . depsxy(nel),depsyz(nel),depszx(nel),
116 . epsxx(nel) ,epsyy(nel) ,epszz(nel) ,
117 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
118 . sigoxx(nel),sigoyy(nel),sigozz(nel),
119 . sigoxy(nel),sigoyz(nel),sigozx(nel),
120 . facyldi(nel),epsp(nel)
121C-----------------------------------------------
122C O U T P U T A r g u m e n t s
123C-----------------------------------------------
124 my_real
125 . signxx(nel),signyy(nel),signzz(nel),
126 . signxy(nel),signyz(nel),signzx(nel),
127 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
128 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
129 . soundsp(nel),viscmax(nel),yld(nel) ,
130 . pla(nel),dpla(nel),etse(nel)
131C-----------------------------------------------
132C I N P U T O U T P U T A r g u m e n t s
133C-----------------------------------------------
134 my_real
135 . uvar(nel,nuvar), off(nel)
136C-----------------------------------------------
137C VARIABLES FOR FUNCTION INTERPOLATION
138C-----------------------------------------------
139 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
140 my_real
141 . finter ,tf(*)
142 EXTERNAL finter
143c Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
144c Y : y = f(x)
145c X : x
146c DYDX : f'(x) = dy/dx
147c IFUNC(J): FUNCTION INDEX
148c J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
149c NPF,TF : FUNCTION PARAMETER
150C-----------------------------------------------
151C L o c a l V a r i a b l e s
152C-----------------------------------------------
153 INTEGER I,II,J,K,NRATE,JJ(NEL),J1,J2,N,NINDX,IFLAG,IFLAGSR,EXPA,EXPAM2,
154 . NMAX,INDX(NEL), IPOS1(NEL),ILEN1(NEL),IAD1(NEL),
155 . IPOS2(NEL),ILEN2(NEL),IAD2(NEL)
156 my_real
157 . y1(nel),y2(nel),dydx(nel),dydx1(nel),dydx2(nel),
158 . xpxx(nel),xpyy(nel),xpxy(nel),xpyz(nel),xpzx(nel),
159 . xppxx(nel),xppyy(nel),xppxy(nel),xppyz(nel),xppzx(nel),
160 . phip(nel),phipp(nel),
161 . deltap(nel), deltapp1(nel), deltapp2(nel), deltapp(nel),
162 . deplxx(nel),deplyy(nel),deplxy(nel),deplzz(nel),
163 . deelzz(nel),xp1(nel),xp2(nel),xpp1(nel),
164 . xpp2(nel),sig(nel),h(nel),sigtrxx(nel),sigtryy(nel),
165 . sigtrxy(nel),rate(nel),sigtryz(nel),sigtrzz(nel),
166 . sigtrzx(nel),deplyz(nel),deplzx(nel),ddeplyz(nel),ddeplzx(nel),
167 . ddeplxx(nel),ddeplyy(nel),ddeplxy(nel),ddeplzz(nel)
168 my_real
169 . e,nu,bulk,a1,a2,g,unsa,du,r,umr,g2,lamda,
170 . al1, al2 , al3 , al4 ,dsigpl1,dspl1,
171 . al5, al6 , al7 , al8 ,dsigpl2,dspl2,
172 . al9, al10, al11, al12,dsigpl3,dspl3,
173 . yscale,ceps,eps0,expn,dav,
174 . lpp11,lpp12,lp13,lpp21,lpp22,lpp23,
175 . lpp34,lpp45,lpp56,lpp13,lpp66,
176 . fct,fctp,bpp,cpp,dpp,fp1,fp2,fp3,fp4,fp5,
177 . kint,fpp1,fpp2,fpp3,dfp1,dfp2,dfp3,fpp4,fpp5,
178 . dfpp1,dfpp2,dfpp3,dp,ap,df1,df2,df3,df,dfp4,dfp5,dfp6,
179 . dsdepl1,dsdepl2,dsdepl3,normef,epspi,dfpp4,dfpp5,dfpp6,
180 . aswift ,epso,qvoce,beta,ko,alpha,nexp,df4,df5,df6,
181 . expv,kswift,kvoce,unsp,unsc,pla_i,yld_i,
182 . dswift,dvoce,dfepsp,yldp, p2,unspt,unsct,dfsr,
183 . tpp1, tpp2,tpp3,tp1,tp2,fiptest,fipptest
184
185 my_real
186 . hk(nel),nu_1mnu(nel),epsf(nel),plap(nel),f_epsp(nel),
187 . svm2(nel),yld2(nel),hi(nel),g3(nel),t_1pnu(nel),
188 . hkin,nu110,nu210,dezz,svm,dpla_i(nel),dr(nel),
189 . aa(nel),bb(nel),dpla_j(nel),u_1mnu(nel),
190 . pp(nel),qq(nel),fail(nel),svmo(nel),f_epspt(nel),
191 . nu11v(nel),nu21v3(nel),aanu11v(nel),bbnu21v3(nel),
192 . hi2(nel), es2(nel)
193 my_real
194 . yfac(nel,2),pui_1,pui_2
195 DATA nmax/5/
196c-----------------------------------------------
197c USER VARIABLES INITIALIZATION
198c-----------------------------------------------
199c Swift/Voce
200c k=alpha*k1(eps,p)+(1-alpha)*k2(eps,p)
201c k1=A*(eps,p+e0)**n
202c k2=Q*(1-exp(-b*eps,p))+sig0
203c-----------------------------------------------
204c UVAR1 =PLA
205c UVAR2 =YLD
206c UVAR3 =DPLA
207c UVAR4 =ETSE
208c UVAR5 =IPOS
209c UVAR6 =H
210c-----------------------------------------------
211C
212 yldp = zero
213 e = uparam(1)
214 nu = uparam(2)
215 g = uparam(23)
216 g2 = two *g
217 bulk = e/three/(one-two*nu)
218
219 al1 = uparam(3)
220 al2 = uparam(4)
221 al3 = uparam(5)
222 al4 = uparam(6)
223 al5 = uparam(7)
224 al6 = uparam(8)
225 al7 = uparam(9)
226 al8 = uparam(10)
227 al9 = uparam(11)
228 al10 = uparam(12)
229 al11 = uparam(13)
230 al12 = uparam(14)
231
232
233c
234 expa = nint(uparam(15) )
235 unsa = one/expa
236 expam2= expa - 2
237c
238 alpha = uparam(16)
239 nexp = uparam(17)
240 epso = uparam(18)
241 aswift = uparam(19)
242
243 nrate = nint(uparam(25))
244 lamda = uparam(25+2*nfunc+9)
245
246
247 qvoce = uparam(25+2*nrate+1)
248 beta = uparam(25+2*nrate+2)
249 ko = uparam(25+2*nrate+3)
250
251 unsp = uparam(25+2*nrate+4)
252 unsc = uparam(25+2*nrate+5)
253
254
255 unspt = uparam(25+2*nrate+6)
256 unsct = uparam(25+2*nrate+7)
257
258 iflagsr = nint(uparam(25+2*nrate+8))
259 !DEPLACER AU STARTER
260 lpp11=(-two* al3 +two*al4 + eight*al5- two*al6)/nine
261 lpp12=(-four*al4 +four*al6+ al3 - four*al5)/nine
262 lpp13=(al3 +two*al4 -four*al5- two*al6 )/nine
263 lpp21=(four*al3 -four*al4-four*al5+ al6)/nine
264 lpp22=(-two* al3+ eight*al4 + two*al5- two*al6)/nine
265 lpp23=(-two* al3 -four*al4 + two*al5+ al6)/nine
266 lpp34=al8
267 lpp45=al11
268 lpp56=al12
269
270
271 iflag = int(uparam(24))
272
273c ---
274 DO i=1,nel
275 dav = depsxx(i)+depsyy(i)+depszz(i)
276 signxx(i) = sigoxx(i) + depsxx(i)*g2 + lamda*dav
277 signyy(i) = sigoyy(i) + depsyy(i)*g2 + lamda*dav
278 signzz(i) = sigozz(i) + depszz(i)*g2 + lamda*dav
279 signxy(i) = sigoxy(i) + depsxy(i)*g
280 signyz(i) = sigoyz(i) + depsyz(i)*g
281 signzx(i) = sigozx(i) + depszx(i)*g
282
283C------------------
284C SIGMA TRIAL
285C------------------
286 sigtrxx(i)=signxx(i)
287 sigtryy(i)=signyy(i)
288 sigtrzz(i)=signzz(i)
289 sigtrxy(i)=signxy(i)
290 sigtryz(i)=signyz(i)
291 sigtrzx(i)=signzx(i)
292c-------------------
293c SOUND SPEED, TANGENT MODULUS
294c-------------------
295 soundsp(i) = sqrt((bulk+four_over_3*g)/rho0(i))
296 viscmax(i) = zero
297 etse(i) = one
298 dpla(i) = zero
299 ENDDO
300C-------------------
301c flag formulation ************************************************
302C-------------------
303
304 IF (iflag == 1) THEN ! ANALYTIC SWIFT-VOCE YIELD
305c ----- Yield ---------------
306 DO i=1,nel
307 expv = exp(-beta*pla(i))
308 IF((pla(i) + epso)>zero) THEN
309 pui_1 = exp(nexp*log((pla(i) + epso)))
310 ELSE
311 pui_1 = zero
312 ENDIF
313! KSWIFT = ASWIFT*(PLA(I) + EPSO)**NEXP
314 kswift = aswift*pui_1
315 kvoce = ko + qvoce*(one - expv)
316 yld(i) = alpha*kswift + (one-alpha)*kvoce
317 f_epspt(i) = one
318 f_epsp(i) = one
319 !!print*, 'yld ', YLD(I)
320 ENDDO
321 IF (unspt /= zero )THEN
322 DO i=1,nel
323 IF(epsp(i)/=zero)
324 . f_epspt(i) =one + exp(unspt * log(unsct*epsp(i)))
325 ENDDO
326 ENDIF
327 DO i=1,nel
328 !--------------------
329 !BARLAT CRITERION
330 !--------------------
331 !XP(5) = LP(5*6) * SIGMA(6)
332 xpxx(i) = (two*al1*signxx(i)-al1*signyy(i)-al1*signzz(i))/three
333 xpyy(i) =(-al2*signxx(i)+two*al2*signyy(i)-al2*signzz(i))/three
334 xpxy(i) = al7*signxy(i)
335 xpyz(i) = al9*signyz(i)
336 xpzx(i) = al10*signzx(i)
337
338 !XPP(5) = LPP(5*6) * SIGMA(6)
339 xppxx(i) = lpp11*signxx(i) + lpp12*signyy(i) +lpp13*signzz(i)
340 xppyy(i) = lpp21*signxx(i) + lpp22*signyy(i) +lpp23*signzz(i)
341 xppxy(i) = lpp34*signxy(i)
342 xppyz(i) = lpp45*signyz(i)
343 xppzx(i) = lpp56*signzx(i)
344
345 deltap(i) = (xpxx(i)-xpyy(i))**2 + four * (xpxy(i)**2 + xpyz(i)**2 + xpzx(i)**2 )
346 phip(i) = sqrt(max(zero, deltap(i)))**expa
347
348
349
350 deltapp(i) = (xppxx(i)-xppyy(i))**2 + four * (xppxy(i)**2 + xppyz(i)**2 + xppzx(i)**2 )
351 deltapp(i) = sqrt(max(zero, deltapp(i)))
352 deltapp1(i) = three_half*(xppxx(i)-xppyy(i))+half*deltapp(i)
353 deltapp2(i) = three_half*(xppxx(i)-xppyy(i))-half*deltapp(i)
354
355 phipp(i)= deltapp1(i)**expa + deltapp2(i)**expa
356
357 !SIG(I) = (HALF*(PHIP(I)+PHIPP(I)))**UNSA
358 IF(((phip(i)+phipp(i)))>zero) THEN
359 sig(i) = exp(unsa*log(half*(phip(i)+phipp(i))))
360 ELSE
361 sig(i) = zero
362 ENDIF
363
364 ENDDO
365
366C Check Yield Condition
367C ----------------------------------------------------
368 nindx = 0
369 DO i=1,nel
370 IF (sig(i)>yld(i).AND.off(i) == one)THEN ! Plastic Loading
371 nindx = nindx + 1
372 indx(nindx) = i
373 ENDIF
374 ENDDO
375 !--------------------
376 !PLASTIC FLOW
377 !--------------------
378 DO ii=1,nindx
379 i = indx(ii)
380 !IF (SIG(I)>YLD(I).AND.OFF(I) == ONE)THEN
381 dpla(i)= uvar(i,3)+ em09
382
383 !DPLA(I)= (SIG(I)-YLD(I))/(THREE*G+DYDX(I))
384 deplxx(i)=zero
385 deplyy(i)=zero
386 deplzz(i)=zero
387 deplxy(i)=zero
388 deplyz(i)=zero
389 deplzx(i)=zero
390
391 !YLD_I=YLD(I)
392 pla_i=pla(i)
393
394 DO k=1,4
395 !algo newton
396 !calcul de dsigma/ddpla=dphi/dsigma*dsigma/ddpla
397 !phi = phip + phipp
398 !dphip/dsigma = dphip/dxp*dxp/dsigma where dxp/dsigma = Lp
399 pla(i) = pla_i + dpla(i)
400 plap(i)= dpla(i)/timestep
401C
402 expv = exp(-beta*pla(i))
403! kswift = aswift*(pla(i) + epso)**nexp
404 IF((pla(i) + epso)>zero) THEN
405 pui_1 = exp(nexp*log(pla(i) + epso))
406 ELSE
407 pui_1 = zero
408 ENDIF
409 kswift = aswift*pui_1
410 kvoce = ko + qvoce*(one - expv)
411 yld(i) = alpha*kswift + (one-alpha)*kvoce
412
413 !===========ADD PL SR DEPENDENCY THROUGH COWPER SEYMONDS
414 !Compute plastic strain rate Cowper seymonds
415 IF( unsp/=zero)THEN
416 IF(plap(i)>zero)THEN
417 f_epsp(i) = one+exp(unsp * log(unsc*plap(i)))
418 ELSE
419 f_epsp(i) = one
420 ENDIF
421 ELSE
422 f_epsp(i) = one
423 ENDIF
424 !DERIVATIVE OF YLD
425 dswift = nexp*kswift/(pla(i) + epso)
426 dvoce = qvoce*beta*expv
427 !DFEPSP = UNSP *F_EPSP(I)/ (DPLA(I)+TIMESTEP/MAX(EM20,UNSC))
428 dfepsp = unsp *(f_epsp(i)-one)/ max(em20,dpla(i))
429 yldp = alpha*dswift + (one-alpha)*dvoce
430 yldp = yldp*f_epsp(i) + yld(i)*dfepsp
431 yld(i) = yld(i) * f_epsp(i) *f_epspt(i)
432 !===========END PL SR DEPENDENCY
433 fct = sig(i)- yld(i) ! residu
434
435 tp1 = expa*(sqrt(max(zero,deltap(i))))**(expa-2)
436 !D_phip/D_XP
437 fp1 = tp1 * (xpxx(i) - xpyy(i))
438 fp2 = -fp1
439 fp3 = three * tp1 * xpxy(i)
440 fp4 = three * tp1 * xpyz(i)
441 fp5 = three * tp1 * xpzx(i)
442
443
444 tpp1 = expa * deltapp1(i) ** (expa - 1)
445 tpp2 = expa * deltapp2(i) ** (expa - 1)
446 tpp3 = two *(tpp1 - tpp2) /deltapp(i)
447
448 !D_phipp/D_XPP
449 fpp1 = tpp1 *(three_half+half* ( xppxx(i) - xppyy(i) )/deltapp(i))
450 . + tpp2 *(three_half-half* ( xppxx(i) - xppyy(i) )/deltapp(i))
451 fpp2 = tpp1 *(-three_half-half* ( xppxx(i) - xppyy(i) )/deltapp(i))
452 . + tpp2 *(-three_half+half* ( xppxx(i) - xppyy(i) )/deltapp(i))
453 fpp3 = tpp3 * xppxy(i)
454 fpp4 = tpp3 * xppyz(i)
455 fpp5 = tpp3 * xppzx(i)
456
457 dfp1=third*(two*al1*fp1-al2*fp2)
458 dfp2=third*(two*al2*fp2-al1*fp1)
459 dfp3=third*(-al1*fp1-al2*fp2)
460 dfp4=al7*fp3
461 dfp5=al9*fp4
462 dfp6=al10*fp5
463
464 dfpp1=fpp1*lpp11+fpp2*lpp21
465 dfpp2=fpp1*lpp12+fpp2*lpp22
466 dfpp3=fpp1*lpp13+fpp2*lpp23
467 dfpp4=fpp3*lpp34
468 dfpp5=fpp4*lpp45
469 dfpp6=fpp5*lpp56
470
471 du = unsa * sig(i) / (max(em20, phip(i)+phipp(i)) )
472
473
474 df1 = dfp1+dfpp1
475 df2 = dfp2+dfpp2
476 df3 = dfp3+dfpp3
477 df4 = dfp4+dfpp4
478 df5 = dfp5+dfpp5
479 df6 = dfp6+dfpp6
480
481 !DF = DF1*DSDEPL1+DF2*DSDEPL2+DSDEPL3*DF3
482 ddeplxx(i)=df1*du
483 ddeplyy(i)=df2*du
484 ddeplzz(i)=df3*du
485 ddeplxy(i)=df4*du
486 ddeplyz(i)=df5*du
487 ddeplzx(i)=df6*du
488
489 df = - g * du**2
490 . * (two*(df1*df1 + df2*df2 + df3*df3 )
491 . + df4*df4 + df5*df5 + df6*df6 )
492
493c DF = -TWO*G * (DU * DF1) **2
494c . -TWO*G * (DU * DF2) **2
495c . -TWO*G * (DU * DF3) **2
496c . -G * (DU * DF4) **2
497c . -G * (DU * DF5) **2
498c . -G * (DU * DF6) **2
499
500 df = df - yldp
501 !CALCUL PAS PLASTIQUE
502 dpla(i)=dpla(i)-fct/df
503 dpla(i)=max(zero,dpla(i))
504
505 !TENSEUR PLASTIQUE
506 deplxx(i)=dpla(i)*df1*du
507 deplyy(i)=dpla(i)*df2*du
508 deplzz(i)=dpla(i)*df3*du
509 deplxy(i)=dpla(i)*df4*du
510 deplyz(i)=dpla(i)*df5*du
511 deplzx(i)=dpla(i)*df6*du
512
513 !ACTUALISATION DES CONTRAINTES
514 signxx(i)=sigtrxx(i)- g2*deplxx(i)
515 signyy(i)=sigtryy(i)- g2*deplyy(i)
516 signzz(i)=sigtrzz(i)- g2*deplzz(i)
517 signxy(i)=sigtrxy(i)- g*deplxy(i)
518 signyz(i)=sigtryz(i)- g*deplyz(i)
519 signzx(i)=sigtrzx(i)- g*deplzx(i)
520 !CALCUL NOUVEAU SIGMA EFFECTIF
521 !--------------------
522 !BARLAT CRITERION
523 !--------------------
524 xpxx(i) = (two*al1*signxx(i)-al1*signyy(i)-al1*signzz(i))/three
525 xpyy(i) =(-al2*signxx(i)+two*al2*signyy(i)-al2*signzz(i))/three
526 xpxy(i) = al7*signxy(i)
527 xpyz(i) = al9*signyz(i)
528 xpzx(i) = al10*signzx(i)
529
530 xppxx(i) = lpp11*signxx(i) + lpp12*signyy(i) +lpp13*signzz(i)
531 xppyy(i) = lpp21*signxx(i) + lpp22*signyy(i) +lpp23*signzz(i)
532 xppxy(i) = lpp34*signxy(i)
533 xppyz(i) = lpp45*signyz(i)
534 xppzx(i) = lpp56*signzx(i)
535
536 deltap(i) = (xpxx(i)-xpyy(i))**2 + four * (xpxy(i)**2 + xpyz(i)**2 + xpzx(i)**2 )
537
538 phip(i) = sqrt(max(zero, deltap(i)))**expa
539
540
541 deltapp(i) = (xppxx(i)-xppyy(i))**2 + four * (xppxy(i)**2 + xppyz(i)**2 + xppzx(i)**2 )
542 deltapp(i) = sqrt(max(em20, deltapp(i)))
543
544 deltapp1(i) = three_half*(xppxx(i)-xppyy(i))+half*deltapp(i)
545 deltapp2(i) = three_half*(xppxx(i)-xppyy(i))-half*deltapp(i)
546
547 phipp(i)= deltapp1(i)**expa + deltapp2(i)**expa
548
549! SIG(I) = (HALF*(PHIP(I)+PHIPP(I)))**UNSA
550 IF(((phip(i)+phipp(i)))>zero) THEN
551 sig(i) = exp(unsa*log(half*(phip(i)+phipp(i))))
552 ELSE
553 sig(i) = zero
554 ENDIF
555
556 ENDDO ! iter
557 uvar(i,3)=dpla(i)
558 ENDDO ! NINDX
559C-------------------
560C-------------------
561 ELSE ! IFLAG *TABULATED YIELD
562C-------------------
563c yield tabulated in below
564C-------------------
565 IF(iflagsr == 0)THEN
566C-------------------
567C-------------------
568 DO i=1,nel
569 jj(i) = 1
570 DO j=2,nrate-1
571 IF (epsp(i) > uparam(25+j)) jj(i) = j
572 ENDDO
573 ENDDO
574 DO i=1,nel
575 epspi = uparam(25+jj(i)) !
576 rate(i)=(epsp(i) - epspi)/(uparam(26+jj(i)) - epspi)
577 yfac(i,1)=uparam(25+nrate+jj(i))*facyldi(i)
578 yfac(i,2)=uparam(25+nrate+jj(i)+1)*facyldi(i)
579 ENDDO
580 DO i=1,nel
581 j1 = jj(i)
582 j2 = j1+1
583 ipos1(i) = nint(uvar(i,j1+5))
584 iad1(i) = npf(ipm(10+j1,mat(1))) / 2 + 1
585 ilen1(i) = npf(ipm(10+j1,mat(1))+1) / 2 - iad1(i)-ipos1(i)
586 ipos2(i) = nint(uvar(i,j2+5))
587 iad2(i) = npf(ipm(10+j2,mat(1))) / 2 + 1
588 ilen2(i) = npf(ipm(10+j2,mat(1))+1) / 2 - iad2(i)-ipos2(i)
589 END DO
590 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
591 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
592 DO i=1,nel
593 j1 = jj(i)
594 j2 = j1+1
595 y1(i)=y1(i)*yfac(i,1)
596 y2(i)=y2(i)*yfac(i,2)
597 yld(i) = (y1(i) + rate(i)*(y2(i)-y1(i)))
598 yld(i) = max(yld(i),em20)
599 dydx1(i)=dydx1(i)*yfac(i,1)
600 dydx2(i)=dydx2(i)*yfac(i,2)
601 dydx(i) = (dydx1(i) + rate(i)*(dydx2(i)-dydx1(i)))
602 uvar(i,j1+5) = ipos1(i)
603 uvar(i,j2+5) = ipos2(i)
604 ENDDO
605c------------------------------
606 ELSE !flagSR
607c------------------------------
608 IF (isigi == 0) THEN
609 IF(time == zero)THEN
610 DO i=1,nel
611 yfac(i,1)= uparam(25+nrate+1)*facyldi(i)
612 yld(i) = yfac(i,1)*finter(ipm(11,mat(1)) ,zero,npf,tf,dydx(i))
613 uvar(i,2)= yld(i)
614 ENDDO
615 ENDIF
616 ENDIF
617 DO i=1,nel
618 pla(i) = uvar(i,1)
619 yld(i) = uvar(i,2)
620 dydx(i) = uvar(i,5)
621 plap(i) = uvar(i,6)
622 jj(i) = 1
623 ENDDO
624 ENDIF !flagSR
625C-------------------
626C-------------------
627c------------------------------
628 DO i=1,nel
629 !--------------------
630 !BARLAT CRITERION
631 !--------------------
632 xpxx(i) = (two*al1*signxx(i)-al1*signyy(i)-al1*signzz(i))/three
633 xpyy(i) =(-al2*signxx(i)+two*al2*signyy(i)-al2*signzz(i))/three
634 xpxy(i) = al7*signxy(i)
635 xpyz(i) = al9*signyz(i)
636 xpzx(i) = al10*signzx(i)
637
638 xppxx(i) = lpp11*signxx(i) + lpp12*signyy(i) +lpp13*signzz(i)
639 xppyy(i) = lpp21*signxx(i) + lpp22*signyy(i) +lpp23*signzz(i)
640 xppxy(i) = lpp34*signxy(i)
641 xppyz(i) = lpp45*signyz(i)
642 xppzx(i) = lpp56*signzx(i)
643
644 deltap(i) = (xpxx(i)-xpyy(i))**2 + four * (xpxy(i)**2 + xpyz(i)**2 + xpzx(i)**2 )
645 phip(i) = sqrt(max(zero, deltap(i)))**expa
646
647
648 deltapp(i) = (xppxx(i)-xppyy(i))**2 + four * (xppxy(i)**2 + xppyz(i)**2 + xppzx(i)**2 )
649 deltapp(i) = sqrt(max(em20, deltapp(i)))
650 deltapp1(i) = three_half*(xppxx(i)-xppyy(i))+half*deltapp(i)
651 deltapp2(i) = three_half*(xppxx(i)-xppyy(i))-half*deltapp(i)
652
653 phipp(i)= deltapp1(i)**expa + deltapp2(i)**expa
654
655! SIG(I) = (HALF*(PHIP(I)+PHIPP(I)))**UNSA
656 IF((phip(i)+phipp(i))>zero) THEN
657 sig(i) = exp(unsa*log(half*(phip(i)+phipp(i))))
658 ELSE
659 sig(i) = zero
660 ENDIF
661 ENDDO
662C ----------------------------------------------------
663C Check Yield Condition
664C ----------------------------------------------------
665 nindx = 0
666 DO i=1,nel
667 IF (sig(i)>yld(i).AND.off(i) == one)THEN ! Plastic Loading
668 nindx = nindx + 1
669 indx(nindx) = i
670 ENDIF
671 ENDDO
672C ----------------------------------------------------
673
674C ----------------------------------------------------
675c flag strain rate option *****************************************
676C ----------------------------------------------------
677 IF(iflagsr == 0)THEN
678 DO ii=1,nindx
679 i = indx(ii)
680 !--------------------
681 !PLASTIC FLOW
682 !--------------------
683 !DPLA(I)= UVAR(I,3)
684 !DPLA(I)= (SIG(I)-YLD(I))/(THREE*G+DYDX(I))
685 dpla(i)= uvar(i,3)+ em09
686 deplxx(i)=zero
687 deplyy(i)=zero
688 deplzz(i)=zero
689 deplxy(i)=zero
690 deplyz(i)=zero
691 deplzx(i)=zero
692
693 yld_i=yld(i)
694 DO k=1,5
695 !algo newton
696 !calcul de dsigma/ddpla=dphi/dsigma*dsigma/ddpla
697 !phi = phip + phipp
698 !dphip/dsigma = dphip/dxp*dxp/dsigma where dxp/dsigma = Lp
699 yld(i) =yld_i+dydx(i)*dpla(i)
700 fct = sig(i)- yld(i) ! residu
701
702 tp1 = expa*(sqrt(max(zero,deltap(i))))**(expa-2)
703
704 !D_phip/D_XP
705 fp1 = tp1 * (xpxx(i) - xpyy(i))
706 fp2 = -fp1
707 fp3 = three * tp1 * xpxy(i)
708 fp4 = three * tp1 * xpyz(i)
709 fp5 = three * tp1 * xpzx(i)
710
711
712 tpp1 = expa * deltapp1(i) ** (expa - 1)
713 tpp2 = expa * deltapp2(i) ** (expa - 1)
714 tpp3 = two *(tpp1 - tpp2) /deltapp(i)
715
716 !D_phipp/D_XPP
717 fpp1 = tpp1 *(three_half+half* ( xppxx(i) - xppyy(i) )/deltapp(i))
718 . +tpp2 *(three_half-half* ( xppxx(i) - xppyy(i) )/deltapp(i))
719 fpp2 = tpp1 *(-three_half-half* ( xppxx(i) - xppyy(i) )/deltapp(i))
720 . +tpp2 *(-three_half+half* ( xppxx(i) - xppyy(i) )/deltapp(i))
721 fpp3 = tpp3 * xppxy(i)
722 fpp4 = tpp3 * xppyz(i)
723 fpp5 = tpp3 * xppzx(i)
724
725 dfp1=third*(two*al1*fp1-al2*fp2)
726 dfp2=third*(two*al2*fp2-al1*fp1)
727 dfp3=third*(-al1*fp1-al2*fp2)
728 dfp4=al7*fp3
729 dfp5=al9*fp4
730 dfp6=al10*fp5
731
732
733 dfpp1=fpp1*lpp11+fpp2*lpp21
734 dfpp2=fpp1*lpp12+fpp2*lpp22
735 dfpp3=fpp1*lpp13+fpp2*lpp23
736 dfpp4=fpp3*lpp34
737 dfpp5=fpp4*lpp45
738 dfpp6=fpp5*lpp56
739
740 du = unsa * sig(i) / (max(em20, phip(i)+phipp(i)) )
741
742 df1 = dfp1+dfpp1
743 df2 = dfp2+dfpp2
744 df3 = dfp3+dfpp3
745 df4 = dfp4+dfpp4
746 df5 = dfp5+dfpp5
747 df6 = dfp6+dfpp6
748
749 df = -two*g * (du * df1) **2
750 . -two*g * (du * df2) **2
751 . -two*g * (du * df3) **2
752 . -g * (du * df4) **2
753 . -g * (du * df5) **2
754 . -g * (du * df6) **2
755
756 df = df - dydx(i)
757
758 !CALCUL PAS PLASTIQUE
759 dpla(i)=dpla(i)-fct/df
760 dpla(i)=max(zero,dpla(i))
761 !TENSEUR PLASTIQUE
762 deplxx(i)=dpla(i)*df1*du
763 deplyy(i)=dpla(i)*df2*du
764 deplzz(i)=dpla(i)*df3*du
765 deplxy(i)=dpla(i)*df4*du
766 deplyz(i)=dpla(i)*df5*du
767 deplzx(i)=dpla(i)*df6*du
768
769 !ACTUALISATION DES CONTRAINTES
770 signxx(i)=sigtrxx(i)- g2*deplxx(i)
771 signyy(i)=sigtryy(i)- g2*deplyy(i)
772 signzz(i)=sigtrzz(i)- g2*deplzz(i)
773 signxy(i)=sigtrxy(i)- g*deplxy(i)
774 signyz(i)=sigtryz(i)- g*deplyz(i)
775 signzx(i)=sigtrzx(i)- g*deplzx(i)
776
777 !CALCUL NOUVEAU SIGMA EFFECTIF
778 !--------------------
779 !BARLAT CRITERION
780 !--------------------
781 xpxx(i) = (two*al1*signxx(i)-al1*signyy(i)-al1*signzz(i))/three
782 xpyy(i) =(-al2*signxx(i)+two*al2*signyy(i)-al2*signzz(i))/three
783 xpxy(i) = al7*signxy(i)
784 xpyz(i) = al9*signyz(i)
785 xpzx(i) = al10*signzx(i)
786
787 xppxx(i) = lpp11*signxx(i) + lpp12*signyy(i) +lpp13*signzz(i)
788 xppyy(i) = lpp21*signxx(i) + lpp22*signyy(i) +lpp23*signzz(i)
789 xppxy(i) = lpp34*signxy(i)
790 xppyz(i) = lpp45*signyz(i)
791 xppzx(i) = lpp56*signzx(i)
792
793 deltap(i) = (xpxx(i)-xpyy(i))**2 + four * (xpxy(i)**2 + xpyz(i)**2 + xpzx(i)**2 )
794 phip(i) = sqrt(max(zero, deltap(i)))**expa
795
796
797 deltapp(i) = (xppxx(i)-xppyy(i))**2 + four * (xppxy(i)**2 + xppyz(i)**2 + xppzx(i)**2 )
798 deltapp(i) = sqrt(max(em20, deltapp(i)))
799 deltapp1(i) = three_half*(xppxx(i)-xppyy(i))+half*deltapp(i)
800 deltapp2(i) = three_half*(xppxx(i)-xppyy(i))-half*deltapp(i)
801
802 phipp(i)= deltapp1(i)**expa + deltapp2(i)**expa
803
804! !SIG(I) = (HALF*(PHIP(I)+PHIPP(I)))**UNSA
805 IF((half*(phip(i)+phipp(i)))>zero) THEN
806 sig(i) = exp(unsa*log(half*(phip(i)+phipp(i))))
807 ELSE
808 sig(i) = zero
809 ENDIF
810
811 ENDDO !iter
812 uvar(i,3)=dpla(i)
813 ENDDO ! NINDX
814C=======================================================================
815 ELSE !IFLAGSR PLASTIC STRAIN RATE ******************************
816C=======================================================================
817 DO ii=1,nindx
818 i = indx(ii)
819 !--------------------
820 !PLASTIC FLOW
821 !--------------------
822 dpla(i) = uvar(i,3) + em09
823 plap(i) = dpla(i)/timestep
824 !--------------
825 !update yield :
826 !--------------
827 jj(i) = 1
828 DO j=2,nrate-1
829 IF(plap(i) >= uparam(25+j)) jj(i) = j
830 ENDDO
831 j1 = jj(i)
832 j2 = j1+1
833 rate(i)=(plap(i) - uparam(25+j1))/(uparam(26+j1) - uparam(25+j1))
834 yfac(i,1)=uparam(25+nrate+j1)*facyldi(i)
835 yfac(i,2)=uparam(25+nrate+j2)*facyldi(i)
836 !Y1 and Y2 for Pl SR dependency
837 y1(i) = yfac(i,1)*finter(ipm(10+j1,mat(1)),pla(i),npf,tf,dydx1(i))
838 y2(i) = yfac(i,2)*finter(ipm(10+j2,mat(1)),pla(i),npf,tf,dydx2(i))
839 yld(i) = y1(i) + rate(i)*(y2(i)-y1(i))
840 yld(i) = max(yld(i),em20)
841 dydx1(i)= dydx1(i)*yfac(i,1)
842 dydx2(i)= dydx2(i)*yfac(i,2)
843 dydx(i) = dydx1(i) + rate(i)*(dydx2(i)-dydx1(i))
844 yldp = dydx(i) + (y2(i)-y1(i))/
845 . (uparam(26+j1)-uparam(25+j1)) /timestep
846
847
848
849 !DPLA(I)= UVAR(I,3)
850 !DPLA(I)= (SIG(I)-YLD(I))/(THREE*G+DYDX(I))
851 deplxx(i)=zero
852 deplyy(i)=zero
853 deplzz(i)=zero
854 deplxy(i)=zero
855 deplyz(i)=zero
856 deplzx(i)=zero
857
858
859 DO k=1,5
860 !algo newton
861 !calcul de dsigma/ddpla=dphi/dsigma*dsigma/ddpla
862 !phi = phip + phipp
863 !dphip/dsigma = dphip/dxp*dxp/dsigma where dxp/dsigma = Lp
864 !---------------------------------------
865 !update yield :
866 !---------------------------------------
867 pla(i) = uvar(i,1) + dpla(i)
868 plap(i) = dpla(i) / timestep
869 jj(i) = 1
870 DO j=2,nrate-1
871 IF(plap(i) >= uparam(25+j)) jj(i) = j
872 ENDDO
873 j1 = jj(i)
874 j2 = j1+1
875 rate(i)=(plap(i) - uparam(25+j1))/(uparam(26+j1) - uparam(25+j1))
876 yfac(i,1)=uparam(25+nrate+j1)*facyldi(i)
877 yfac(i,2)=uparam(25+nrate+j2)*facyldi(i)
878 !Y1 and Y2 for Pl SR dependency
879 y1(i) = yfac(i,1)*finter(ipm(10+j1,mat(1)),pla(i),npf,tf,dydx1(i))
880 y2(i) = yfac(i,2)*finter(ipm(10+j2,mat(1)),pla(i),npf,tf,dydx2(i))
881 yld(i) = y1(i) + rate(i)*(y2(i)-y1(i))
882 yld(i) = max(yld(i),em20)
883 dydx1(i)=dydx1(i)*yfac(i,1)
884 dydx2(i)=dydx2(i)*yfac(i,2)
885 dydx(i) = (dydx1(i) + rate(i)*(dydx2(i)-dydx1(i)))
886 !---------------------------------------
887 yldp = dydx(i)+ (y2(i)-y1(i))
888 . /(uparam(26+j1)-uparam(25+j1)) /timestep
889
890 !---------------------------------------
891
892 !YLD(I) =YLD_I+DYDX(I)*DPLA(I)
893
894 fct = sig(i)- yld(i)
895
896 tp1 = expa*(sqrt(max(zero,deltap(i))))**(expa-2)
897
898 !d_phip/d_xp
899 fp1 = tp1 * (xpxx(i) - xpyy(i))
900 fp2 = -fp1
901 fp3 = three * tp1 * xpxy(i)
902 fp4 = three * tp1 * xpyz(i)
903 fp5 = three * tp1 * xpzx(i)
904
905
906 tpp1 = expa * deltapp1(i) ** (expa - 1)
907 tpp2 = expa * deltapp2(i) ** (expa - 1)
908 tpp3 = two *(tpp1 - tpp2) /deltapp(i)
909
910 !D_phipp/D_XPP
911 fpp1 = tpp1 *(three_half+half* ( xppxx(i) - xppyy(i) )/deltapp(i))
912 . +tpp2 *(three_half-half* ( xppxx(i) - xppyy(i) )/deltapp(i))
913 fpp2 = tpp1 *(-three_half-half* ( xppxx(i) - xppyy(i) )/deltapp(i))
914 . +tpp2 *(-three_half+half* ( xppxx(i) - xppyy(i) )/deltapp(i))
915 fpp3 = tpp3 * xppxy(i)
916 fpp4 = tpp3 * xppyz(i)
917 fpp5 = tpp3 * xppzx(i)
918
919 dfp1=third*(two*al1*fp1-al2*fp2)
920 dfp2=third*(two*al2*fp2-al1*fp1)
921 dfp3=third*(-al1*fp1-al2*fp2)
922 dfp4=al7*fp3
923 dfp5=al9*fp4
924 dfp6=al10*fp5
925
926
927 dfpp1=fpp1*lpp11+fpp2*lpp21
928 dfpp2=fpp1*lpp12+fpp2*lpp22
929 dfpp3=fpp1*lpp13+fpp2*lpp23
930 dfpp4=fpp3*lpp34
931 dfpp5=fpp4*lpp45
932 dfpp6=fpp5*lpp56
933
934 du = unsa * sig(i) / (max(em20, phip(i)+phipp(i)) )
935
936 df1 = dfp1+dfpp1
937 df2 = dfp2+dfpp2
938 df3 = dfp3+dfpp3
939 df4 = dfp4+dfpp4
940 df5 = dfp5+dfpp5
941 df6 = dfp6+dfpp6
942
943 df = -two*g * (du * df1) **2
944 . -two*g * (du * df2) **2
945 . -two*g * (du * df3) **2
946 . -g * (du * df4) **2
947 . -g * (du * df5) **2
948 . -g * (du * df6) **2
949
950 df = df - yldp
951
952 !CALCUL PAS PLASTIQUE
953 dpla(i)=dpla(i)-fct/df
954 dpla(i)=max(zero,dpla(i))
955 !TENSEUR PLASTIQUE
956 deplxx(i)=dpla(i)*df1*du
957 deplyy(i)=dpla(i)*df2*du
958 deplzz(i)=dpla(i)*df3*du
959 deplxy(i)=dpla(i)*df4*du
960 deplyz(i)=dpla(i)*df5*du
961 deplzx(i)=dpla(i)*df6*du
962
963 !ACTUALISATION DES CONTRAINTES
964 signxx(i)=sigtrxx(i)- g2*deplxx(i)
965 signyy(i)=sigtryy(i)- g2*deplyy(i)
966 signzz(i)=sigtrzz(i)- g2*deplzz(i)
967 signxy(i)=sigtrxy(i)- g*deplxy(i)
968 signyz(i)=sigtryz(i)- g*deplyz(i)
969 signzx(i)=sigtrzx(i)- g*deplzx(i)
970
971 !CALCUL NOUVEAU SIGMA EFFECTIF
972 !--------------------
973 !BARLAT CRITERION
974 !--------------------
975 xpxx(i) = (two*al1*signxx(i)-al1*signyy(i)-al1*signzz(i))/three
976 xpyy(i) =(-al2*signxx(i)+two*al2*signyy(i)-al2*signzz(i))/three
977 xpxy(i) = al7*signxy(i)
978 xpyz(i) = al9*signyz(i)
979 xpzx(i) = al10*signzx(i)
980
981 xppxx(i) = lpp11*signxx(i) + lpp12*signyy(i) +lpp13*signzz(i)
982 xppyy(i) = lpp21*signxx(i) + lpp22*signyy(i) +lpp23*signzz(i)
983 xppxy(i) = lpp34*signxy(i)
984 xppyz(i) = lpp45*signyz(i)
985 xppzx(i) = lpp56*signzx(i)
986
987 deltap(i) = (xpxx(i)-xpyy(i))**2 + four * (xpxy(i)**2 + xpyz(i)**2 + xpzx(i)**2 )
988 phip(i) = sqrt(max(zero, deltap(i)))**expa
989
990
991 deltapp(i) = (xppxx(i)-xppyy(i))**2 + four * (xppxy(i)**2 + xppyz(i)**2 + xppzx(i)**2 )
992 deltapp(i) = sqrt(max(em20, deltapp(i)))
993 deltapp1(i) = three_half*(xppxx(i)-xppyy(i))+half*deltapp(i)
994 deltapp2(i) = three_half*(xppxx(i)-xppyy(i))-half*deltapp(i)
995
996 phipp(i)= deltapp1(i)**expa + deltapp2(i)**expa
997
998! !SIG(I) = (HALF*(PHIP(I)+PHIPP(I)))**UNSA
999 IF((half*(phip(i)+phipp(i)))>zero) THEN
1000 sig(i) = exp(unsa*log(half*(phip(i)+phipp(i))))
1001 ELSE
1002 sig(i) = zero
1003 ENDIF
1004 ENDDO ! iter k
1005 uvar(i,3) =dpla(i)
1006 ENDDO !ndx
1007c----------------------------
1008 ENDIF !flag formulation
1009 ENDIF
1010C=======================================================================
1011 DO i=1,nel
1012C
1013 IF( dpla(i) > zero)THEN
1014 h(i)=(yld(i)-uvar(i,2))/dpla(i)
1015 h(i)=max(em10,h(i))
1016 uvar(i,2)=yld(i)
1017 uvar(i,4)=h(i)/g2!(H(I)+E)
1018 uvar(i,5)=h(i)
1019 ENDIF
1020 pla(i)=uvar(i,1) + dpla(i)
1021 uvar(i,1) = pla(i)
1022 etse(i)= uvar(i,4)
1023 yld(i) = uvar(i,2)
1024 h(i) = uvar(i,5)
1025 ENDDO
1026 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
max
#define max(a, b)
Definition macros.h:21
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72