1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine resultsmech_se(co,kon,ipkon,lakon,ne,v,
20     &  stx,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,
21     &  ielmat,ielorien,norien,orab,ntmat_,t0,t1,ithermal,prestr,
22     &  iprestr,eme,iperturb,fn,iout,vold,nmethod,
23     &  veold,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon,
24     &  xstateini,xstiff,xstate,npmat_,matname,mi,ielas,icmd,
25     &  ncmat_,nstate_,stiini,vini,ener,eei,enerini,istep,iinc,
26     &  springarea,reltime,calcul_fn,calcul_cauchy,nener,
27     &  ikin,ne0,thicke,emeini,pslavsurf,
28     &  pmastsurf,mortar,clearini,nea,neb,ielprop,prop,dfn,
29     &  idesvar,nodedesi,fn0,sti,icoordinate,dxstiff,
30     &  ialdesi,xdesi)
31!
32!     calculates stresses and the material tangent at the integration
33!     points and the internal forces at the nodes
34!
35      implicit none
36!
37      integer cauchy
38!
39      character*8 lakon(*),lakonl
40      character*80 amat,matname(*)
41!
42      integer kon(*),konl(26),nea,neb,mi(*),mint2d,nopes,
43     &  nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),
44     &  ielorien(mi(3),*),ntmat_,ipkon(*),ne0,iflag,null,
45     &  istep,iinc,mt,ne,mattyp,ithermal(*),iprestr,i,j,k,m1,m2,jj,
46     &  i1,m3,m4,kk,nener,indexe,nope,norien,iperturb(*),iout,
47     &  icmd,ihyper,nmethod,kode,imat,mint3d,iorien,ielas,
48     &  istiff,ncmat_,nstate_,ikin,ilayer,nlayer,ki,kl,ielprop(*),
49     &  nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,calcul_fn,
50     &  calcul_cauchy,nopered,mortar,jfaces,igauss,
51     &  idesvar,node,nodedesi(*),kscale,idir,nlgeom_undo,
52     &  iactive,icoordinate,ialdesi(*),ii,
53     &  node1,node2,ifaceqexp(2,20),ifacewexp(2,15)
54!
55      real*8 co(3,*),v(0:mi(2),*),shp(4,26),stiini(6,mi(1),*),
56     &  stx(6,mi(1),*),xl(3,26),vl(0:mi(2),26),stre(6),prop(*),
57     &  elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),xs2(3,7),
58     &  alcon(0:6,ntmat_,*),vini(0:mi(2),*),thickness,
59     &  alzero(*),orab(7,*),elas(21),rho,fn(0:mi(2),*),
60     &  fnl(3,10),skl(3,3),beta(6),xl2(3,8),qa(4),
61     &  vkl(0:3,3),t0(*),t1(*),prestr(6,mi(1),*),eme(6,mi(1),*),
62     &  ckl(3,3),vold(0:mi(2),*),eloc(9),veold(0:mi(2),*),
63     &  springarea(2,*),elconloc(21),eth(6),xkl(3,3),voldl(0:mi(2),26),
64     &  xikl(3,3),ener(mi(1),*),emec(6),eei(6,mi(1),*),enerini(mi(1),*),
65     &  emec0(6),veoldl(0:mi(2),26),xsj2(3),shp2(7,8),
66     &  e,un,al,um,am1,xi,et,ze,tt,exx,eyy,ezz,exy,exz,eyz,
67     &  xsj,vj,t0l,t1l,dtime,weight,pgauss(3),vij,time,ttime,
68     &  plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
69     &  xstiff(27,mi(1),*),xstate(nstate_,mi(1),*),plconloc(802),
70     &  vokl(3,3),xstateini(nstate_,mi(1),*),vikl(3,3),
71     &  gs(8,4),a,reltime,tlayer(4),dlayer(4),xlayer(mi(3),4),
72     &  thicke(mi(3),*),emeini(6,mi(1),*),clearini(3,9,*),
73     &  pslavsurf(3,*),pmastsurf(6,*),sti(6,mi(1),*),
74     &  fn0(0:mi(2),*),dfn(0:mi(2),*),hglf(3,4),ahr,
75     &  dxstiff(27,mi(1),ne,*),xdesi(3,*)
76!
77!
78!
79      include "gauss.f"
80!
81      iflag=3
82      null=0
83      qa(3)=-1.d0
84      qa(4)=0.d0
85!
86      mt=mi(2)+1
87!
88!     nodes in expansion direction of hex element
89!
90      data ifaceqexp /5,17,
91     &                6,18,
92     &                7,19,
93     &                8,20,
94     &                1,17,
95     &                2,18,
96     &                3,19,
97     &                4,20,
98     &                13,0,
99     &                14,0,
100     &                15,0,
101     &                16,0,
102     &                9,0,
103     &                10,0,
104     &                11,0,
105     &                12,0,
106     &                0,0,
107     &                0,0,
108     &                0,0,
109     &                0,0/
110!
111!     nodes in expansion direction of wedge element
112!
113      data ifacewexp /4,13,
114     &                5,14,
115     &                6,16,
116     &                1,13,
117     &                2,14,
118     &                3,15,
119     &                10,0,
120     &                11,0,
121     &                12,0,
122     &                7,0,
123     &                8,0,
124     &                9,0,
125     &                0,0,
126     &                0,0,
127     &                0,0/
128!
129!     -------------------------------------------------------------
130!     Initialisation of the loop for the variation of
131!     the internal forces
132!     -------------------------------------------------------------
133!
134!
135!     Loop over all elements in thread
136      do ii=nea,neb
137         if(idesvar.gt.0) then
138            i=ialdesi(ii)
139         else
140            i=ii
141         endif
142!
143         lakonl=lakon(i)
144!
145         if(ipkon(i).lt.0) cycle
146!
147!        no 3D-fluid elements
148!
149         if(lakonl(1:1).eq.'F') cycle
150         if(lakonl(1:7).eq.'DCOUP3D') cycle
151!
152         if(lakonl(7:8).ne.'LC') then
153!
154            imat=ielmat(1,i)
155            amat=matname(imat)
156            if(norien.gt.0) then
157               iorien=max(0,ielorien(1,i))
158            else
159               iorien=0
160            endif
161!
162            if(nelcon(1,imat).lt.0) then
163               ihyper=1
164            else
165               ihyper=0
166            endif
167         elseif(lakonl(4:5).eq.'20') then
168!
169!     composite materials
170!
171            mint2d=4
172            nopes=8
173!     determining the number of layers
174!
175            nlayer=0
176            do k=1,mi(3)
177               if(ielmat(k,i).ne.0) then
178                  nlayer=nlayer+1
179               endif
180            enddo
181!
182!     determining the layer thickness and global thickness
183!     at the shell integration points
184!
185            iflag=1
186            indexe=ipkon(i)
187            do kk=1,mint2d
188               xi=gauss3d2(1,kk)
189               et=gauss3d2(2,kk)
190               call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
191               tlayer(kk)=0.d0
192               do k=1,nlayer
193                  thickness=0.d0
194                  do j=1,nopes
195                     thickness=thickness+thicke(k,indexe+j)*shp2(4,j)
196                  enddo
197                  tlayer(kk)=tlayer(kk)+thickness
198                  xlayer(k,kk)=thickness
199               enddo
200            enddo
201            iflag=3
202!
203            ilayer=0
204            do k=1,4
205               dlayer(k)=0.d0
206            enddo
207         elseif(lakonl(4:5).eq.'15') then
208!
209!     composite materials
210!
211            mint2d=3
212            nopes=6
213!     determining the number of layers
214!
215            nlayer=0
216            do k=1,mi(3)
217               if(ielmat(k,i).ne.0) then
218                  nlayer=nlayer+1
219               endif
220            enddo
221!
222!     determining the layer thickness and global thickness
223!     at the shell integration points
224!
225            iflag=1
226            indexe=ipkon(i)
227            do kk=1,mint2d
228               xi=gauss3d10(1,kk)
229               et=gauss3d10(2,kk)
230               call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
231               tlayer(kk)=0.d0
232               do k=1,nlayer
233                  thickness=0.d0
234                  do j=1,nopes
235                     thickness=thickness+thicke(k,indexe+j)*shp2(4,j)
236                  enddo
237                  tlayer(kk)=tlayer(kk)+thickness
238                  xlayer(k,kk)=thickness
239               enddo
240            enddo
241            iflag=3
242!
243            ilayer=0
244            do k=1,3
245               dlayer(k)=0.d0
246            enddo
247!
248         endif
249!
250         indexe=ipkon(i)
251c     Bernhardi start
252         if(lakonl(1:5).eq.'C3D8I') then
253            nope=11
254         elseif(lakonl(4:5).eq.'20') then
255c     Bernhardi end
256            nope=20
257         elseif(lakonl(4:4).eq.'8') then
258            nope=8
259         elseif(lakonl(4:5).eq.'10') then
260            nope=10
261         elseif(lakonl(4:4).eq.'4') then
262            nope=4
263         elseif(lakonl(4:5).eq.'15') then
264            nope=15
265         elseif(lakonl(4:4).eq.'6') then
266            nope=6
267         elseif((lakonl(1:1).eq.'E').and.(lakonl(7:7).ne.'F')) then
268!
269!           spring elements, contact spring elements and
270!           dashpot elements
271!
272            if(lakonl(7:7).eq.'C') then
273!
274!              contact spring elements
275!
276               if(mortar.eq.1) then
277!
278!                 face-to-face penalty
279!
280                  nope=kon(ipkon(i))
281               elseif(mortar.eq.0) then
282!
283!                 node-to-face penalty
284!
285                  nope=ichar(lakonl(8:8))-47
286                  konl(nope+1)=kon(indexe+nope+1)
287               endif
288            else
289!
290!              genuine spring elements and dashpot elements
291!
292               nope=ichar(lakonl(8:8))-47
293            endif
294         else
295            cycle
296         endif
297!
298!        computation of the coordinates of the local nodes w/ variation
299!        if the designnode belongs to the considered element
300!
301         if(idesvar.gt.0) then
302!
303            if(icoordinate.eq.1) then
304!
305!              the coordinates of the nodes are the design variables
306!
307               do j=1,nope
308                  node=kon(indexe+j)
309                  if(node.eq.nodedesi(idesvar)) then
310                     iactive=j
311                     exit
312                  endif
313               enddo
314            endif
315!
316         endif
317!
318         if(lakonl(4:5).eq.'8R') then
319            mint3d=1
320         elseif(lakonl(4:7).eq.'20RB') then
321            if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
322               mint3d=50
323            else
324               call beamintscheme(lakonl,mint3d,ielprop(i),prop,
325     &              null,xi,et,ze,weight)
326            endif
327         elseif((lakonl(4:4).eq.'8').or.
328     &          (lakonl(4:6).eq.'20R')) then
329            if(lakonl(7:8).eq.'LC') then
330               mint3d=8*nlayer
331            else
332               mint3d=8
333            endif
334         elseif(lakonl(4:4).eq.'2') then
335            mint3d=27
336         elseif(lakonl(4:5).eq.'10') then
337            mint3d=4
338         elseif(lakonl(4:4).eq.'4') then
339            mint3d=1
340         elseif(lakonl(4:5).eq.'15') then
341            if(lakonl(7:8).eq.'LC') then
342               mint3d=6*nlayer
343            else
344               mint3d=9
345            endif
346         elseif(lakonl(4:4).eq.'6') then
347            mint3d=2
348         elseif(lakonl(1:1).eq.'E') then
349            mint3d=0
350         endif
351
352         do j=1,nope
353            konl(j)=kon(indexe+j)
354            do k=1,3
355               xl(k,j)=co(k,konl(j))
356               vl(k,j)=v(k,konl(j))
357               voldl(k,j)=vold(k,konl(j))
358            enddo
359         enddo
360!
361!        computation of the coordinates of the local nodes w/ variation
362!        if the designnode belongs to the considered element
363!
364         if((idesvar.gt.0).and.(icoordinate.eq.1)) then
365            do j=1,3
366               xl(j,iactive)=xl(j,iactive)+xdesi(j,idesvar)
367            enddo
368            if(lakonl(1:5).eq.'C3D20') then
369               if((lakonl(7:7).eq.'A').or.
370     &            (lakonl(7:7).eq.'L').or.
371     &            (lakonl(7:7).eq.'S').or.
372     &            (lakonl(7:7).eq.'E')) then
373                  node1=ifaceqexp(1,iactive)
374                  node2=ifaceqexp(2,iactive)
375                  do j=1,3
376                     if(iactive.le.8) then
377                        xl(j,node1)=xl(j,node1)+xdesi(j,idesvar)
378                        xl(j,node2)=xl(j,node2)+xdesi(j,idesvar)
379                     else
380                        xl(j,node1)=xl(j,node1)+xdesi(j,idesvar)
381                     endif
382                  enddo
383               endif
384            elseif(lakonl(1:5).eq.'C3D15') then
385               if((lakonl(7:7).eq.'A').or.
386     &            (lakonl(7:7).eq.'L').or.
387     &            (lakonl(7:7).eq.'S').or.
388     &            (lakonl(7:7).eq.'E')) then
389                  node1=ifacewexp(1,iactive)
390                  node2=ifacewexp(2,iactive)
391                  do j=1,3
392                     if(iactive.le.6) then
393                        xl(j,node1)=xl(j,node1)+xdesi(j,idesvar)
394                        xl(j,node2)=xl(j,node2)+xdesi(j,idesvar)
395                     else
396                        xl(j,node1)=xl(j,node1)+xdesi(j,idesvar)
397                     endif
398                  enddo
399               endif
400            elseif(lakonl(1:4).eq.'C3D8') then
401               if((lakonl(7:7).eq.'A').or.
402     &            (lakonl(7:7).eq.'L').or.
403     &            (lakonl(7:7).eq.'S').or.
404     &            (lakonl(7:7).eq.'E')) then
405                  node1=ifaceqexp(1,iactive)
406                  do j=1,3
407                     xl(j,node1)=xl(j,node1)+xdesi(j,idesvar)
408                  enddo
409               endif
410            elseif(lakonl(1:4).eq.'C3D6') then
411               if((lakonl(7:7).eq.'A').or.
412     &              (lakonl(7:7).eq.'L').or.
413     &              (lakonl(7:7).eq.'S').or.
414     &              (lakonl(7:7).eq.'E')) then
415                  node1=ifaceqexp(1,iactive)
416                  do j=1,3
417                     xl(j,node1)=xl(j,node1)+xdesi(j,idesvar)
418                  enddo
419               endif
420            endif
421         endif
422!
423!        calculating the forces for the contact elements
424!
425         if(mint3d.eq.0) then
426!
427!           "normal" spring and dashpot elements
428!
429            kode=nelcon(1,imat)
430            if(lakonl(7:7).eq.'A') then
431               t0l=0.d0
432               t1l=0.d0
433               if(ithermal(1).eq.1) then
434                  t0l=(t0(konl(1))+t0(konl(2)))/2.d0
435                  t1l=(t1(konl(1))+t1(konl(2)))/2.d0
436               elseif(ithermal(1).ge.2) then
437                  t0l=(t0(konl(1))+t0(konl(2)))/2.d0
438                  t1l=(vold(0,konl(1))+vold(0,konl(2)))/2.d0
439               endif
440            endif
441!
442!           spring elements (including contact springs)
443!
444            if(lakonl(2:2).eq.'S') then
445               if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'1').or.
446     &            (lakonl(7:7).eq.'2').or.((mortar.eq.0).and.
447     &          ((nmethod.ne.1).or.(iperturb(1).ge.2).or.(iout.ne.-1))))
448     &               then
449                  call springforc_n2f(xl,konl,vl,imat,elcon,nelcon,elas,
450     &              fnl,ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc,
451     &              plicon,nplicon,npmat_,ener(1,i),nener,
452     &              stx(1,1,i),mi,springarea(1,konl(nope+1)),nmethod,
453     &              ne0,nstate_,xstateini,xstate,reltime,
454     &              ielas,ener(1,i+ne),ielorien,orab,norien,i)
455               elseif((mortar.eq.1).and.
456     &           ((nmethod.ne.1).or.(iperturb(1).ge.2).or.(iout.ne.-1)))
457     &               then
458                  jfaces=kon(indexe+nope+2)
459                  igauss=kon(indexe+nope+1)
460                  call springforc_f2f(xl,vl,imat,elcon,nelcon,elas,
461     &              fnl,ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc,
462     &              plicon,nplicon,npmat_,ener(1,i),nener,
463     &              stx(1,1,i),mi,springarea(1,igauss),nmethod,
464     &              ne0,nstate_,xstateini,xstate,reltime,
465     &              ielas,jfaces,igauss,pslavsurf,pmastsurf,
466     &              clearini,ener(1,i+ne),kscale,konl,iout,i)
467               endif
468!
469!              next lines are not executed in linstatic.c before the
470!              setup of the stiffness matrix (i.e. nmethod=1 and
471!              iperturb(1)<1 and iout=-1).
472!
473               if((lakonl(7:7).eq.'A').or.
474     &           ((nmethod.ne.1).or.(iperturb(1).ge.2).or.(iout.ne.-1)))
475     &              then
476                  if(idesvar.eq.0) then
477                     do j=1,nope
478                        do k=1,3
479                           fn0(k,indexe+j)=fn0(k,indexe+j)+fnl(k,j)
480                        enddo
481                     enddo
482                  else
483                     do j=1,nope
484                        do k=1,3
485                           dfn(k,konl(j))=dfn(k,konl(j))+fnl(k,j)
486                        enddo
487                     enddo
488                  endif
489               endif
490!
491!              dashpot elements (including contact dashpots)
492!
493            elseif((nmethod.eq.4).or.(nmethod.eq.5).or.
494     &             ((abs(nmethod).eq.1).and.(iperturb(1).ge.2))) then
495            endif
496         elseif(ikin.eq.1) then
497            do j=1,nope
498               do k=1,3
499                  veoldl(k,j)=veold(k,konl(j))
500               enddo
501            enddo
502         endif
503!
504         do jj=1,mint3d
505            if(lakonl(4:5).eq.'8R') then
506               xi=gauss3d1(1,jj)
507               et=gauss3d1(2,jj)
508               ze=gauss3d1(3,jj)
509               weight=weight3d1(jj)
510            elseif(lakonl(4:7).eq.'20RB') then
511               if((lakonl(8:8).eq.'R').or.(lakonl(8:8).eq.'C')) then
512                  xi=gauss3d13(1,jj)
513                  et=gauss3d13(2,jj)
514                  ze=gauss3d13(3,jj)
515                  weight=weight3d13(jj)
516               else
517                  call beamintscheme(lakonl,mint3d,ielprop(i),prop,
518     &                 jj,xi,et,ze,weight)
519               endif
520            elseif((lakonl(4:4).eq.'8').or.
521     &             (lakonl(4:6).eq.'20R'))
522     &        then
523               if(lakonl(7:8).ne.'LC') then
524                  xi=gauss3d2(1,jj)
525                  et=gauss3d2(2,jj)
526                  ze=gauss3d2(3,jj)
527                  weight=weight3d2(jj)
528               else
529                  kl=mod(jj,8)
530                  if(kl.eq.0) kl=8
531!
532                  xi=gauss3d2(1,kl)
533                  et=gauss3d2(2,kl)
534                  ze=gauss3d2(3,kl)
535                  weight=weight3d2(kl)
536!
537                  ki=mod(jj,4)
538                  if(ki.eq.0) ki=4
539!
540                  if(kl.eq.1) then
541                     ilayer=ilayer+1
542                     if(ilayer.gt.1) then
543                        do k=1,4
544                           dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
545                        enddo
546                     endif
547                  endif
548                  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
549     &                 tlayer(ki)-1.d0
550                  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
551!
552!                 material and orientation
553!
554                  imat=ielmat(ilayer,i)
555                  amat=matname(imat)
556                  if(norien.gt.0) then
557                     iorien=max(0,ielorien(ilayer,i))
558                  else
559                     iorien=0
560                  endif
561!
562                  if(nelcon(1,imat).lt.0) then
563                     ihyper=1
564                  else
565                     ihyper=0
566                  endif
567               endif
568            elseif(lakonl(4:4).eq.'2') then
569               xi=gauss3d3(1,jj)
570               et=gauss3d3(2,jj)
571               ze=gauss3d3(3,jj)
572               weight=weight3d3(jj)
573            elseif(lakonl(4:5).eq.'10') then
574               xi=gauss3d5(1,jj)
575               et=gauss3d5(2,jj)
576               ze=gauss3d5(3,jj)
577               weight=weight3d5(jj)
578            elseif(lakonl(4:4).eq.'4') then
579               xi=gauss3d4(1,jj)
580               et=gauss3d4(2,jj)
581               ze=gauss3d4(3,jj)
582               weight=weight3d4(jj)
583            elseif(lakonl(4:5).eq.'15') then
584               if(lakonl(7:8).ne.'LC') then
585                  xi=gauss3d8(1,jj)
586                  et=gauss3d8(2,jj)
587                  ze=gauss3d8(3,jj)
588                  weight=weight3d8(jj)
589               else
590                  kl=mod(jj,6)
591                  if(kl.eq.0) kl=6
592!
593                  xi=gauss3d10(1,kl)
594                  et=gauss3d10(2,kl)
595                  ze=gauss3d10(3,kl)
596                  weight=weight3d10(kl)
597!
598                  ki=mod(jj,3)
599                  if(ki.eq.0) ki=3
600!
601                  if(kl.eq.1) then
602                     ilayer=ilayer+1
603                     if(ilayer.gt.1) then
604                        do k=1,3
605                           dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
606                        enddo
607                     endif
608                  endif
609                  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
610     &                 tlayer(ki)-1.d0
611                  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
612!
613!                 material and orientation
614!
615                  imat=ielmat(ilayer,i)
616                  amat=matname(imat)
617                  if(norien.gt.0) then
618                     iorien=max(0,ielorien(ilayer,i))
619                  else
620                     iorien=0
621                  endif
622!
623                  if(nelcon(1,imat).lt.0) then
624                     ihyper=1
625                  else
626                     ihyper=0
627                  endif
628               endif
629            elseif(lakonl(4:4).eq.'6') then
630               xi=gauss3d7(1,jj)
631               et=gauss3d7(2,jj)
632               ze=gauss3d7(3,jj)
633               weight=weight3d7(jj)
634            endif
635!
636c     Bernhardi start
637            if(lakonl(1:5).eq.'C3D8R') then
638               call shape8hr(xl,xsj,shp,gs,a)
639            elseif(lakonl(1:5).eq.'C3D8I') then
640               call shape8hu(xi,et,ze,xl,xsj,shp,iflag)
641            elseif(nope.eq.20) then
642c     Bernhardi end
643               if(lakonl(7:7).eq.'A') then
644                  call shape20h_ax(xi,et,ze,xl,xsj,shp,iflag)
645               elseif((lakonl(7:7).eq.'E').or.
646     &                (lakonl(7:7).eq.'S')) then
647                  call shape20h_pl(xi,et,ze,xl,xsj,shp,iflag)
648               else
649                  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
650               endif
651            elseif(nope.eq.8) then
652               call shape8h(xi,et,ze,xl,xsj,shp,iflag)
653            elseif(nope.eq.10) then
654               call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
655            elseif(nope.eq.4) then
656               call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
657            elseif(nope.eq.15) then
658               call shape15w(xi,et,ze,xl,xsj,shp,iflag)
659            else
660               call shape6w(xi,et,ze,xl,xsj,shp,iflag)
661            endif
662!
663!                 vkl(m2,m3) contains the derivative of the m2-
664!                 component of the displacement with respect to
665!                 direction m3
666!
667            do m2=1,3
668               do m3=1,3
669                  vkl(m2,m3)=0.d0
670               enddo
671            enddo
672!
673            do m1=1,nope
674               do m2=1,3
675                  do m3=1,3
676                     vkl(m2,m3)=vkl(m2,m3)+shp(m3,m1)*vl(m2,m1)
677                  enddo
678c                  write(*,*) 'vnoeie',i,konl(m1),(vkl(m2,k),k=1,3)
679               enddo
680            enddo
681!
682!           for frequency analysis or buckling with preload the
683!           strains are calculated with respect to the deformed
684!           configuration
685!           for a linear iteration within a nonlinear increment:
686!           the tangent matrix is calculated at strain at the end
687!           of the previous increment
688!
689            if((iperturb(1).eq.1).or.(iperturb(1).eq.-1))then
690               do m2=1,3
691                  do m3=1,3
692                     vokl(m2,m3)=0.d0
693                  enddo
694               enddo
695!
696               do m1=1,nope
697                  do m2=1,3
698                     do m3=1,3
699                        vokl(m2,m3)=vokl(m2,m3)+
700     &                       shp(m3,m1)*voldl(m2,m1)
701                     enddo
702                  enddo
703               enddo
704            endif
705!
706            kode=nelcon(1,imat)
707!
708!           calculating the strain
709!
710!           attention! exy,exz and eyz are engineering strains!
711!
712            exx=vkl(1,1)
713            eyy=vkl(2,2)
714            ezz=vkl(3,3)
715            exy=vkl(1,2)+vkl(2,1)
716            exz=vkl(1,3)+vkl(3,1)
717            eyz=vkl(2,3)+vkl(3,2)
718!
719            if(iperturb(2).eq.1) then
720!
721!                 Lagrangian strain
722!
723               exx=exx+(vkl(1,1)**2+vkl(2,1)**2+vkl(3,1)**2)/2.d0
724               eyy=eyy+(vkl(1,2)**2+vkl(2,2)**2+vkl(3,2)**2)/2.d0
725               ezz=ezz+(vkl(1,3)**2+vkl(2,3)**2+vkl(3,3)**2)/2.d0
726               exy=exy+vkl(1,1)*vkl(1,2)+vkl(2,1)*vkl(2,2)+
727     &              vkl(3,1)*vkl(3,2)
728               exz=exz+vkl(1,1)*vkl(1,3)+vkl(2,1)*vkl(2,3)+
729     &              vkl(3,1)*vkl(3,3)
730               eyz=eyz+vkl(1,2)*vkl(1,3)+vkl(2,2)*vkl(2,3)+
731     &              vkl(3,2)*vkl(3,3)
732!
733!           for frequency analysis or buckling with preload the
734!           strains are calculated with respect to the deformed
735!           configuration
736!
737            elseif(iperturb(1).eq.1) then
738               exx=exx+vokl(1,1)*vkl(1,1)+vokl(2,1)*vkl(2,1)+
739     &              vokl(3,1)*vkl(3,1)
740               eyy=eyy+vokl(1,2)*vkl(1,2)+vokl(2,2)*vkl(2,2)+
741     &              vokl(3,2)*vkl(3,2)
742               ezz=ezz+vokl(1,3)*vkl(1,3)+vokl(2,3)*vkl(2,3)+
743     &              vokl(3,3)*vkl(3,3)
744               exy=exy+vokl(1,1)*vkl(1,2)+vokl(1,2)*vkl(1,1)+
745     &              vokl(2,1)*vkl(2,2)+vokl(2,2)*vkl(2,1)+
746     &              vokl(3,1)*vkl(3,2)+vokl(3,2)*vkl(3,1)
747               exz=exz+vokl(1,1)*vkl(1,3)+vokl(1,3)*vkl(1,1)+
748     &              vokl(2,1)*vkl(2,3)+vokl(2,3)*vkl(2,1)+
749     &              vokl(3,1)*vkl(3,3)+vokl(3,3)*vkl(3,1)
750               eyz=eyz+vokl(1,2)*vkl(1,3)+vokl(1,3)*vkl(1,2)+
751     &              vokl(2,2)*vkl(2,3)+vokl(2,3)*vkl(2,2)+
752     &              vokl(3,2)*vkl(3,3)+vokl(3,3)*vkl(3,2)
753            endif
754!
755!              storing the local strains
756!
757            if(iperturb(1).ne.-1) then
758               eloc(1)=exx
759               eloc(2)=eyy
760               eloc(3)=ezz
761               eloc(4)=exy/2.d0
762               eloc(5)=exz/2.d0
763               eloc(6)=eyz/2.d0
764            else
765!
766!              linear iteration within a nonlinear increment:
767!
768               eloc(1)=vokl(1,1)+
769     &           (vokl(1,1)**2+vokl(2,1)**2+vokl(3,1)**2)/2.d0
770               eloc(2)=vokl(2,2)+
771     &           (vokl(1,2)**2+vokl(2,2)**2+vokl(3,2)**2)/2.d0
772               eloc(3)=vokl(3,3)+
773     &           (vokl(1,3)**2+vokl(2,3)**2+vokl(3,3)**2)/2.d0
774               eloc(4)=(vokl(1,2)+vokl(2,1)+vokl(1,1)*vokl(1,2)+
775     &              vokl(2,1)*vokl(2,2)+vokl(3,1)*vokl(3,2))/2.d0
776               eloc(5)=(vokl(1,3)+vokl(3,1)+vokl(1,1)*vokl(1,3)+
777     &              vokl(2,1)*vokl(2,3)+vokl(3,1)*vokl(3,3))/2.d0
778               eloc(6)=(vokl(2,3)+vokl(3,2)+vokl(1,2)*vokl(1,3)+
779     &              vokl(2,2)*vokl(2,3)+vokl(3,2)*vokl(3,3))/2.d0
780            endif
781!
782!                 calculating the deformation gradient (needed to
783!                 convert the element stiffness matrix from spatial
784!                 coordinates to material coordinates
785!                 deformation plasticity)
786!
787            if((kode.eq.-50).or.(kode.le.-100)) then
788!
789!                    calculating the deformation gradient
790!
791c     Bernhardi start
792               xkl(1,1)=vkl(1,1)+1.0d0
793               xkl(2,2)=vkl(2,2)+1.0d0
794               xkl(3,3)=vkl(3,3)+1.0d0
795c     Bernhardi end
796               xkl(1,2)=vkl(1,2)
797               xkl(1,3)=vkl(1,3)
798               xkl(2,3)=vkl(2,3)
799               xkl(2,1)=vkl(2,1)
800               xkl(3,1)=vkl(3,1)
801               xkl(3,2)=vkl(3,2)
802!
803!                    calculating the Jacobian
804!
805               vj=xkl(1,1)*(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))
806     &              -xkl(1,2)*(xkl(2,1)*xkl(3,3)-xkl(2,3)*xkl(3,1))
807     &              +xkl(1,3)*(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))
808!
809!              inversion of the deformation gradient (only for
810!              deformation plasticity)
811!
812               if(kode.eq.-50) then
813!
814                  ckl(1,1)=(xkl(2,2)*xkl(3,3)-xkl(2,3)*xkl(3,2))/vj
815                  ckl(2,2)=(xkl(1,1)*xkl(3,3)-xkl(1,3)*xkl(3,1))/vj
816                  ckl(3,3)=(xkl(1,1)*xkl(2,2)-xkl(1,2)*xkl(2,1))/vj
817                  ckl(1,2)=(xkl(1,3)*xkl(3,2)-xkl(1,2)*xkl(3,3))/vj
818                  ckl(1,3)=(xkl(1,2)*xkl(2,3)-xkl(2,2)*xkl(1,3))/vj
819                  ckl(2,3)=(xkl(2,1)*xkl(1,3)-xkl(1,1)*xkl(2,3))/vj
820                  ckl(2,1)=(xkl(3,1)*xkl(2,3)-xkl(2,1)*xkl(3,3))/vj
821                  ckl(3,1)=(xkl(2,1)*xkl(3,2)-xkl(2,2)*xkl(3,1))/vj
822                  ckl(3,2)=(xkl(3,1)*xkl(1,2)-xkl(1,1)*xkl(3,2))/vj
823!
824!                 converting the Lagrangian strain into Eulerian
825!                 strain (only for deformation plasticity)
826!
827                  cauchy=0
828                  call str2mat(eloc,ckl,vj,cauchy)
829               endif
830!
831            endif
832!
833!                 calculating fields for incremental plasticity
834!
835            if(kode.le.-100) then
836!
837!              calculating the deformation gradient at the
838!              start of the increment
839!
840!              calculating the displacement gradient at the
841!              start of the increment
842!
843               do m2=1,3
844                  do m3=1,3
845                     vikl(m2,m3)=0.d0
846                  enddo
847               enddo
848!
849               do m1=1,nope
850                  do m2=1,3
851                     do m3=1,3
852                        vikl(m2,m3)=vikl(m2,m3)
853     &                       +shp(m3,m1)*vini(m2,konl(m1))
854                     enddo
855                  enddo
856               enddo
857!
858!              calculating the deformation gradient of the old
859!              fields
860!
861               xikl(1,1)=vikl(1,1)+1.d0
862               xikl(2,2)=vikl(2,2)+1.d0
863               xikl(3,3)=vikl(3,3)+1.d0
864               xikl(1,2)=vikl(1,2)
865               xikl(1,3)=vikl(1,3)
866               xikl(2,3)=vikl(2,3)
867               xikl(2,1)=vikl(2,1)
868               xikl(3,1)=vikl(3,1)
869               xikl(3,2)=vikl(3,2)
870!
871!              calculating the Jacobian
872!
873               vij=xikl(1,1)*(xikl(2,2)*xikl(3,3)
874     &              -xikl(2,3)*xikl(3,2))
875     &              -xikl(1,2)*(xikl(2,1)*xikl(3,3)
876     &              -xikl(2,3)*xikl(3,1))
877     &              +xikl(1,3)*(xikl(2,1)*xikl(3,2)
878     &              -xikl(2,2)*xikl(3,1))
879!
880!              stresses at the start of the increment
881!
882               do m1=1,6
883                  stre(m1)=sti(m1,jj,i)
884               enddo
885!
886            endif
887!
888!                 prestress values
889!
890            if(iprestr.eq.1) then
891               do kk=1,6
892                  beta(kk)=-prestr(kk,jj,i)
893               enddo
894            else
895               do kk=1,6
896                  beta(kk)=0.d0
897               enddo
898            endif
899!
900            if(ithermal(1).ge.1) then
901!
902!              calculating the temperature difference in
903!              the integration point
904!
905               t0l=0.d0
906               t1l=0.d0
907               if(ithermal(1).eq.1) then
908                  if((lakonl(4:5).eq.'8 ').or.
909     &               (lakonl(4:5).eq.'8I')) then
910                     do i1=1,8
911                        t0l=t0l+t0(konl(i1))/8.d0
912                        t1l=t1l+t1(konl(i1))/8.d0
913                     enddo
914                  elseif(lakonl(4:6).eq.'20 ') then
915                     nopered=20
916                     call lintemp(t0,konl,nopered,jj,t0l)
917                     call lintemp(t1,konl,nopered,jj,t1l)
918                  elseif(lakonl(4:6).eq.'10T') then
919                     call linscal10(t0,konl,t0l,null,shp)
920                     call linscal10(t1,konl,t1l,null,shp)
921                  else
922                     do i1=1,nope
923                        t0l=t0l+shp(4,i1)*t0(konl(i1))
924                        t1l=t1l+shp(4,i1)*t1(konl(i1))
925                     enddo
926                  endif
927               elseif(ithermal(1).ge.2) then
928                  if((lakonl(4:5).eq.'8 ').or.
929     &               (lakonl(4:5).eq.'8I')) then
930                     do i1=1,8
931                        t0l=t0l+t0(konl(i1))/8.d0
932                        t1l=t1l+vold(0,konl(i1))/8.d0
933                     enddo
934                  elseif(lakonl(4:6).eq.'20 ') then
935                     nopered=20
936                     call lintemp_th0(t0,konl,nopered,jj,t0l,mi)
937                     call lintemp_th1(vold,konl,nopered,jj,t1l,mi)
938                  elseif(lakonl(4:6).eq.'10T') then
939                     call linscal10(t0,konl,t0l,null,shp)
940                     call linscal10(vold,konl,t1l,mi(2),shp)
941                  else
942                     do i1=1,nope
943                        t0l=t0l+shp(4,i1)*t0(konl(i1))
944                        t1l=t1l+shp(4,i1)*vold(0,konl(i1))
945                     enddo
946                  endif
947               endif
948               tt=t1l-t0l
949            endif
950!
951!                 calculating the coordinates of the integration point
952!                 for material orientation purposes (for cylindrical
953!                 coordinate systems)
954!
955            if((iorien.gt.0).or.(kode.le.-100)) then
956               do j=1,3
957                  pgauss(j)=0.d0
958                  do i1=1,nope
959                     pgauss(j)=pgauss(j)+shp(4,i1)*co(j,konl(i1))
960                  enddo
961               enddo
962            endif
963!
964!                 material data; for linear elastic materials
965!                 this includes the calculation of the stiffness
966!                 matrix
967!
968            istiff=0
969!
970            call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,
971     &           nalcon,imat,amat,iorien,pgauss,orab,ntmat_,
972     &           elas,rho,i,ithermal,alzero,mattyp,t0l,t1l,ihyper,
973     &           istiff,elconloc,eth,kode,plicon,nplicon,
974     &           plkcon,nplkcon,npmat_,plconloc,mi(1),dtime,jj,
975     &           xstiff,ncmat_)
976!
977!           determining the mechanical strain
978!
979            if(ithermal(1).ne.0) then
980               do m1=1,6
981                  emec(m1)=eloc(m1)-eth(m1)
982c                 emec0(m1)=emeini(m1,jj,i)
983               enddo
984            else
985               do m1=1,6
986                  emec(m1)=eloc(m1)
987c                 emec0(m1)=emeini(m1,jj,i)
988               enddo
989            endif
990            if(kode.le.-100) then
991               do m1=1,6
992                  emec0(m1)=emeini(m1,jj,i)
993               enddo
994            endif
995!
996!           subtracting the plastic initial strains
997!
998            if(iprestr.eq.2) then
999               do m1=1,6
1000                  emec(m1)=emec(m1)-prestr(m1,jj,i)
1001               enddo
1002            endif
1003!
1004!           calculating the local stiffness and stress
1005!
1006            nlgeom_undo=0
1007            call mechmodel(elconloc,elas,emec,kode,emec0,ithermal,
1008     &           icmd,beta,stre,xkl,ckl,vj,xikl,vij,
1009     &           plconloc,xstate,xstateini,ielas,
1010     &           amat,t1l,dtime,time,ttime,i,jj,nstate_,mi(1),
1011     &           iorien,pgauss,orab,eloc,mattyp,qa(3),istep,iinc,
1012     &           ipkon,nmethod,iperturb,qa(4),nlgeom_undo)
1013!
1014            if(((nmethod.ne.4).or.(iperturb(1).ne.0)).and.
1015     &         (nmethod.ne.5)) then
1016               if(idesvar.eq.0) then
1017                  do m1=1,21
1018                     xstiff(m1,jj,i)=elas(m1)
1019                  enddo
1020               elseif(icoordinate.ne.1) then
1021!
1022!                 for orientation design variables: store the modified
1023!                 stiffness matrix for use in mafillsmse
1024!
1025                  idir=idesvar-3*((idesvar-1)/3)
1026                  do m1=1,21
1027                     dxstiff(m1,jj,i,idir)=elas(m1)
1028                  enddo
1029               endif
1030            endif
1031!
1032            if((iperturb(1).eq.-1).and.(nlgeom_undo.eq.0)) then
1033!
1034!                    if the forced displacements were changed at
1035!                    the start of a nonlinear step, the nodal
1036!                    forces due do this displacements are
1037!                    calculated in a purely linear way, and
1038!                    the first iteration is purely linear in order
1039!                    to allow the displacements to redistribute
1040!                    in a quasi-static way (only applies to
1041!                    quasi-static analyses (*STATIC))
1042!
1043               eloc(1)=exx-vokl(1,1)
1044               eloc(2)=eyy-vokl(2,2)
1045               eloc(3)=ezz-vokl(3,3)
1046               eloc(4)=exy-(vokl(1,2)+vokl(2,1))
1047               eloc(5)=exz-(vokl(1,3)+vokl(3,1))
1048               eloc(6)=eyz-(vokl(2,3)+vokl(3,2))
1049!
1050               if(mattyp.eq.1) then
1051                  e=elas(1)
1052                  un=elas(2)
1053                  um=e/(1.d0+un)
1054                  al=un*um/(1.d0-2.d0*un)
1055                  um=um/2.d0
1056                  am1=al*(eloc(1)+eloc(2)+eloc(3))
1057                  stre(1)=am1+2.d0*um*eloc(1)
1058                  stre(2)=am1+2.d0*um*eloc(2)
1059                  stre(3)=am1+2.d0*um*eloc(3)
1060                  stre(4)=um*eloc(4)
1061                  stre(5)=um*eloc(5)
1062                  stre(6)=um*eloc(6)
1063               elseif(mattyp.eq.2) then
1064                  stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2)
1065     &                 +eloc(3)*elas(4)
1066                  stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3)
1067     &                 +eloc(3)*elas(5)
1068                  stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5)
1069     &                 +eloc(3)*elas(6)
1070                  stre(4)=eloc(4)*elas(7)
1071                  stre(5)=eloc(5)*elas(8)
1072                  stre(6)=eloc(6)*elas(9)
1073               elseif(mattyp.eq.3) then
1074                  stre(1)=eloc(1)*elas(1)+eloc(2)*elas(2)+
1075     &                 eloc(3)*elas(4)+eloc(4)*elas(7)+
1076     &                 eloc(5)*elas(11)+eloc(6)*elas(16)
1077                  stre(2)=eloc(1)*elas(2)+eloc(2)*elas(3)+
1078     &                 eloc(3)*elas(5)+eloc(4)*elas(8)+
1079     &                 eloc(5)*elas(12)+eloc(6)*elas(17)
1080                  stre(3)=eloc(1)*elas(4)+eloc(2)*elas(5)+
1081     &                 eloc(3)*elas(6)+eloc(4)*elas(9)+
1082     &                 eloc(5)*elas(13)+eloc(6)*elas(18)
1083                  stre(4)=eloc(1)*elas(7)+eloc(2)*elas(8)+
1084     &                 eloc(3)*elas(9)+eloc(4)*elas(10)+
1085     &                 eloc(5)*elas(14)+eloc(6)*elas(19)
1086                  stre(5)=eloc(1)*elas(11)+eloc(2)*elas(12)+
1087     &                 eloc(3)*elas(13)+eloc(4)*elas(14)+
1088     &                 eloc(5)*elas(15)+eloc(6)*elas(20)
1089                  stre(6)=eloc(1)*elas(16)+eloc(2)*elas(17)+
1090     &                 eloc(3)*elas(18)+eloc(4)*elas(19)+
1091     &                 eloc(5)*elas(20)+eloc(6)*elas(21)
1092               endif
1093            endif
1094!
1095            skl(1,1)=stre(1)
1096            skl(2,2)=stre(2)
1097            skl(3,3)=stre(3)
1098            skl(2,1)=stre(4)
1099            skl(3,1)=stre(5)
1100            skl(3,2)=stre(6)
1101!
1102            skl(1,2)=skl(2,1)
1103            skl(1,3)=skl(3,1)
1104            skl(2,3)=skl(3,2)
1105!
1106!                 calculation of the nodal forces
1107!
1108            if(calcul_fn.eq.1)then
1109               if(idesvar.eq.0) then
1110!
1111!                 unperturbed state
1112!
1113                  do m1=1,nope
1114                     do m2=1,3
1115!
1116!                       linear elastic part
1117!
1118                        do m3=1,3
1119                           fn0(m2,indexe+m1)=fn0(m2,indexe+m1)+
1120     &                          xsj*skl(m2,m3)*shp(m3,m1)*weight
1121                        enddo
1122!
1123!                       nonlinear geometric part
1124!
1125                        if((iperturb(2).eq.1).and.(nlgeom_undo.eq.0))
1126     &                     then
1127                           do m3=1,3
1128                              do m4=1,3
1129                                 fn0(m2,indexe+m1)=fn0(m2,indexe+m1)+
1130     &                                xsj*skl(m4,m3)*weight*
1131     &                                (vkl(m2,m4)*shp(m3,m1)+
1132     &                                vkl(m2,m3)*shp(m4,m1))/2.d0
1133                              enddo
1134                           enddo
1135                        endif
1136!
1137                     enddo
1138                  enddo
1139c     Bernhardi start
1140                  if(lakonl(1:5).eq.'C3D8R') then
1141                     ahr=elas(1)*a
1142!
1143                     do i1=1,3
1144                        do k=1,4
1145                           hglf(i1,k)=0.0d0
1146                           do j=1,8
1147                              hglf(i1,k)=hglf(i1,k)+gs(j,k)*vl(i1,j)
1148                           enddo
1149                           hglf(i1,k)=hglf(i1,k)*ahr
1150                        enddo
1151                     enddo
1152                     do i1=1,3
1153                        do j=1,8
1154                           do k=1,4
1155                              fn0(i1,indexe+j)=fn0(i1,indexe+j)+
1156     &                                      hglf(i1,k)*gs(j,k)
1157                           enddo
1158                        enddo
1159                     enddo
1160                  endif
1161c     Bernhardi end
1162               else
1163!
1164!                 perturbed state
1165!
1166                  do m1=1,nope
1167                     do m2=1,3
1168!
1169!                       linear elastic part
1170!
1171                        do m3=1,3
1172                           dfn(m2,konl(m1))=dfn(m2,konl(m1))+
1173     &                          xsj*skl(m2,m3)*shp(m3,m1)*weight
1174                        enddo
1175!
1176!                       nonlinear geometric part
1177!
1178                        if((iperturb(2).eq.1).and.(nlgeom_undo.eq.0))
1179     &                     then
1180                           do m3=1,3
1181                              do m4=1,3
1182                                 dfn(m2,konl(m1))=dfn(m2,konl(m1))+
1183     &                                xsj*skl(m4,m3)*weight*
1184     &                                (vkl(m2,m4)*shp(m3,m1)+
1185     &                                vkl(m2,m3)*shp(m4,m1))/2.d0
1186                              enddo
1187                           enddo
1188                        endif
1189!
1190                     enddo
1191                  enddo
1192c     Bernhardi start
1193                  if(lakonl(1:5).eq.'C3D8R') then
1194!
1195                     ahr=elas(1)*a
1196!
1197                     do i1=1,3
1198                        do k=1,4
1199                           hglf(i1,k)=0.0d0
1200                           do j=1,8
1201                              hglf(i1,k)=hglf(i1,k)+gs(j,k)*vl(i1,j)
1202                           enddo
1203                           hglf(i1,k)=hglf(i1,k)*ahr
1204                        enddo
1205                     enddo
1206                     do i1=1,3
1207                        do j=1,8
1208                           do k=1,4
1209                              dfn(i1,konl(j))=dfn(i1,konl(j))
1210     &                             +hglf(i1,k)*gs(j,k)
1211                           enddo
1212                        enddo
1213                     enddo
1214                  endif
1215c     Bernhardi end
1216               endif
1217            endif
1218!
1219         enddo   ! enddo loop over the integration points
1220!
1221!     subtracting the unperturbed state of the element at stake
1222!
1223         if(idesvar.gt.0) then
1224            do m1=1,nope
1225               do m2=1,3
1226                  dfn(m2,konl(m1))=dfn(m2,konl(m1))
1227     &                 -fn0(m2,indexe+m1)
1228               enddo
1229            enddo
1230         endif
1231!
1232!     end of loop over all elements in thread
1233      enddo
1234!
1235      return
1236      end
1237