51
52
53
54 USE matparam_def_mod
56
57
58
59#include "implicit_f.inc"
60#include "comlock.inc"
61
62
63
64#include "mvsiz_p.inc"
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102#include "param_c.inc"
103#include "sms_c.inc"
104#include "impl1_c.inc"
105
106 INTEGER ,INTENT(IN) :: NVARTMP
107 INTEGER ,INTENT(INOUT) :: VARTMP(NEL,NVARTMP)
108 INTEGER, INTENT(IN) :: ISMSTR, IHET, JSMS
109 INTEGER NEL, NUPARAM, NUVAR,
110 . NGL(NEL),MAT(NEL),IPM(NPROPMI,*)
112 . time,timestep,uparam(*),
113 . rho(nel),rho0(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 . sig0xx(nel),sig0yy(nel),sig0zz(nel),
119 . sig0xy(nel),sig0yz(nel),sig0zx(nel),
120 . epsp(nel)
121 TYPE(MATPARAM_STRUCT_) ,INTENT(IN) :: MATPARAM
122
123
124
126 . signxx(nel),signyy(nel),signzz(nel),
127 . signxy(nel),signyz(nel),signzx(nel),
128 . soundsp(nel),viscmax(nel),et(nel),
129 . rhosp(mvsiz)
130
131
132
134 . uvar(nel,nuvar), off(nel), pla(nel), offg(nel), rhoref(mvsiz)
135
136
137
138 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
140 . tf(*),finter
141 EXTERNAL finter
142
143
144
145
146
147
148
149
150
151
152
153 INTEGER :: I,II,J,J1,J2,IADBUF,NFUNC,
154 . NINDEX_EPSSM,NINDEX_EPST,NINDEX_MU,NLOAD,NUNLOAD,MX,
155 . ILOAD0,IDAMAGE,ITENS,NUMTABL
156 INTEGER, DIMENSION(100) :: IFUNC
157 INTEGER, DIMENSION(NEL) :: JJ,IE_CST,ILOAD
158 INTEGER, DIMENSION(NEL) :: INDEX_EPSSM,INDEX_EPST,INDEX_MU
159 INTEGER, DIMENSION(NEL) :: IPOS1,IADP1,ILENP1
160 INTEGER, DIMENSION(NEL,2) :: IPOS
161 my_real :: yld_emax,r,fac,yp1,yp2,epsp0,epspmax,e,delta,svm2,svm
162 . aa,e0,nu,emax,epssmax,dsig,p,expo,hys,de,yfacc,mu,
alpha
163 my_real,
DIMENSION(NEL) :: alpha_1,aa1,aa2,epst,yldmin,yldmax,yldelas,
164 . g,sig0,eps0,epss,tab_loc1,tab_loc2,yld,dfc
165 my_real,
DIMENSION(NEL,2) :: xvec
166
167 mx = mat(1)
168 nfunc = ipm(10,mx)
169 DO j=1,nfunc
170 ifunc(j)=ipm(10+j,mx)
171 ENDDO
172
173 iadbuf = ipm(7,mx)-1
174 e0 = uparam(iadbuf + 2)
175 aa = uparam(iadbuf + 3)
176 epssmax = uparam(iadbuf + 4)
177 g(:) = uparam(iadbuf + 5)
178 nu = uparam(iadbuf + 6)
179 epsp0 = uparam(iadbuf + 9)
180 nload = nint(uparam(iadbuf + 7))
181 nunload = nint(uparam(iadbuf + 8))
182 idamage = nint(uparam(iadbuf + 2*nfunc + 9))
183 itens = nint(uparam(iadbuf + 2*nfunc + 13))
184 yfacc = uparam(iadbuf + 2*nfunc + 8)
185 expo = uparam(iadbuf + 2*nfunc + 10)
186 hys = uparam(iadbuf + 2*nfunc + 11)
187 emax = uparam(iadbuf + 2*nfunc + 12)
188 yld_emax= uparam(iadbuf + 2*nfunc + 15)
189 numtabl = matparam%NTABLE
190
191
192
193 DO i=1,nel
194 epst(i) = epsxx(i)**2 + epsyy(i)**2 + epszz(i)**2
195 . + half*(epsxy(i)**2+epsyz(i)**2+epszx(i)**2)
196 epst(i) = sqrt(epst(i))
197 eps0(i) = uvar(i,1)
198 sig0(i) = uvar(i,2)
199 ENDDO
200
201
202
203 nindex_epst = 0
204 nindex_epssm = 0
205
206
207 DO i=1,nel
208 ipos(i,1) = vartmp(i,1)
209 ipos(i,2) = 1
210 xvec(i,1) = epst(i)
211 xvec(i,2) = epsp0
212 ENDDO
213
215 vartmp(1:nel,1)=ipos(1:nel,1)
216
217 DO i=1,nel
218 IF (epst(i) >= epssmax) THEN
219 nindex_epssm = nindex_epssm + 1
220 index_epssm(nindex_epssm) = i
221 yldelas(i) = yld_emax + emax*(epst(i) - epssmax)
222 ELSE
223 nindex_epst = nindex_epst + 1
224 index_epst(nindex_epst) = i
225 yldelas(i) = yld(i)
226 ENDIF
227 ENDDO
228
229
230
231 IF (nload > 1) THEN
232 epspmax = matparam%TABLE(1)%X(2)%VALUES(nload)
233#include "vectorize.inc"
234 DO i=1,nel
235 ipos(i,1) = vartmp(i,2)
236 ipos(i,2) = vartmp(i,5)
237 xvec(i,1) =
min(epst(i), epssmax)
238 xvec(i,2) = epsp(i)
239
240
241 ENDDO
242
244
245 ELSE IF (nload == 1) THEN
246
247 ipos(:,1) = vartmp(:,2)
248 ipos(:,2) = vartmp(:,5)
249 DO i=1,nel
250 xvec(i,1) =
min(epst(i), epssmax)
251 ENDDO
252
254
255 END IF
256
257 vartmp(1:nel,2) = ipos(1:nel,1)
258 vartmp(1:nel,5) = ipos(1:nel,2)
259
260 yldmax(:) = yld(:)
261
262
263
264 IF (nindex_epssm > 0) THEN
265#include "vectorize.inc"
266 DO ii=1,nindex_epssm
267 i = index_epssm(ii)
268 yldmax(i) = yldmax(i) + emax*(epst(i) - epssmax)
269 ENDDO
270 END IF
271
272
273
274 IF (nunload > 0) THEN
275 IF (nunload > 1) THEN
276 epspmax = matparam%TABLE(1)%X(2)%VALUES(nunload)
277#include "vectorize.inc"
278 DO i=1,nel
279 ipos(i,1) = vartmp(i,3)
280 ipos(i,2) = vartmp(i,6)
281 xvec(i,1) =
min(epst(i), epssmax)
282 xvec(i,2) = epsp(i)
283 ENDDO
284
286
287 ELSE IF (nunload == 1) THEN
288 ipos(:,1) = vartmp(:,3)
289 ipos(:,2) = vartmp(:,6)
290 DO i=1,nel
291 xvec(i,1) = epst(i)
292 ENDDO
293
295
296 END IF
297
298 vartmp(1:nel,3)=ipos(1:nel,1)
299 vartmp(1:nel,6)=ipos(1:nel,2)
300
301 yldmin(:) = yld(:)
302
303
304
305 IF (nindex_epssm > 0) THEN
306#include "vectorize.inc"
307 DO ii=1,nindex_epssm
308 i = index_epssm(ii)
309 yldmin(i) = yldmin(i) + emax*(epst(i) - epssmax)
310 ENDDO
311 END IF
312
313 END IF
314
315
316
317 IF (itens > 0) THEN
318 DO i=1,nel
320 sig0xx(i) =
alpha*sig0xx(i)
321 sig0yy(i) =
alpha*sig0yy(i)
322 sig0zz(i) =
alpha*sig0zz(i)
323 sig0xy(i) =
alpha*sig0xy(i)
324 sig0zx(i) =
alpha*sig0zx(i)
325 sig0yz(i) =
alpha*sig0yz(i)
326 ENDDO
327 ENDIF
328
329 DO i = 1,nel
330 iload0 = int(uvar(i,5))
331 ie_cst(i)= 0
332 delta = epst(i) - uvar(i,4)
333 e = uvar(i,3)
334
335 IF (delta >= zero ) THEN
336 yld(i) = yldmax(i)
337 iload(i) = 1
338 ELSEIF (delta < zero) THEN
339 yld(i) = yldmin(i)
340 iload(i) = -1
341 IF (idamage /= 0 ) THEN
342 yld(i) = yldmax(i)
343 ENDIF
344 ENDIF
345
346 epss(i) = epst(i) - yld(i)/ e
347 epss(i) =
max(zero, epss(i))
348 de = aa*(epss(i) - uvar(i,1))
349 IF (iload(i) == 1) THEN
350 e = e +
max(de, zero)
351 IF(iload0 == -1) e= uvar(i,3)
352 uvar(i,1) =
max(uvar(i,1), epss(i))
353 ELSE
354 e = e +
min(de ,zero)
355 IF(iload0 == 1) e= uvar(i,3)
356 uvar(i,1) =
min(epss(i),uvar(i,1))
357 ENDIF
360 uvar(i,3) = e
361
362 aa1(i) = e*(one-nu)/(one + nu)/(one - two*nu)
363 aa2(i) = aa1(i)*nu/(one - nu)
364 g(i) = half*e/(one + nu)
365
366 signxx(i)= aa1(i)*depsxx(i) + aa2(i)*(depsyy(i) + depszz(i))
367 signyy(i)= aa1(i)*depsyy(i) + aa2(i)*(depsxx(i) + depszz(i))
368 signzz(i)= aa1(i)*depszz(i) + aa2(i)*(depsxx(i) + depsyy(i))
369 signxy(i)= g(i) *depsxy(i)
370 signyz(i)= g(i) *depsyz(i)
371 signzx(i)= g(i) *depszx(i)
372 dsig = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
373 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
374 dsig =sqrt(dsig)
375
376 signxx(i)=sig0xx(i) + aa1(i)*depsxx(i)
377 . + aa2(i)*(depsyy(i) + depszz(i))
378 signyy(i)=sig0yy(i) + aa1(i)*depsyy(i)
379 . + aa2(i)*(depsxx(i) + depszz(i))
380 signzz(i)=sig0zz(i) + aa1(i)*depszz(i)
381 . + aa2(i)*(depsxx(i) + depsyy(i))
382 signxy(i)=sig0xy(i) + g(i) *depsxy(i)
383 signyz(i)=sig0yz(i) + g(i) *depsyz(i)
384 signzx(i)=sig0zx(i) + g(i) *depszx(i)
385
386 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
387 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
388 svm =sqrt(svm2)
389
390 IF (idamage == 0 ) THEN
391 IF(delta >= zero ) THEN
392 yld(i) =
min(yldmax(i),svm)
393 ELSEIF(delta < zero) THEN
394 yld(i) =
max(yldmin(i),svm)
395 yld(i) =
min(yld(i),uvar(i,6))
396 ie_cst(i) = 1
397 IF (dsig > yldmin(i) .AND. dsig > svm)THEN
398 yld(i) = yldmin(i)
399 ie_cst(i) = 0
400 ENDIF
401 ENDIF
402 ELSE
403 yld(i) = yldmax(i)
404 uvar(i,8) = uvar(i,8) + half*(yldelas(i) + uvar(i,9))*delta
405 uvar(i,8) =
max(zero, uvar(i,8))
406 uvar(i,2) =
max(uvar(i,2) , uvar(i,8))
407 ENDIF
408 uvar(i,9) = yldelas(i)
409 ENDDO
410
411
412
413
414 IF (idtmins/=2.OR.jsms==0)THEN
415 DO i = 1,nel
416
417
418
419
420 soundsp(i) = sqrt(aa1(i)/rho0(i))
421 rhosp(i) = rho0(i)
422 viscmax(i) = zero
423 ENDDO
424 ELSEIF (ismstr==1 .OR. ismstr==11)THEN
425 DO i = 1,nel
426
427
428
429
430 soundsp(i) = sqrt(aa1(i)/rhoref(i))
431 rhosp(i) = rhoref(i)
432 viscmax(i) = zero
433 ENDDO
434 ELSE
435 DO i = 1,nel
436 IF (abs(offg(i)) <= one)THEN
437
438
439
440
441 soundsp(i) = sqrt(aa1(i)/rho(i))
442 rhosp(i) = rho(i)
443 ELSE
444
445
446
447
448 soundsp(i) = sqrt(aa1(i)/rhoref(i))
449 rhosp(i) = rhoref(i)
450 END IF
451 viscmax(i) = zero
452 ENDDO
453 END IF
454
455
456
457
458 IF (idamage == 0 ) THEN
459 DO i=1,nel
460 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
461 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
462 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
463 signxy(i)= g(i) *epsxy(i)
464 signyz(i)= g(i) *epsyz(i)
465 signzx(i)= g(i) *epszx(i)
466
467 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
468 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
469 svm =sqrt(svm2)
470 r = yld(i)/
max(em20,svm)
471 et(i) = r
472 signxx(i)=signxx(i)*r
473 signyy(i)=signyy(i)*r
474 signzz(i)=signzz(i)*r
475 signxy(i)=signxy(i)*r
476 signyz(i)=signyz(i)*r
477 signzx(i)=signzx(i)*r
478
479
480
481
482
483
484
485 IF (ie_cst(i) == 1) THEN
486 iload0 = int(uvar(i,5))
487 IF(iload0 /= iload(i))THEN
488 iload(i) = iload0
489 uvar(i,1) = eps0(i)
490 ENDIF
491 ENDIF
492 uvar(i,4) = epst(i)
493 uvar(i,2) = svm*r
494 uvar(i,5) = iload(i)
495 uvar(i,6) = yld(i)
496 uvar(i,7) = epsp(i)
497 ENDDO
498
499 ELSEIF (idamage == 1) THEN
500
501 DO i=1,nel
502 r = one
503 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
504 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
505 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i
506 signxy(i)= g(i) *epsxy(i)
507 signyz(i)= g(i) *epsyz(i)
508 signzx(i)= g(i) *epszx(i)
509
510 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
511 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
512 svm =sqrt(svm2)
513 r = (yld(i)/
max(em20,svm))
514 et(i) = r
515 signxx(i)=signxx(i)*r
516 signyy(i)=signyy(i)*r
517 signzz(i)=signzz(i)*r
518 signxy(i)=signxy(i)*r
519 signyz(i)=signyz(i)*r
520 signzx(i)=signzx(i)*r
521
522 IF (iload(i) == -1) THEN
523 r = (yldmin(i)/
max(em20,yldelas(i)))
524 p = third*(signxx(i) + signyy(i) + signzz(i))
525 signxx(i)=(signxx(i) - p)*r + p
526 signyy(i)=(signyy(i) - p)*r + p
527 signzz(i)=(signzz(i) - p)*r + p
528 signxy(i)=signxy(i)*r
529 signyz(i)=signyz(i)*r
530 signzx(i)=signzx(i)*r
531 ENDIF
532
533 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
534 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
535 svm =sqrt(svm2)
536 uvar(i,4) = epst(i)
537 uvar(i,5) = iload(i)
538 uvar(i,6) = yld(i)
539 uvar(i,7) = epsp(i)
540 ENDDO
541
542 ELSEIF(idamage == 2 ) THEN
543
544 DO i=1,nel
545 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
546 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
547 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
548 signxy(i)= g(i) *epsxy(i)
549 signyz(i)= g(i) *epsyz(i)
550 signzx(i)= g(i) *epszx(i)
551
552 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
553 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
554 svm =sqrt(svm2)
555 r = yld(i)/
max(em20,svm)
556 et(i) = r
557 signxx
558 signyy(i)=signyy(i)*r
559 signzz(i)=signzz(i)*r
560 signxy(i)=signxy(i)*r
561 signyz(i)=signyz(i)*r
562 signzx(i)=signzx(i)*r
563
564 IF (iload(i) == -1) THEN
565 r = yldmin(i)/
max(em20,yldelas(i))
566 signxx(i)=signxx(i)*r
567 signyy(i)=signyy(i)*r
568 signzz(i)=signzz(i)*r
569 signxy(i)=signxy(i)*r
570 signyz(i)=signyz(i)*r
571 signzx(i)=signzx(i)*r
572 et(i) = r*et(i)
573 ENDIF
574
575 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
576 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
577 svm =sqrt(svm2)
578 uvar(i,4) = epst(i)
579 uvar(i,5) = iload(i)
580 uvar(i,6) = yld(i)
581 uvar(i,7) = epsp(i)
582 ENDDO
583
584 ELSEIF (idamage == 3) THEN
585
586 DO i=1,nel
587 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
588 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
589 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
590 signxy(i)= g(i) *epsxy(i)
591 signyz(i)= g(i) *epsyz(i)
592 signzx(i)= g(i) *epszx(i)
593
594 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
595 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
596 svm =sqrt(svm2)
597 r = (yld(i)/
max(em20,svm))
598 signxx(i)=signxx(i)*r
599 signyy(i)=signyy(i)*r
600 signzz(i)=signzz(i)*r
601 signxy(i)=signxy(i)*r
602 signyz(i)=signyz(i)*r
603 signzx(i)=signzx(i)*r
604 et(i) = r
605
606 IF (iload(i)==-1 .AND.uvar(i,2)/=zero)THEN
607 r = one - (uvar(i,8)/uvar(i,2))**expo
608 r = one - (one - hys)*r
609 p = third*(signxx(i) + signyy(i) + signzz(i))
610 signxx(i)=(signxx(i) - p)*r + p
611 signyy(i)=(signyy(i) - p)*r + p
612 signzz(i)=(signzz(i) - p)*r + p
613 signxy(i)=signxy(i)*r
614 signyz(i)=signyz(i)*r
615 signzx(i)=signzx(i)*r
616 et(i) = r*et(i)
617 ENDIF
618
619 uvar(i,4) = epst(i)
620 uvar(i,5) = iload(i)
621 uvar(i,6) = yld(i)
622 uvar(i,7) = epsp(i)
623 ENDDO
624
625 ELSEIF (idamage == 4) THEN
626
627 DO i=1,nel
628 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
629 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
630 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
631 signxy(i)= g(i) *epsxy(i)
632 signyz(i)= g(i) *epsyz(i)
633 signzx(i)= g(i) *epszx(i)
634
635 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
636 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
637 svm =sqrt(svm2)
638 r = (yld(i)/
max(em20,svm))
639 et(i) = r
640 signxx(i)=signxx(i)*r
641 signyy(i)=signyy(i)*r
642 signzz(i)=signzz(i)*r
643 signxy(i)=signxy(i)*r
644 signyz(i)=signyz(i)*r
645 signzx(i)=signzx(i)*r
646
647 IF (iload(i)==-1 .AND.uvar(i,2)/= zero)THEN
648 r = one - (uvar(i,8)/uvar(i,2))**expo
649 r = one - (one - hys)*r
650 signxx(i)=signxx(i)*r
651 signyy(i)=signyy(i)*r
652 signzz(i)=signzz(i)*r
653 signxy(i)=signxy(i)*r
654 signyz(i)=signyz(i)*r
655 signzx(i)=signzx(i)*r
656 ENDIF
657 uvar(i,4) = epst(i)
658 uvar(i,5) = iload(i)
659 uvar(i,6) = yld(i)
660 uvar(i,7) = epsp(i)
661 ENDDO
662
663 ENDIF
664
665 IF (itens > 0 ) THEN
666 alpha_1(1:nel) = one
667 nindex_mu = 0
668
669 DO i=1,nel
670 mu = 1 - rho(i)/rho0(i)
671 IF (mu > 0) THEN
672 nindex_mu = nindex_mu + 1
673 index_mu(nindex_mu) = i
674 ipos1(nindex_mu) = vartmp(i,4)
675 iadp1(nindex_mu) = npf(ifunc(nfunc)) / 2 + 1
676 ilenp1(nindex_mu) = npf(ifunc(nfunc)+1) / 2 - iadp1(nindex_mu) - ipos1(nindex_mu)
677 tab_loc1(nindex_mu) = mu
678 ENDIF
679 ENDDO
680
681 IF ( nindex_mu > 0 ) THEN
682
683 CALL vinter2(tf,iadp1,ipos1,ilenp1,nindex_mu,tab_loc1,dfc,yld)
684
685 DO ii=1,nindex_mu
686 i = index_mu(ii)
687 alpha_1(i) = yfacc*yld(ii)
688 alpha_1(i) =
max(zero, alpha_1(i))
689 vartmp(i,4) = ipos1(ii)
690 ENDDO
691 ENDIF
692
693 DO i=1,nel
694 signxx(i) = alpha_1(i)*signxx(i)
695 signyy(i) = alpha_1(i)*signyy(i)
696 signzz(i) = alpha_1(i)*signzz(i)
697 signxy(i) = alpha_1(i)*signxy(i)
698 signzx(i) = alpha_1(i)*signzx(i)
699 signyz(i) = alpha_1(i)*signyz(i)
700 uvar(i,10) = alpha_1(i)
701 et(i) = alpha_1(i)*et(i)
702 ENDDO
703 ENDIF
704
705 IF (impl_s > 0 .OR. ihet > 1) THEN
706 DO i=1,nel
707 et(i) =
min(one , et(i))
708 ENDDO
709 ENDIF
710
711 RETURN
subroutine table_mat_vinterp(table, dimx, nel, ipos, xx, yy, dydx)
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)