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!     construction of the b matrix
20!
21      subroutine resultnet(itg,ieg,ntg,
22     &     bc,nload,sideload,nelemload,xloadact,lakon,ntmat_,
23     &     v,shcon,nshcon,ipkon,kon,co,nflow,
24     &     iinc,istep,dtime,ttime,time,
25     &     ikforc,ilforc,xforcact,nforc,ielmat,nteq,prop,ielprop,
26     &     nactdog,nacteq,iin,physcon,camt,camf,camp,rhcon,nrhcon,
27     &     ipobody,ibody,xbodyact,nbody,dtheta,vold,xloadold,
28     &     reltime,nmethod,set,mi,ineighe,cama,vamt,vamf,vamp,vama,
29     &     nmpc,nodempc,ipompc,coefmpc,labmpc,iaxial,qat,qaf,ramt,
30     &     ramf,ramp,cocon,ncocon,iponoel,inoel,iplausi)
31!
32      implicit none
33!
34      logical identity
35      character*8 lakonl,lakon(*)
36      character*20 sideload(*),labmpc(*)
37      character*81 set(*)
38!
39      integer mi(*),itg(*),ieg(*),ntg,nteq,nflow,nload,
40     &     ielmat(mi(3),*),iflag,ider,iaxial,nalt,nalf,
41     &     nelemload(2,*),nope,nopes,mint2d,i,j,k,m,nrhcon(*),
42     &     node,imat,ntmat_,id,ifaceq(8,6),ifacet(6,4),numf,
43     &     ifacew(8,5),node1,node2,nshcon(*),nelem,ig,index,konl(20),
44     &     ipkon(*),kon(*),idof,ineighe(*),idir,ncocon(2,*),
45     &     iinc,istep,jltyp,nfield,ikforc(*),ipobody(2,*),
46     &     ilforc(*),nforc,nodem,idirf(8),ieq,nactdog(0:3,*),nbody,
47     &     nacteq(0:3,*),ielprop(*),nodef(8),iin,kflag,ibody(3,*),icase,
48     &     inv, index2,nmethod,nelem0,nodem0,nelem1,nodem1,nelem2,
49     &     nodem2,nelemswirl,nmpc,nodempc(3,*),ipompc(*),
50     &     iponoel(*),inoel(2,*),iplausi,indexe
51!
52      real*8 bc(nteq),xloadact(2,*),cp,h(2),physcon(*),r,dvi,rho,
53     &     xl2(3,8),coords(3),dxsj2,temp,xi,et,weight,xsj2(3),
54     &     gastemp,v(0:mi(2),*),shcon(0:3,ntmat_,*),co(3,*),shp2(7,8),
55     &     field,prop(*),tg1,tg2,dtime,ttime,time,g(3),eta,
56     &     xforcact(*),areaj,xflow,tvar(2),f,df(8),camt(*),camf(*),
57     &     camp(*),tl2(8),cama(*),vamt,vamf,vamp,vama,term,
58     &     rhcon(0:1,ntmat_,*),xbodyact(7,*),sinktemp,kappa,A,T,Tt,pt,
59     &     dtheta,ts1,ts2,xs2(3,7),pt1,pt2,Qred_crit,xflow_crit,
60     &     xflow0,xflow1,reltime,coefmpc(*),qat,qaf,ramt,ramf,ramp,
61     &     xflow2,heat,cocon(0:6,ntmat_,*),
62     &     Cp_cor,U,Ct,vold(0:mi(2),*),xloadold(2,*),bnac,
63     &     xflow360,Tt1,Tt2,T1,T2,heatnod,heatfac,A2,d,l,s
64!
65      include "gauss.f"
66!
67      data ifaceq /4,3,2,1,11,10,9,12,
68     &     5,6,7,8,13,14,15,16,
69     &     1,2,6,5,9,18,13,17,
70     &     2,3,7,6,10,19,14,18,
71     &     3,4,8,7,11,20,15,19,
72     &     4,1,5,8,12,17,16,20/
73      data ifacet /1,3,2,7,6,5,
74     &     1,2,4,5,9,8,
75     &     2,3,4,6,10,9,
76     &     1,4,3,8,10,7/
77      data ifacew /1,3,2,9,8,7,0,0,
78     &     4,5,6,10,11,12,0,0,
79     &     1,2,5,4,7,14,10,13,
80     &     2,3,6,5,8,15,11,14,
81     &     4,6,3,1,12,15,9,13/
82      data iflag /2/
83!
84      kflag=2
85      ider=0
86!
87      tvar(1)=time
88      tvar(2)=ttime+time
89!
90      heatnod=0.d0
91!
92!     calculating the maximum correction to the solution in the
93!     present network iteration
94!
95      camt(1)=0.d0
96      camf(1)=0.d0
97      camp(1)=0.d0
98      cama(1)=0.d0
99!
100      camt(2)=0.5d0
101      camf(2)=0.5d0
102      camp(2)=0.5d0
103      cama(2)=0.5d0
104!
105!     maximum change in the solution since the start of the
106!     network iterations (= entering radflowload)
107!
108      vamt=0.d0
109      vamf=0.d0
110      vamp=0.d0
111      vama=0.d0
112!
113      do i=1,ntg
114        node=itg(i)
115        do j=0,3
116          if(nactdog(j,node).eq.0) cycle
117          idof=nactdog(j,node)
118!
119          if(dabs(v(j,node)).gt.1.d-30) then
120            if(dabs(bc(idof)/v(j,node)).lt.1.d-13) bc(idof)=0.d0
121          endif
122          bnac=bc(idof)
123!
124          if(j.eq.0) then
125            if(dabs(bnac).gt.dabs(camt(1))) then
126              camt(1)=bnac
127              camt(2)=node+0.5d0
128            endif
129          elseif(j.eq.1) then
130            if(dabs(bnac).gt.dabs(camf(1))) then
131              camf(1)=bnac
132              camf(2)=node+0.5d0
133            endif
134          elseif(j.eq.2) then
135            if(dabs(bnac).gt.dabs(camp(1))) then
136              camp(1)=bnac
137              camp(2)=node+0.5d0
138            endif
139          else
140            if(dabs(bnac).gt.dabs(cama(1))) then
141              cama(1)=bnac
142              cama(2)=node+0.5d0
143            endif
144          endif
145        enddo
146      enddo
147!
148!     updating v
149!     vold is the solution at the entry of radflowload (= at
150!     the start of the network iterations)
151!
152      do i=1,ntg
153        node=itg(i)
154        do j=0,2
155          if(nactdog(j,node).eq.0) cycle
156          v(j,node)=v(j,node)+bc(nactdog(j,node))*dtheta
157!
158          if((j.eq.0).and.(dabs(v(j,node)-vold(j,node)).gt.vamt)) then
159            vamt=dabs(v(j,node)-vold(j,node))
160          elseif((j.eq.1).and.
161     &           (dabs(v(j,node)-vold(j,node)).gt.vamf)) then
162            vamf=dabs(v(j,node)-vold(j,node))
163          elseif((j.eq.2).and.
164     &           (dabs(v(j,node)-vold(j,node)).gt.vamp)) then
165            vamp=dabs(v(j,node)-vold(j,node))
166          endif
167!
168        enddo
169      enddo
170c     write(*,*) 'resultnet'
171c     write(*,*) 'mass flow in node 3: ',v(1,3)
172c     write(*,*) 'pressure in node 4: ',v(2,4)
173c     write(*,*) 'mass flow in node 5: ',v(1,5)
174!
175!     update geometry changes
176!
177      do i=1,nflow
178        if(lakon(ieg(i))(6:7).eq.'GV') then
179          index=ipkon(ieg(i))
180          node=kon(index+2)
181          if(nactdog(3,node).eq.0) cycle
182          index=ielprop(ieg(i))
183          v(3,node)=v(3,node)+bc(nactdog(3,node))*dtheta
184          if(v(3,node).gt.0.99999d0) then
185            v(3,node)=0.99999d0
186          elseif(v(3,node).lt.0.12501d0) then
187            v(3,node)=0.12501d0
188          endif
189          if(dabs(v(3,node)).gt.vama) vama=dabs(v(3,node))
190!
191!     Geometry factor of an ACC tube
192!     for all versions of the ACC tube
193!
194        elseif(lakon(ieg(i))(2:7).eq.'ACCTUB') then
195          index=ipkon(ieg(i))
196          node=kon(index+2)
197          if(nactdog(3,node).eq.0) cycle
198          index=ielprop(ieg(i))
199!     Interval factor unknown
200          if(nint(prop(index+1)).eq.2) then
201!     Using a smaller step gets better convergence(factor 0.5)
202            v(3,node)=v(3,node)+bc(nactdog(3,node))*
203     &           dtheta
204            if(v(3,node).lt.0.1d0)then
205              v(3,node) = 0.1d0
206            endif
207!     Hole diameter factor unknown
208          elseif(nint(prop(index+1)).eq.3) then
209            v(3,node)=v(3,node)+bc(nactdog(3,node))*dtheta
210            if(v(3,node).lt.0.1d0)then
211              v(3,node) = 0.1d0
212            endif
213          endif
214          if(dabs(v(3,node)).gt.vama) vama=dabs(v(3,node))
215!
216!     update location of hydraulic jump
217!
218        elseif(lakon(ieg(i))(2:5).eq.'LICH') then
219          if((lakon(ieg(i))(6:7).eq.'SG').or.
220     &         (lakon(ieg(i))(6:7).eq.'WE').or.
221     &         (lakon(ieg(i))(6:7).eq.'DS')) then
222            index=ipkon(ieg(i))
223            node=kon(index+2)
224            if(nactdog(3,node).eq.0) cycle
225            index=ielprop(ieg(i))
226            if(lakon(ieg(i))(6:7).eq.'SG') then
227              eta=prop(index+4)+bc(nactdog(3,node))*dtheta
228              prop(index+4)=eta
229              nelem=nint(prop(index+7))
230            elseif(lakon(ieg(i))(6:7).eq.'WE') then
231              eta=prop(index+4)+bc(nactdog(3,node))*dtheta
232              prop(index+4)=eta
233              nelem=nint(prop(index+7))
234            elseif(lakon(ieg(i))(6:7).eq.'DS') then
235              eta=prop(index+7)+bc(nactdog(3,node))*dtheta
236              prop(index+7)=eta
237              nelem=nint(prop(index+9))
238            endif
239            v(3,node)=eta
240            vama=eta
241!
242!     check whether 0<=eta<=1. If not, the hydraulic jump
243!     does not take place in the element itself and has to
244!     be forced out of the element by adjusting the
245!     water depth of one of the end nodes
246!
247            if((eta.lt.0.d0).or.(eta.gt.1.d0)) then
248              index=ipkon(nelem)
249              node1=kon(index+1)
250              nodem=kon(index+2)
251              node2=kon(index+3)
252              xflow=v(1,nodem)
253!
254!     determining the temperature for the
255!     material properties
256!
257              if(xflow.gt.0) then
258                if(node1.eq.0) then
259                  gastemp=v(0,node2)
260                else
261                  gastemp=v(0,node1)
262                endif
263              else
264                if(node2.eq.0) then
265                  gastemp=v(0,node1)
266                else
267                  gastemp=v(0,node2)
268                endif
269              endif
270!
271              imat=ielmat(1,nelem)
272!
273              call materialdata_tg(imat,ntmat_,gastemp,shcon,nshcon,
274     &             cp,r,dvi,rhcon,nrhcon,rho)
275!
276              do j=1,3
277                g(j)=0.d0
278              enddo
279              if(nbody.gt.0) then
280                index=nelem
281                do
282                  j=ipobody(1,index)
283                  if(j.eq.0) exit
284                  if(ibody(1,j).eq.2) then
285                    g(1)=g(1)+xbodyact(1,j)*xbodyact(2,j)
286                    g(2)=g(2)+xbodyact(1,j)*xbodyact(3,j)
287                    g(3)=g(3)+xbodyact(1,j)*xbodyact(4,j)
288                  endif
289                  index=ipobody(2,index)
290                  if(index.eq.0) exit
291                enddo
292              endif
293!
294              kflag=3
295              call flux(node1,node2,nodem,nelem,lakon,kon,ipkon,
296     &             nactdog,identity,
297     &             ielprop,prop,kflag,v,xflow,f,nodef,idirf,df,
298     &             cp,r,rho,physcon,g,co,dvi,numf,vold,set,shcon,
299     &             nshcon,rhcon,nrhcon,ntmat_,mi,ider,ttime,time,
300     &             iaxial,iplausi)
301              kflag=2
302!
303            endif
304          endif
305        endif
306      enddo
307!
308c     write(*,*) 'resultnet'
309c     do i=1,ntg
310c     node=itg(i)
311c     write(*,'(i10,3(1x,e11.4))') node,(v(j,node),j=0,2)
312c     enddo
313!
314!     testing the validity of the pressures
315!
316      do i=1,ntg
317        node=itg(i)
318        if(v(2,node).lt.0) then
319          write(*,*) 'wrong pressure; node ',node,v(2,node)
320          iin=0
321          return
322        endif
323      enddo
324!
325!     testing validity of temperatures
326!
327      do i=1,ntg
328        node=itg(i)
329        if(v(0,node).lt.0) then
330          write(*,*) 'wrong temperature; node ',node,v(0,node)
331          iin=0
332          return
333        endif
334      enddo
335!
336!     testing the validity of the solution for branches elements
337!     and restrictor. Since the element properties are dependent on
338!     a predefined flow direction a change of this will lead to
339!     wrong head losses
340!
341      do i=1,nflow
342        nelem=ieg(i)
343        if ((lakon(nelem)(4:5).eq.'ATR').or.
344     &       (lakon(nelem)(4:5).eq.'RTA')) then
345          xflow=v(1,kon(ipkon(nelem)+2))
346          if(xflow.lt.0d0)then
347            Write(*,*)'*WARNING in resultnet.f'
348            write(*,*)'Element',nelem,'of TYPE ABSOLUTE TO RELATIVE'
349            write(*,*)'The flow direction is no more conform '
350            write(*,*)'to element definition'
351            write(*,*)'Check the pertinence of the results'
352          endif
353        elseif(lakon(nelem)(2:3).eq.'RE') then
354!
355          if(lakon(nelem)(4:5).ne.'BR') then
356            nodem=kon(ipkon(nelem)+2)
357            xflow=v(1,nodem)
358            if (xflow.lt.0) then
359              Write(*,*)'*WARNING in resultnet.f'
360              write(*,*)'Element',nelem,'of TYPE RESTRICTOR'
361              write(*,*)'The flow direction is no more conform '
362              write(*,*)'to element definition'
363              write(*,*)'Check the pertinence of the results'
364            endif
365!
366          elseif(lakon(nelem)(4:5).eq.'BR') then
367            index=ielprop(nelem)
368!
369            nelem0=nint(prop(index+1))
370            nodem0=kon(ipkon(nelem0)+2)
371            xflow0=v(1,nodem0)
372!
373            nelem1=nint(prop(index+2))
374            nodem1=kon(ipkon(nelem1)+2)
375            xflow1=v(1,nodem1)
376!
377            nelem2=nint(prop(index+3))
378            nodem2=kon(ipkon(nelem2)+2)
379            xflow2=v(1,nodem2)
380!
381            if((xflow0.lt.0).or.(xflow1.lt.0).or.(xflow2.lt.0)) then
382              Write(*,*)'*WARNING in resultnet.f'
383              write(*,*)'Element',nelem,'of TYPE BRANCH'
384              write(*,*)'The flow direction is no more conform '
385              write(*,*)'to element definition'
386              write(*,*)'Check the pertinence of the results'
387            endif
388          endif
389        endif
390      enddo
391!
392!     determining the static temperature and checking for
393!     critical conditions (only for non-chambers, i.e. for
394!     nodes purely belonging to gas pipes or restrictors except
395!     restrictor wall orifice)
396!
397      iplausi=1
398      do i=1,ntg
399        node=itg(i)
400        nelem=ineighe(i)
401        if(nelem.eq.-1) then
402          v(3,node)=v(0,node)
403        elseif(nelem.gt.0) then
404!
405          nodem=kon(ipkon(nelem)+2)
406          T=v(3,node)
407          Tt=v(0,node)
408          pt=v(2,node)
409          xflow=v(1,nodem)
410!
411          icase=0
412          inv=1
413          imat=ielmat(1,nelem)
414          call materialdata_tg(imat,ntmat_,v(3,node),
415     &         shcon,nshcon,cp,r,dvi,rhcon,nrhcon,rho)
416!
417          index=ielprop(nelem)
418          indexe=ipkon(nelem)
419!
420          kappa=(cp/(cp-R))
421!
422          if(lakon(nelem)(2:5).eq.'GAPF') then
423!
424!     distinguish between standard pipes and flexible
425!     pipes
426!
427            if((lakon(nelem)(6:8).eq.'AFR').or.
428     &           (lakon(nelem)(6:8).eq.'ARL').or.
429     &           (lakon(nelem)(6:8).eq.'ARG').or.
430     &           (lakon(nelem)(6:8).eq.'IFR').or.
431     &           (lakon(nelem)(6:8).eq.'IRL')) then
432              call calcgeomelemnet(vold,co,prop,lakon,nelem,
433     &             ttime,time,ielprop,mi,A,A2,d,l,s)
434            else
435              A=prop(index+1)
436            endif
437!
438            if(lakon(nelem)(2:6).eq.'GAPFA') then
439              icase=0
440            elseif(lakon(nelem)(2:6).eq.'GAPFI') then
441              icase=1
442            endif
443!
444          elseif(lakon(nelem)(2:5).eq.'GAPR') then
445!
446!     rotating adiabatic gas pipe with variable cross section:
447!     check whether section 1 or 2
448!
449            if(kon(ipkon(nelem)+1).eq.node) then
450              A=prop(index+1)
451            else
452              A=prop(index+2)
453            endif
454            icase=0
455!
456          elseif(lakon(nelem)(2:3).eq.'RE') then
457            node1=kon(indexe+1)
458            node2=kon(indexe+3)
459!
460            if(lakon(nelem)(4:5).eq.'EX') then
461              if(lakon(nint(prop(index+4)))(2:6).eq.'GAPFA') then
462                icase=0
463              elseif(lakon(nint(prop(index+4)))(2:6).eq.'GAPFI')then
464                icase=1
465              endif
466            else
467              icase=0
468            endif
469!
470!     defining the sections
471!
472            if(lakon(nelem)(4:5).eq.'BE') then
473              A=prop(index+1)
474!
475            elseif(lakon(nelem)(4:5).eq.'BR') then
476              if(lakon(nelem)(4:6).eq.'BRJ') then
477                if(nelem.eq.nint(prop(index+2)))then
478                  A=prop(index+5)
479                elseif(nelem.eq.nint(prop(index+3))) then
480                  A=prop(index+6)
481                endif
482              elseif(lakon(nelem)(4:6).eq.'BRS') then
483                if(nelem.eq.nint(prop(index+2)))then
484                  A=prop(index+5)
485                elseif(nelem.eq.nint(prop(index+3))) then
486                  A=prop(index+6)
487                endif
488              endif
489!
490            else
491!
492              if((lakon(nelem)(2:5).eq.'RBSF')) then
493                call calcgeomelemnet(vold,co,prop,lakon,nelem,
494     &               ttime,time,ielprop,mi,A,A2,d,l,s)
495                if(node.eq.node2) then
496                  A=A2
497                endif
498              else
499                if(node.eq.node1) then
500                  A=prop(index+1)
501                elseif(node.eq.node2) then
502                  A=prop(index+2)
503                endif
504              endif
505            endif
506          elseif(lakon(nelem)(2:3).eq.'UP') then
507!
508!     user elements whose names start with UP are assumed
509!     to be Pipe-like (static and total temperatures differ)
510!
511            call calcgeomelemnet(vold,co,prop,lakon,nelem,
512     &           ttime,time,ielprop,mi,A,A2,d,l,s)
513c     A=prop(index+1)
514            icase=0
515          endif
516!
517          xflow360=dabs(xflow*iaxial)
518          call ts_calc(xflow360,Tt,pt,kappa,r,A,T,icase)
519!
520          v(3,node)=T
521!
522          if(lakon(nelem)(2:3).eq.'RE') then
523!
524!     check whether the mass flow does not exceed
525!     critical conditions (only for restrictor elements)
526!
527            if(xflow.lt.0d0) then
528              inv=-1
529            else
530              inv=1
531            endif
532!
533            if(icase.eq.0) then
534              Qred_crit=dsqrt(kappa/R)*(1.d0+0.5d0*(kappa-1.d0))
535     &             **(-0.5d0*(kappa+1.d0)/(kappa-1.d0))
536            else
537              Qred_crit=dsqrt(1/R)*(1.d0+0.5d0*(kappa-1.d0)/kappa)
538     &             **(-0.5d0*(kappa+1.d0)/(kappa-1.d0))
539            endif
540            xflow_crit=Qred_crit*pt*A/dsqrt(Tt)
541!
542            if(xflow360.ge.xflow_crit) then
543!
544!     next 3 lines: change on 20200929 (Thomas Petzke)
545!
546              if(dabs((xflow_crit-xflow360)/xflow_crit).gt.1.d-5)
547     &             then
548                iplausi=0
549              endif
550!
551              v(1,nodem)=inv*xflow_crit/iaxial
552              if(icase.eq.1) then
553!
554                if(nactdog(0,node2).ne.0) then
555                  node1=kon(indexe+1)
556                  node2=kon(indexe+3)
557                  v(3,node2)=v(3,node1)
558                  v(0,node2)=v(3,node2)
559     &                 *(1+0.5d0*(kappa-1)/kappa)
560
561                endif
562              endif
563            endif
564          endif
565!
566        endif
567!
568      enddo
569!
570!     testing if the flow direction is conform to the pressure
571!     gradient
572!
573c      iplausi=1
574      do i=1,nflow
575        nelem=ieg(i)
576        indexe=ipkon(nelem)
577        node1=kon(indexe+1)
578        nodem=kon(indexe+2)
579        node2=kon(indexe+3)
580!
581        if((node1.ne.0).and.(node2.ne.0)) then
582!
583!     the mass flow rates must always be oriented in the direction
584!     of the decreasing pressure gradient except for ACCTUBE,
585!     CROSS SPLIT, MOEHRING, ROTATING CAVITY, RADIAL OUTFLOW and
586!     INLET, RIMSEAL, S-PUMP and VORTEX
587!
588          if((lakon(nelem)(2:8).ne.'ACCTUBE').and.
589     &         (lakon(nelem)(2:5).ne.'ATR').and.
590     &         (lakon(nelem)(2:5).ne.'RTA').and.
591     &         (lakon(nelem)(2:5).ne.'CROS').and.
592     &         (lakon(nelem)(2:3).ne.'LI').and.
593     &         (lakon(nelem)(2:3).ne.'LP').and.
594     &         (lakon(nelem)(2:4).ne.'MRG').and.
595     &         (lakon(nelem)(2:4).ne.'RCV').and.
596     &         (lakon(nelem)(2:3).ne.'RO').and.
597     &         (lakon(nelem)(2:5).ne.'RIMS').and.
598     &         (lakon(nelem)(2:5).ne.'FDPF').and.
599     &         (lakon(nelem)(2:5).ne.'FCVF').and.
600     &         (lakon(nelem)(2:5).ne.'MFPC').and.
601     &         (lakon(nelem)(2:6).ne.'SPUMP').and.
602     &         (lakon(nelem)(2:5).ne.'GAPR').and.
603     &         (lakon(nelem)(2:3).ne.'VO')) then
604            if(nactdog(1,nodem).ne.0) then
605              if(((v(1,nodem).gt.0.d0).and.
606     &             (v(2,node1)/v(2,node2).lt.(1.d0-1.d-10))).or.
607     &             ((v(1,nodem).lt.0.d0).and.
608     &             (v(2,node2)/v(2,node1).lt.(1.d0-1.d-10)))) then
609c              if(((v(1,nodem).gt.0.d0).and.
610c     &             (v(2,node1).lt.v(2,node2))).or.
611c     &             ((v(1,nodem).lt.0.d0).and.
612c     &             (v(2,node1).gt.v(2,node2)))) then
613                iplausi=0
614                write(*,*) '*WARNING in resultnet: the flow'
615                write(*,*) '         direction is not conform'
616                write(*,*) '         to the pressure gradient'
617                write(*,*) '         element',nelem
618              endif
619            endif
620          endif
621        endif
622      enddo
623!
624!     reinitialisation of the Bc matrix
625!
626      do i=1,nteq
627        bc(i)=0.d0
628      enddo
629!
630      qat=0.d0
631      qaf=0.d0
632      nalt=0
633      nalf=0
634!
635!     determining the residual
636!
637      do i=1,nflow
638        nelem=ieg(i)
639        index=ipkon(nelem)
640        node1=kon(index+1)
641        nodem=kon(index+2)
642        node2=kon(index+3)
643        xflow=v(1,nodem)
644!
645!     gas: the property temperature is the static temperature
646!
647        if((lakon(nelem)(2:3).ne.'LP').and.
648     &       (lakon(nelem)(2:3).ne.'LI')) then
649          if(node1.eq.0) then
650            tg1=v(0,node2)
651            tg2=tg1
652            ts1=v(3,node2)
653            ts2=ts1
654          elseif(node2.eq.0) then
655            tg1=v(0,node1)
656            tg2=tg1
657            ts1=v(3,node1)
658            ts2=ts1
659          else
660            tg1=v(0,node1)
661            tg2=v(0,node2)
662            ts1=v(3,node1)
663            ts2=v(3,node2)
664          endif
665          gastemp=(ts1+ts2)/2.d0
666        else
667!
668!     liquid: only one temperature
669!
670          if(xflow.gt.0) then
671            if(node1.eq.0) then
672              gastemp=v(0,node2)
673            else
674              gastemp=v(0,node1)
675            endif
676          else
677            if(node2.eq.0) then
678              gastemp=v(0,node1)
679            else
680              gastemp=v(0,node2)
681            endif
682          endif
683!
684          if(node1.eq.0) then
685            tg2=v(0,node2)
686            tg1=tg2
687          elseif(node2.eq.0) then
688            tg1=v(0,node1)
689            tg2=tg1
690          else
691            tg1=v(0,node1)
692            tg2=v(0,node2)
693          endif
694        endif
695!
696        imat=ielmat(1,nelem)
697!
698        call materialdata_tg(imat,ntmat_,gastemp,shcon,nshcon,cp,r,dvi,
699     &       rhcon,nrhcon,rho)
700        kappa=Cp/(Cp-R)
701!
702!     Definitions of the constant for isothermal flow elements
703!
704        if(lakon(nelem)(2:6).eq.'GAPFI') then
705          if((node1.ne.0).and.(node2.ne.0)) then
706!
707            if((lakon(nelem)(7:8).eq.'FR').or.
708     &           (lakon(nelem)(7:8).eq.'RL')) then
709              call calcgeomelemnet(vold,co,prop,lakon,nelem,
710     &             ttime,time,ielprop,mi,A,A2,d,l,s)
711            else
712              A=prop(ielprop(nelem)+1)
713            endif
714!
715            pt1=v(2,node1)
716            pt2=v(2,node2)
717!
718            xflow360=xflow*iaxial
719            icase=1
720!
721            if(pt1.ge.pt2)then
722              if(dabs(tg2/ts2-(1.d0+0.5d0*(kappa-1)/kappa)).lt.1d-5)
723     &             then
724                pt2=dabs(xflow360)*dsqrt(Tg2*R)/A
725     &               *(1.d0+0.5d0*(kappa-1.d0)/kappa)
726     &               **(0.5d0*(kappa+1.d0)/(kappa-1.d0))
727!
728              endif
729              Tt1=v(0,node1)-physcon(1)
730              Tt2=v(0,node2)-physcon(1)
731              T1=v(3,node1)
732              call ts_calc(xflow360,Tt1,pt1,kappa,r,A,T1,icase)
733              call ts_calc(xflow360,Tt2,pt2,kappa,r,A,T2,icase)
734              v(3,node1)=T1
735              v(3,node2)=T2
736            else
737              pt1=v(2,node2)
738              pt2=v(2,node1)
739              if(v(2,nodem).ge.(pt2/pt1))then
740                pt2=v(2,nodem)*pt1
741              endif
742!
743              Tt1=v(0,node2)-physcon(1)
744              call ts_calc(xflow360,Tt1,pt1,kappa,r,A,T1,icase)
745              Tt2=v(0,node1)-physcon(1)
746              call ts_calc(xflow360,Tt2,pt2,kappa,r,A,T2,icase)
747            endif
748!
749          endif
750        endif
751!
752        if(node1.ne.0) then
753!
754!     energy equation contribution node1
755!
756          if (nacteq(0,node1).ne.0) then
757            ieq=nacteq(0,node1)
758!
759            if(nacteq(3,node1).eq.0) then
760              if (xflow.lt.0d0)then
761                term=cp*(tg1-tg2)*xflow
762                bc(ieq)=bc(ieq)+term
763                qat=qat+dabs(term)
764                nalt=nalt+1
765              endif
766!
767            elseif(lakon(nelem)(2:6).eq.'GAPFI') then
768              if((nacteq(3,node1).eq.node2)) then
769!
770                bc(ieq)=(T2-T1)
771!
772              endif
773            endif
774          endif
775!
776!     mass equation contribution node1
777!
778          if (nacteq(1,node1).ne.0) then
779            ieq=nacteq(1,node1)
780            bc(ieq)=bc(ieq)-xflow
781            qaf=qaf+dabs(xflow)
782            nalf=nalf+1
783          endif
784        endif
785!
786        if(node2.ne.0) then
787!
788!     energy equation contribution node2
789!
790          if (nacteq(0,node2).ne.0) then
791            ieq=nacteq(0,node2)
792!
793            if(nacteq(3,node2).eq.0) then
794              if (xflow.gt.0d0)then
795                term=cp*(tg2-tg1)*xflow
796                bc(ieq)=bc(ieq)-term
797                qat=qat+dabs(term)
798                nalt=nalt+1
799              endif
800!
801            elseif(lakon(nelem)(2:6).eq.'GAPFI') then
802              if(nacteq(3,node2).eq.node1) then
803!
804                bc(ieq)=(T2-T1)
805!
806              endif
807            endif
808          endif
809!
810!     mass equation contribution node2
811!     only in the case of an ACCTUBE but not in the case of an ACCTUBO
812!
813          if(lakon(nelem)(2:8).ne.'ACCTUBE') then
814            if (nacteq(1,node2).ne.0) then
815              ieq=nacteq(1,node2)
816              bc(ieq)=bc(ieq)+xflow
817              qaf=qaf+dabs(xflow)
818              nalf=nalf+1
819            endif
820          else
821            if (nacteq(1,node2).ne.0) then
822              if(nelem.ne.nint(prop(ielprop(nelem)+14))) then
823                ieq=nacteq(1,node2)
824                bc(ieq)=bc(ieq)+xflow-v(0,nodem)
825                qaf=qaf+dabs(xflow)
826                nalf=nalf+1
827              else
828                ieq=nacteq(1,node2)
829                bc(ieq)=bc(ieq)+xflow
830                qaf=qaf+dabs(xflow)
831                nalf=nalf+1
832              endif
833            endif
834          endif
835        endif
836!
837!     element equation
838!
839        if (nacteq(2,nodem).ne.0) then
840          ieq=nacteq (2,nodem)
841!
842!     for liquids: determine the gravity vector
843!
844          if((lakon(nelem)(2:3).eq.'LI').or.
845     &         (lakon(nelem)(2:3).eq.'LP')) then
846            do j=1,3
847              g(j)=0.d0
848            enddo
849            if(nbody.gt.0) then
850              index=nelem
851              do
852                j=ipobody(1,index)
853                if(j.eq.0) exit
854                if(ibody(1,j).eq.2) then
855                  g(1)=g(1)+xbodyact(1,j)*xbodyact(2,j)
856                  g(2)=g(2)+xbodyact(1,j)*xbodyact(3,j)
857                  g(3)=g(3)+xbodyact(1,j)*xbodyact(4,j)
858                endif
859                index=ipobody(2,index)
860                if(index.eq.0) exit
861              enddo
862            endif
863          endif
864!
865          call flux(node1,node2,nodem,nelem,lakon,kon,ipkon,
866     &         nactdog,identity,
867     &         ielprop,prop,kflag,v,xflow,f,nodef,idirf,df,
868     &         cp,r,rho,physcon,g,co,dvi,numf,vold,set,shcon,
869     &         nshcon,rhcon,nrhcon,ntmat_,mi,ider,ttime,time,
870     &         iaxial,iplausi)
871          bc(ieq)=-f
872        endif
873      enddo
874!
875!     convection with the walls: contribution to the energy equations
876!
877      do i=1,nload
878        if(sideload(i)(3:4).eq.'FC') then
879          nelem=nelemload(1,i)
880          index=ipkon(nelem)
881          if(index.lt.0) cycle
882          lakonl=lakon(nelem)
883          node=nelemload(2,i)
884          ieq=nacteq(0,node)
885          if(ieq.eq.0) then
886            cycle
887          endif
888!
889          call nident(itg,node,ntg,id)
890!
891!     calculate the area
892!
893          read(sideload(i)(2:2),'(i1)') ig
894!
895!     number of nodes and integration points in the face
896!
897          if(lakonl(4:4).eq.'2') then
898            nope=20
899            nopes=8
900          elseif(lakonl(4:4).eq.'8') then
901            nope=8
902            nopes=4
903          elseif(lakonl(4:5).eq.'10') then
904            nope=10
905            nopes=6
906          elseif(lakonl(4:4).eq.'4') then
907            nope=4
908            nopes=3
909          elseif(lakonl(4:5).eq.'15') then
910            nope=15
911          else
912            nope=6
913          endif
914!
915          if(lakonl(4:5).eq.'8R') then
916            mint2d=1
917          elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
918     &           then
919            if(lakonl(7:7).eq.'A') then
920              mint2d=2
921            else
922              mint2d=4
923            endif
924          elseif(lakonl(4:4).eq.'2') then
925            mint2d=9
926          elseif(lakonl(4:5).eq.'10') then
927            mint2d=3
928          elseif(lakonl(4:4).eq.'4') then
929            mint2d=1
930          endif
931!
932          if(lakonl(4:4).eq.'6') then
933            mint2d=1
934            if(ig.le.2) then
935              nopes=3
936            else
937              nopes=4
938            endif
939          endif
940          if(lakonl(4:5).eq.'15') then
941            if(ig.le.2) then
942              mint2d=3
943              nopes=6
944            else
945              mint2d=4
946              nopes=8
947            endif
948          endif
949!
950!     connectivity of the element
951!
952          do k=1,nope
953            konl(k)=kon(index+k)
954          enddo
955!
956!     coordinates of the nodes belonging to the face
957!
958          if((nope.eq.20).or.(nope.eq.8)) then
959            do k=1,nopes
960              tl2(k)=v(0,konl(ifaceq(k,ig)))
961              do j=1,3
962                xl2(j,k)=co(j,konl(ifaceq(k,ig)))+
963     &               v(j,konl(ifaceq(k,ig)))
964              enddo
965            enddo
966          elseif((nope.eq.10).or.(nope.eq.4)) then
967            do k=1,nopes
968              tl2(k)=v(0,konl(ifacet(k,ig)))
969              do j=1,3
970                xl2(j,k)=co(j,konl(ifacet(k,ig)))+
971     &               v(j,konl(ifacet(k,ig)))
972              enddo
973            enddo
974          else
975            do k=1,nopes
976              tl2(k)=v(0,konl(ifacew(k,ig)))
977              do j=1,3
978                xl2(j,k)=co(j,konl(ifacew(k,ig)))+
979     &               v(j,konl(ifacew(k,ig)))
980              enddo
981            enddo
982          endif
983!
984!     integration to obtain the area and the mean
985!     temperature
986!
987          do m=1,mint2d
988            if((lakonl(4:5).eq.'8R').or.
989     &           ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
990              xi=gauss2d1(1,m)
991              et=gauss2d1(2,m)
992              weight=weight2d1(m)
993            elseif((lakonl(4:4).eq.'8').or.
994     &             (lakonl(4:6).eq.'20R').or.
995     &             ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
996              xi=gauss2d2(1,m)
997              et=gauss2d2(2,m)
998              weight=weight2d2(m)
999            elseif(lakonl(4:4).eq.'2') then
1000              xi=gauss2d3(1,m)
1001              et=gauss2d3(2,m)
1002              weight=weight2d3(m)
1003            elseif((lakonl(4:5).eq.'10').or.
1004     &             ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
1005              xi=gauss2d5(1,m)
1006              et=gauss2d5(2,m)
1007              weight=weight2d5(m)
1008            elseif((lakonl(4:4).eq.'4').or.
1009     &             ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
1010              xi=gauss2d4(1,m)
1011              et=gauss2d4(2,m)
1012              weight=weight2d4(m)
1013            endif
1014!
1015            if(nopes.eq.8) then
1016              call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
1017            elseif(nopes.eq.4) then
1018              call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
1019            elseif(nopes.eq.6) then
1020              call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
1021            else
1022              call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
1023            endif
1024!
1025            dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
1026     &           xsj2(3)*xsj2(3))
1027            areaj=dxsj2*weight
1028!
1029            temp=0.d0
1030            do k=1,3
1031              coords(k)=0.d0
1032            enddo
1033            do j=1,nopes
1034              temp=temp+tl2(j)*shp2(4,j)
1035              do k=1,3
1036                coords(k)=coords(k)+xl2(k,j)*shp2(4,j)
1037              enddo
1038            enddo
1039!
1040            sinktemp=v(0,node)
1041            if(sideload(i)(5:6).ne.'NU') then
1042              h(1)=xloadact(1,i)
1043            else
1044              read(sideload(i)(2:2),'(i1)') jltyp
1045              jltyp=jltyp+10
1046              call film(h,sinktemp,temp,istep,
1047     &             iinc,tvar,nelem,m,coords,jltyp,field,nfield,
1048     &             sideload(i),node,areaj,v,mi,ipkon,kon,lakon,
1049     &             iponoel,inoel,ielprop,prop,ielmat,shcon,nshcon,
1050     &             rhcon,nrhcon,ntmat_,cocon,ncocon,
1051     &             ipobody,xbodyact,ibody,heatnod,heatfac)
1052            endif
1053            if(lakonl(5:7).eq.'0RA') then
1054              term=2.d0*((temp-sinktemp)*h(1)+heatnod)*dxsj2*weight
1055              bc(ieq)=bc(ieq)+term
1056              qat=qat+dabs(term)
1057              nalt=nalt+1
1058            else
1059              term=((temp-sinktemp)*h(1)+heatnod)*dxsj2*weight
1060              bc(ieq)=bc(ieq)+term
1061              qat=qat+dabs(term)
1062              nalt=nalt+1
1063            endif
1064          enddo
1065        endif
1066      enddo
1067!
1068!     prescribed heat generation: contribution to the energy equations
1069!
1070      do i=1,ntg
1071        node=itg(i)
1072        idof=8*(node-1)
1073        call nident(ikforc,idof,nforc,id)
1074        if(id.gt.0) then
1075          if(ikforc(id).eq.idof) then
1076            ieq=nacteq(0,node)
1077            if(ieq.ne.0) then
1078              term=xforcact(ilforc(id))
1079              bc(ieq)=bc(ieq)+term
1080              qat=qat+dabs(term)
1081              nalt=nalt+1
1082            endif
1083            cycle
1084          endif
1085        endif
1086      enddo
1087!
1088!     in the case of forced vortices, when temperature change
1089!     is required, additional heat input is added in the energy
1090!     equation for the downstream node
1091!
1092      do i=1,nflow
1093        nelem=ieg(i)
1094!
1095!     free vortex and no temperature change
1096!
1097        if(lakon(nelem)(2:3).eq.'VO') then
1098          if(lakon(nelem)(4:5).eq.'FR') then
1099            if((nint(prop(ielprop(nelem)+8)).eq.0).or.
1100     &           (nint(prop(ielprop(nelem)+8)).eq.1)) cycle
1101          elseif(lakon(nelem)(4:5).eq.'FO') then
1102            if(nint(prop(ielprop(nelem)+6)).eq.0) cycle
1103          endif
1104!
1105          call calcheatnet(nelem,lakon,ipkon,kon,v,ielprop,prop,
1106     &         ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,ipobody,ibody,
1107     &         xbodyact,mi,nacteq,bc,qat,nalt)
1108!
1109        elseif(lakon(nelem)(2:5).eq.'GAPR') then
1110!
1111          call calcheatnet(nelem,lakon,ipkon,kon,v,ielprop,prop,
1112     &         ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,ipobody,ibody,
1113     &         xbodyact,mi,nacteq,bc,qat,nalt)
1114!
1115        elseif(lakon(nelem)(2:2).eq.'U') then
1116          call calcheatnet(nelem,lakon,ipkon,kon,v,ielprop,prop,
1117     &         ielmat,ntmat_,shcon,nshcon,rhcon,nrhcon,ipobody,ibody,
1118     &         xbodyact,mi,nacteq,bc,qat,nalt)
1119        else
1120          cycle
1121        endif
1122      enddo
1123!
1124!     transfer element ABSOLUTE TO RELATIVE / RELATIVE TO ABSOLUTE
1125!
1126      do i=1,nflow
1127        nelem=ieg(i)
1128!
1129        if((lakon(nelem)(2:4).eq.'ATR').or.
1130     &       (lakon(nelem)(2:4).eq.'RTA'))  then
1131!
1132          nodem=kon(ipkon(nelem)+2)
1133          xflow=v(1,nodem)
1134          if(xflow.gt.0d0) then
1135            node1=kon(ipkon(nelem)+1)
1136            node2=kon(ipkon(nelem)+3)
1137          else
1138            node1=kon(ipkon(nelem)+1)
1139            node2=kon(ipkon(nelem)+3)
1140          endif
1141!
1142!     computing temperature corrected Cp=Cp(T) coefficient
1143!
1144          Tg1=v(0,node1)
1145          Tg2=v(0,node2)
1146          if((lakon(nelem)(2:3).ne.'LP').and.
1147     &         (lakon(nelem)(2:3).ne.'LI')) then
1148            gastemp=(tg1+tg2)/2.d0
1149          else
1150            if(xflow.gt.0) then
1151              gastemp=tg1
1152            else
1153              gastemp=tg2
1154            endif
1155          endif
1156!
1157          imat=ielmat(1,nelem)
1158          call materialdata_tg(imat,ntmat_,gastemp,
1159     &         shcon,nshcon,cp,r,dvi,rhcon,nrhcon,rho)
1160!
1161          call cp_corrected(cp,Tg1,Tg2,cp_cor)
1162!
1163          index=ielprop(nelem)
1164          U=prop(index+1)
1165          ct=prop(index+2)
1166!
1167!     if a swirl element was given by the user, this takes
1168!     precendence to a value of ct
1169!
1170          if(nint(prop(index+3)).ne.0) then
1171            nelemswirl=nint(prop(index+3))
1172            index2=ielprop(nelemswirl)
1173!
1174!     previous element is a preswirl nozzle
1175!
1176            if(lakon(nelemswirl)(2:5).eq.'ORPN') then
1177              ct=prop(index2+4)
1178!
1179!     previous element is a forced vortex
1180!
1181            elseif(lakon(nelemswirl)(2:5).eq.'VOFO') then
1182              ct=prop(index2+7)
1183!
1184!     previous element is a free vortex
1185!
1186            elseif(lakon(nelemswirl)(2:5).eq.'VOFR') then
1187              ct=prop(index2+9)
1188            endif
1189          endif
1190!
1191          if(lakon(nelem)(2:4).eq.'ATR') then
1192            heat=Cp/Cp_cor*(0.5d0*(U**2-2d0*U*Ct)*xflow)
1193!
1194          elseif(lakon(nelem)(2:4).eq.'RTA') then
1195            heat=Cp/Cp_cor*(-0.5d0*(U**2-2d0*U*Ct)*xflow)
1196          endif
1197!
1198!     including the resulting additional heat flux in the energy equation
1199!
1200          if(xflow.gt.0d0)then
1201            ieq=nacteq(0,node2)
1202            if(ieq.ne.0) then
1203              bc(ieq)=bc(ieq)+heat
1204              qat=qat+dabs(heat)
1205              nalt=nalt+1
1206            endif
1207          else
1208            ieq=nacteq(0,node1)
1209            if(ieq.ne.0) then
1210              bc(ieq)=bc(ieq)+heat
1211              qat=qat+dabs(heat)
1212              nalt=nalt+1
1213            endif
1214          endif
1215        endif
1216      enddo
1217!
1218!     additional multiple point constraints
1219!
1220      j=nteq+1
1221      do i=nmpc,1,-1
1222        if(labmpc(i)(1:7).ne.'NETWORK') cycle
1223        j=j-1
1224        if(labmpc(i)(8:8).eq.' ') then
1225!
1226!     linear equation
1227!
1228          index=ipompc(i)
1229!
1230          do
1231            node=nodempc(1,index)
1232            idir=nodempc(2,index)
1233            bc(j)=bc(j)-v(idir,node)*coefmpc(index)
1234            index=nodempc(3,index)
1235            if(index.eq.0) exit
1236          enddo
1237        else
1238!
1239!     user-defined network equation
1240!
1241          call networkmpc_rhs(i,ipompc,nodempc,coefmpc,labmpc,
1242     &         v,bc,j,mi,ipkon,kon,lakon,iponoel,
1243     &         inoel,ielprop,prop,ielmat,
1244     &         shcon,nshcon,rhcon,nrhcon,ntmat_,cocon,ncocon)
1245        endif
1246      enddo
1247!
1248!     determining the typical energy flow
1249!
1250      if(nalt.gt.0) qat=qat/nalt
1251!
1252!     determining the typical mass flow
1253!
1254      if(nalf.gt.0) qaf=qaf/nalf
1255!
1256!     max. residual energy flow, residual mass flow
1257!     and residuals from the element equations
1258!
1259      ramt=0.d0
1260      ramf=0.d0
1261      ramp=0.d0
1262      do i=1,ntg
1263        node=itg(i)
1264        if(nacteq(0,node).ne.0) then
1265          ramt=max(ramt,dabs(bc(nacteq(0,node))))
1266        endif
1267        if(nacteq(1,node).ne.0) then
1268          ramf=max(ramf,dabs(bc(nacteq(1,node))))
1269        endif
1270        if(nacteq(2,node).ne.0) then
1271          ramp=max(ramp,dabs(bc(nacteq(2,node))))
1272        endif
1273      enddo
1274!
1275c     write(*,*) 'bc in resultnet'
1276c     do i=1,3
1277c     write(*,'(1x,e11.4)') bc(i)
1278c     enddo
1279!
1280      return
1281      end
1282