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 envtemp(itg,ieg,ntg,ntr,sideload,nelemload,
20     &     ipkon,kon,lakon,ielmat,ne,nload,kontri,ntri,nloadtr,
21     &     nflow,ndirboun,nactdog,nodeboun,nacteq,nboun,
22     &     ielprop,prop,nteq,v,network,physcon,shcon,ntmat_,
23     &     co,vold,set,nshcon,rhcon,nrhcon,mi,
24     &     nmpc,nodempc,ipompc,labmpc,ikboun,nasym,ttime,time,iaxial)
25!
26!     determines the number of gas temperatures and radiation
27!     temperatures
28!
29      implicit none
30!
31      logical identity,walltemp,temperaturebc,pressurebc,massflowbcall,
32     &     pressurebcall
33!
34      character*8 lakon(*)
35      character*20 labmpc(*)
36      character*20 sideload(*)
37      character*81 set(*)
38!
39      integer itg(*),ntg,ntr,nelemload(2,*),ipkon(*),network,mi(*),
40     &     kon(*),ielmat(mi(3),*),ne,i,j,k,l,index,id,node,nload,
41     &     ifaceq(8,6),ider,nasym,indexe,iaxial,networkmpcs,
42     &     ifacet(6,4),ifacew(8,5),kontri3(3,1),kontri4(3,2),
43     &     kontri6(3,4),kontri8(3,6),kontri(4,*),ntri,
44     &     konf(8),nloadtr(*),nelem,nope,nopes,ig,nflow,ieg(*),
45     &     ndirboun(*),nactdog(0:3,*),nboun,nodeboun(*),ntmat_,
46     &     idir,ntq,nteq,nacteq(0:3,*),node1,node2,nodem,
47     &     ielprop(*),idirf(8),iflag,imat,numf,nrhcon(*),nshcon(*),
48     &     nmpc,nodempc(3,*),ipompc(*),ikboun(*),iplausi
49!
50      real*8 prop(*),f,xflow,nodef(8),df(8),v(0:mi(2),*),g(3),
51     &     cp,r,physcon(*),shcon(0:3,ntmat_,*),rho,ttime,time,
52     &     co(3,*),dvi,vold(0:mi(2),*),rhcon(*)
53!
54      data ifaceq /4,3,2,1,11,10,9,12,
55     &            5,6,7,8,13,14,15,16,
56     &            1,2,6,5,9,18,13,17,
57     &            2,3,7,6,10,19,14,18,
58     &            3,4,8,7,11,20,15,19,
59     &            4,1,5,8,12,17,16,20/
60      data ifacet /1,3,2,7,6,5,
61     &     1,2,4,5,9,8,
62     &     2,3,4,6,10,9,
63     &     1,4,3,8,10,7/
64      data ifacew /1,3,2,9,8,7,0,0,
65     &     4,5,6,10,11,12,0,0,
66     &     1,2,5,4,7,14,10,13,
67     &     2,3,6,5,8,15,11,14,
68     &     4,6,3,1,12,15,9,13/
69      data kontri3 /1,2,3/
70      data kontri4 /1,2,4,2,3,4/
71      data kontri6 /1,4,6,4,5,6,4,2,5,6,5,3/
72      data kontri8 /1,5,8,8,5,7,8,7,4,5,2,6,5,6,7,7,6,3/
73!
74      ntg=0
75      ntr=0
76      ntri=0
77!
78      walltemp=.false.
79      temperaturebc=.false.
80      pressurebc=.false.
81      massflowbcall=.true.
82      pressurebcall=.true.
83!
84      networkmpcs=0
85!
86!     ordering the gas temperature nodes and counting them
87!     counting the radiation temperatures
88!
89      do i=1,nload
90         if(sideload(i)(3:4).eq.'FC') then
91            walltemp=.true.
92            call nident(itg,nelemload(2,i),ntg,id)
93            if(id.gt.0) then
94               if(itg(id).eq.nelemload(2,i)) then
95                  nactdog(0,nelemload(2,i))=1
96                  cycle
97               endif
98            endif
99            ntg=ntg+1
100            do j=ntg,id+2,-1
101               itg(j)=itg(j-1)
102            enddo
103            itg(id+1)=nelemload(2,i)
104            nactdog(0,nelemload(2,i))=1
105!
106         elseif(sideload(i)(3:4).eq.'NP') then
107            call nident(itg,nelemload(2,i),ntg,id)
108            if(id.gt.0) then
109               if(itg(id).eq.nelemload(2,i)) then
110                  nactdog(2,nelemload(2,i))=1
111                  cycle
112               endif
113            endif
114            ntg=ntg+1
115            do j=ntg,id+2,-1
116               itg(j)=itg(j-1)
117            enddo
118            itg(id+1)=nelemload(2,i)
119            nactdog(2,nelemload(2,i))=1
120!
121         elseif(sideload(i)(3:4).eq.'CR') then
122            ntr=ntr+1
123            nelem=nelemload(1,i)
124            read(sideload(i)(2:2),'(i1)') ig
125!
126!     number of nodes in the face
127!
128            if(lakon(nelem)(4:4).eq.'2') then
129               nope=20
130               nopes=8
131            elseif(lakon(nelem)(4:4).eq.'8') then
132               nope=8
133               nopes=4
134            elseif(lakon(nelem)(4:5).eq.'10') then
135               nope=10
136               nopes=6
137            elseif(lakon(nelem)(4:4).eq.'4') then
138               nope=4
139               nopes=3
140            elseif(lakon(nelem)(4:4).eq.'6') then
141               nope=6
142               if(ig.le.2) then
143                  nopes=3
144               else
145                  nopes=4
146               endif
147            elseif(lakon(nelem)(4:5).eq.'15') then
148               nope=15
149               if(ig.le.2) then
150                  nopes=6
151               else
152                  nopes=8
153               endif
154            endif
155!
156!     nodes in the face
157!
158            if((nope.eq.20).or.(nope.eq.8)) then
159               do k=1,nopes
160                  konf(k)=kon(ipkon(nelem)+ifaceq(k,ig))
161               enddo
162            elseif((nope.eq.10).or.(nope.eq.4)) then
163               do k=1,nopes
164                  konf(k)=kon(ipkon(nelem)+ifacet(k,ig))
165               enddo
166            else
167               do k=1,nopes
168                  konf(k)=kon(ipkon(nelem)+ifacew(k,ig))
169               enddo
170            endif
171!
172!     triangulation of the face
173!
174            nloadtr(ntr)=i
175            if((lakon(nelem)(4:4).eq.'2').or.
176     &           ((lakon(nelem)(4:5).eq.'15').and.(ig.gt.2))) then
177               do k=1,6
178                  ntri=ntri+1
179                  do l=1,3
180                     kontri(l,ntri)=konf(kontri8(l,k))
181                  enddo
182                  kontri(4,ntri)=ntr
183               enddo
184            elseif((lakon(nelem)(4:4).eq.'8').or.
185     &              ((lakon(nelem)(4:4).eq.'6').and.(ig.gt.2))) then
186               do k=1,2
187                  ntri=ntri+1
188                  do l=1,3
189                     kontri(l,ntri)=konf(kontri4(l,k))
190                  enddo
191                  kontri(4,ntri)=ntr
192               enddo
193            elseif((lakon(nelem)(4:5).eq.'10').or.
194     &              ((lakon(nelem)(4:5).eq.'15').and.(ig.le.2))) then
195               do k=1,4
196                  ntri=ntri+1
197                  do l=1,3
198                     kontri(l,ntri)=konf(kontri6(l,k))
199                  enddo
200                  kontri(4,ntri)=ntr
201               enddo
202            elseif((lakon(nelem)(4:4).eq.'4').or.
203     &              ((lakon(nelem)(4:4).eq.'6').and.(ig.le.2))) then
204               do k=1,1
205                  ntri=ntri+1
206                  do l=1,3
207                     kontri(l,ntri)=konf(kontri3(l,k))
208                  enddo
209                  kontri(4,ntri)=ntr
210               enddo
211            endif
212         endif
213      enddo
214!
215!     storing the gas elements in a dedicated array
216!
217      nflow=0
218!
219      do i=1,ne
220         if(lakon(i)(1:1).eq.'D') then
221            if((lakon(i)(2:2).ne.' ').and.(network.ne.1)) then
222               nflow=nflow+1
223               ieg(nflow)=i
224            else
225               nasym=1
226!
227!              removing gas nodes belonging to 'D '-elements
228!              in which a 'FC'-film condition was applied from the
229!              itg vector
230!
231               indexe=ipkon(i)
232               do j=1,3,2
233                  node=kon(indexe+j)
234                  call nident(itg,node,ntg,id)
235                  if(id.gt.0) then
236                     if(itg(id).eq.node) then
237                        ntg=ntg-1
238                        do k=id,ntg
239                           itg(k)=itg(k+1)
240                        enddo
241                        nactdog(0,node)=0
242                        nactdog(2,node)=0
243                     endif
244                  endif
245               enddo
246!
247            endif
248         endif
249!
250!        removing gas nodes belonging to advective elements
251!        (last node in the element topology)
252!        these are usually gas nodes not belonging to any
253!        "D"-element
254!
255         if(lakon(i)(1:7).eq.'ESPRNGF') then
256            read(lakon(i)(8:8),'(i1)') nopes
257            nope=nopes+1
258            node=kon(ipkon(i)+nope)
259!
260            call nident(itg,node,ntg,id)
261            if(id.gt.0) then
262               if(itg(id).eq.node) then
263                  ntg=ntg-1
264                  do k=id,ntg
265                     itg(k)=itg(k+1)
266                  enddo
267                  nactdog(0,node)=0
268                  nactdog(2,node)=0
269               endif
270            endif
271         endif
272      enddo
273!
274!     mass flux nodes are also taken as unknowns in the
275!     gas temperature system; determining the active
276!     degrees of freedom
277!
278!     first node of the flow element
279!
280      do i=1,nflow
281         index=ipkon(ieg(i))
282         node=kon(index+1)
283         if (node.eq.0) cycle
284         call nident(itg,node,ntg,id)
285         if(id.gt.0) then
286            if(itg(id).eq.node) then
287!
288!              upstream depth of SO,WO and DO is known
289!
290               if((lakon(ieg(i))(4:7).eq.'CHSO').or.
291     &              (lakon(ieg(i))(4:7).eq.'CHWO').or.
292     &              (lakon(ieg(i))(4:7).eq.'CHDO')) cycle
293               nactdog(0,node)=1
294               nactdog(2,node)=1
295               cycle
296           endif
297         endif
298         ntg=ntg+1
299         do j=ntg,id+2,-1
300            itg(j)=itg(j-1)
301         enddo
302         itg(id+1)=node
303!
304!        upstream depth of SO,WO and DO is known
305!
306         if((lakon(ieg(i))(4:7).eq.'CHSO').or.
307     &      (lakon(ieg(i))(4:7).eq.'CHWO').or.
308     &      (lakon(ieg(i))(4:7).eq.'CHDO')) cycle
309         nactdog(0,node)=1
310         nactdog(2,node)=1
311      enddo
312!
313!     middle node of the flow element :flux
314!
315      do i=1,nflow
316         index=ipkon(ieg(i))
317         node=kon(index+2)
318         call nident(itg,node,ntg,id)
319         if(id.gt.0) then
320            if(itg(id).eq.node) cycle
321         endif
322         ntg=ntg+1
323         do j=ntg,id+2,-1
324            itg(j)=itg(j-1)
325        enddo
326        itg(id+1)=node
327        nactdog(1,node)=1
328!
329!        variable geometric property
330!
331        if(lakon(ieg(i))(6:7).eq.'GV') then
332           index=ielprop(ieg(i))
333           if(prop(index+2).le.0.d0) nactdog(3,node)=1
334        elseif((lakon(ieg(i))(4:7).eq.'CHSG').or.
335     &         (lakon(ieg(i))(4:7).eq.'CHWE').or.
336     &         (lakon(ieg(i))(4:7).eq.'CHDS')) then
337           nactdog(3,node)=1
338        elseif(lakon(ieg(i))(2:7).eq.'ACCTUB') then
339!
340            index=ielprop(ieg(i))
341            if(prop(index+1).eq.2) then
342!              Interval factor unknown,setting a DOF for geometry
343               nactdog(3,node)=1
344            elseif(prop(index+1).eq.3) then
345!              Hole diameter factor unknown,setting a DOF for geometry
346               nactdog(3,node)=1
347            endif
348        endif
349      enddo
350!
351!     third node of the flow element
352!
353      do i=1,nflow
354         index=ipkon(ieg(i))
355         node=kon(index+3)
356         if (node.eq.0) cycle
357         call nident(itg,node,ntg,id)
358         if(id.gt.0) then
359            if(itg(id).eq.node) then
360!
361!              downstream depth of SG,WE and DS is known
362!
363               if((lakon(ieg(i))(4:7).eq.'CHSG').or.
364     &            (lakon(ieg(i))(4:7).eq.'CHWE').or.
365     &            (lakon(ieg(i))(4:7).eq.'CHDS')) cycle
366               nactdog(0,node)=1
367               nactdog(2,node)=1
368               cycle
369            endif
370         endif
371         ntg=ntg+1
372         do j=ntg,id+2,-1
373            itg(j)=itg(j-1)
374         enddo
375         itg(id+1)=node
376!
377!        downstream depth of SG,WE and DS is known
378!
379         if((lakon(ieg(i))(4:7).eq.'CHSG').or.
380     &        (lakon(ieg(i))(4:7).eq.'CHWE').or.
381     &        (lakon(ieg(i))(4:7).eq.'CHDS')) cycle
382         nactdog(0,node)=1
383         nactdog(2,node)=1
384      enddo
385!
386!     tagging the network MPC's
387!
388      do i=1,nmpc
389!
390!        check whether network MPC
391!
392         index=ipompc(i)
393         do
394            node=nodempc(1,index)
395            call nident(itg,node,ntg,id)
396            if(id.gt.0) then
397               if(itg(id).eq.node) then
398                  labmpc(i)(1:7)='NETWORK'
399                  networkmpcs=1
400                  exit
401               endif
402            endif
403            index=nodempc(3,index)
404            if(index.eq.0) exit
405         enddo
406      enddo
407!
408!     taking the MPC network nodes into account
409!
410      if(networkmpcs.eq.1) then
411         do i=1,nmpc
412            if(labmpc(i)(1:7).ne.'NETWORK') cycle
413!
414            index=ipompc(i)
415            do
416               node=nodempc(1,index)
417               call nident(itg,node,ntg,id)
418               if(id.gt.0) then
419                  if(itg(id).eq.node) then
420                     nactdog(nodempc(2,index),node)=1
421                     index=nodempc(3,index)
422                     if(index.eq.0) then
423                        exit
424                     else
425                        cycle
426                     endif
427                  endif
428               endif
429!
430!              adding a node to itg
431!
432               ntg=ntg+1
433               do j=ntg,id+2,-1
434                  itg(j)=itg(j-1)
435               enddo
436               itg(id+1)=node
437               nactdog(nodempc(2,index),node)=1
438               index=nodempc(3,index)
439               if(index.eq.0) exit
440            enddo
441         enddo
442      endif
443!
444!     subtracting the SPC conditions
445!
446      do i=1,nboun
447         node=nodeboun(i)
448         call nident(itg,node,ntg,id)
449         if (id.gt.0) then
450            if (itg(id).eq.node) then
451               idir=ndirboun(i)
452               nactdog(idir,node)=0
453               if(idir.eq.0) then
454                  temperaturebc=.true.
455               elseif(idir.eq.2) then
456                  pressurebc=.true.
457               endif
458            endif
459         endif
460      enddo
461!
462!     temporarily removing the dependent nodes of the MPC's
463!     only for mass flow and pressure
464!
465!     these are the only dofs for which the corresponding equations
466!     (mass equilibrium in end node and element equation in middle
467!     node) are not applied in the nodes where the dofs are lacking
468!     (mass flow in middle node and pressure in end node)
469!
470      if(networkmpcs.eq.1) then
471         do i=1,nmpc
472            if(labmpc(i)(1:7).ne.'NETWORK') cycle
473            index=ipompc(i)
474            idir=nodempc(2,index)
475            if((idir.eq.1).or.(idir.eq.2)) then
476               node=nodempc(1,index)
477               call nident(itg,node,ntg,id)
478               if(id.gt.0) then
479                  if(itg(id).eq.node) then
480                     nactdog(idir,node)=0
481                  endif
482               endif
483            endif
484         enddo
485      endif
486!
487!     determining the active equations
488!
489!     element contributions
490!
491      iflag=0
492      do i=1,nflow
493         nelem=ieg(i)
494         index=ipkon(nelem)
495         node1=kon(index+1)
496         nodem=kon(index+2)
497         node2=kon(index+3)
498!
499!     "end of network" element ---X---O
500!                             1   m   2
501!
502         if(node1.eq.0)then
503            if ((nactdog(1,nodem).ne.0).or.(nactdog(3,nodem).ne.0))then
504               nacteq(1,node2)=1                      ! mass equation
505            endif
506            if (nactdog(0,node2).ne.0)then
507               nacteq(0,node2)=1                      ! energy equation
508            endif
509!
510!     "end of network" element node O---X---
511!     1   m   2
512         elseif (node2.eq.0) then
513            if ((nactdog(1,nodem).ne.0).or.(nactdog(3,nodem).ne.0))then
514               nacteq(1,node1)=1                      ! mass equation
515            endif
516            if (nactdog(0,node1).ne.0)then
517               nacteq(0,node1)=1                      ! energy equation
518            endif
519!
520!     "flow element" O---X---O
521!                    1   m   2
522!
523         else
524            if((nactdog(1,nodem).ne.0).or.(nactdog(3,nodem).ne.0)) then
525               nacteq(1,node2)=1                       ! mass equation
526               nacteq(1,node1)=1                       ! mass equation
527            endif
528!
529            ider=0
530            call flux(node1,node2,nodem,nelem,lakon,kon,ipkon,
531     &           nactdog,identity,ielprop,prop,iflag,v,xflow,f,
532     &           nodef,idirf,df,cp,r,rho,physcon,g,co,dvi,numf,
533     &           vold,set,shcon,nshcon,rhcon,nrhcon,ntmat_,mi,ider,
534     &           ttime,time,iaxial,iplausi)
535!
536            if (.not.identity) then
537               nacteq(2,nodem)=1                       ! momentum equation
538            endif
539!
540            if (nactdog(0,node1).ne.0)then
541               nacteq(0,node1)=1                       ! energy equation
542            endif
543!
544            if (nactdog(0,node2).ne.0)then
545               nacteq(0,node2)=1                       ! energy equation
546            endif
547         endif
548      enddo
549!
550!     wall convective contributions
551!
552      do i=1, nload
553         if((sideload(i)(3:4)).eq.'FC')then
554            node=nelemload(2,i)
555            if (nactdog(0,node).ne.0) then
556               nacteq(0,node)=1
557            endif
558         endif
559      enddo
560!
561!     restoring the dependent nodes of the MPC's
562!     only for mass flow and pressure
563!
564      if(networkmpcs.eq.1) then
565         do i=1,nmpc
566            if(labmpc(i)(1:7).ne.'NETWORK') cycle
567            index=ipompc(i)
568            idir=nodempc(2,index)
569            if((idir.eq.1).or.(idir.eq.2)) then
570               node=nodempc(1,index)
571               call nident(itg,node,ntg,id)
572               if(id.gt.0) then
573                  if(itg(id).eq.node) then
574                     nactdog(idir,node)=1
575                  endif
576               endif
577            endif
578         enddo
579      endif
580!
581!     removing the energy equation from those end nodes for which
582!     the temperature constitutes the first term in a network MPC
583!
584      if(networkmpcs.eq.1) then
585         do i=1,nmpc
586            if(labmpc(i)(1:7).ne.'NETWORK') cycle
587            index=ipompc(i)
588            idir=nodempc(2,index)
589            if(idir.eq.0) then
590               node=nodempc(1,index)
591               call nident(itg,node,ntg,id)
592               if(id.gt.0) then
593                  if(itg(id).eq.node) then
594                     nacteq(0,node)=0
595                  endif
596               endif
597            endif
598         enddo
599      endif
600!
601!     check whether all mass flow is known
602!
603      do i=1,ntg
604         node=itg(i)
605         if((nactdog(1,node).ne.0).or.(nactdog(3,node).ne.0)) then
606            massflowbcall=.false.
607            exit
608         endif
609      enddo
610!
611!     check whether all pressures are known
612!
613      do i=1,ntg
614         node=itg(i)
615         if(nactdog(2,node).ne.0) then
616            pressurebcall=.false.
617            exit
618         endif
619      enddo
620!
621!     check for special cases
622!
623      if(massflowbcall.and.((.not.pressurebc).or.(pressurebcall))) then
624!
625!        purely thermal (only set to 1 if D-type elements are present)
626!
627         if(network.gt.1) network=2
628         do i=1,nflow
629            nelem=ieg(i)
630            index=ipkon(nelem)
631            node1=kon(index+1)
632            nodem=kon(index+2)
633            node2=kon(index+3)
634            nacteq(2,nodem)=0
635            if(node1.ne.0) nactdog(2,node1)=0
636            if(node2.ne.0) nactdog(2,node2)=0
637         enddo
638      elseif((.not.temperaturebc).and.(.not.walltemp)) then
639!
640!        pure liquid dynamics
641!
642         write(*,*) '*INFO in envtemp: no thermal boundary conditions'
643         write(*,*) '      detected; the network is considered to be'
644         write(*,*) '      athermal and no gas temperatures will be'
645         write(*,*) '      calculated'
646         network=4
647         do i=1,ntg
648            node=itg(i)
649            nactdog(0,node)=0
650            nacteq(0,node)=0
651         enddo
652      elseif((.not.temperaturebc).and.walltemp) then
653         write(*,*) '*ERROR in envtemp: at least one temperature'
654         write(*,*) '       boundary condition must be given'
655         call exit(201)
656      elseif(.not.pressurebc) then
657         write(*,*) '*ERROR in envtemp: at least one pressure'
658         write(*,*) '       boundary condition must be given'
659         call exit(201)
660      endif
661!
662!     check whether a specific gas constant was defined for all fluid
663!     elements (except for purely thermal calculations)
664!
665      if(network.gt.2) then
666         do i=1,nflow
667            nelem=ieg(i)
668            if((lakon(nelem)(2:3).eq.'LI').or.
669     &         (lakon(nelem)(2:3).eq.'LP').or.
670     &         (lakon(nelem)(2:3).eq.'  ')) cycle
671            imat=ielmat(1,nelem)
672            r=shcon(3,1,imat)
673            if(r.lt.1.d-10) then
674               write(*,*)'*ERROR in envtemp: specific gas',
675     &              'constant is close to zero'
676               call exit(201)
677            endif
678         enddo
679      endif
680!
681!     check whether the temperature at each inlet or outlet node
682!     is given
683!
684      if(network.lt.4) then
685         do i=1,nflow
686            nelem=ieg(i)
687            index=ipkon(nelem)
688            node1=kon(index+1)
689            node2=kon(index+3)
690            if(node1.eq.0) then
691               if(nactdog(0,node2).ne.0) then
692               write(*,*) '*WARNING in envtemp: it is advised to'
693               write(*,*) '         define the temperature at all'
694               write(*,*) '         inlets and outlets by a boundary'
695               write(*,*) '         condition. This is lacking for'
696               write(*,*) '         node ',node2,' of element ',nelem
697               endif
698            endif
699            if(node2.eq.0) then
700               if(nactdog(0,node1).ne.0) then
701               write(*,*) '*WARNING in envtemp: it is advised to'
702               write(*,*) '         define the temperature at all'
703               write(*,*) '         inlets and outlets by a boundary'
704               write(*,*) '         condition. This is lacking for'
705               write(*,*) '         node ',node1,' of element ',nelem
706               endif
707            endif
708         enddo
709      endif
710!
711!     numbering the active equations
712!
713      nteq=0
714      do i=1,ntg
715         node=itg(i)
716         do j=0,2
717            if (nacteq(j,node).ne.0) then
718               nteq=nteq+1
719               nacteq(j,node)=nteq
720            endif
721         enddo
722!         write(30,*) 'unknowns ',node,(nactdog(j,node),j=0,3)
723      enddo
724!      do i=1,ntg
725!         node=itg(i)
726!         write(30,*) 'equations',node,(nacteq(j,node),j=0,2)
727!      enddo
728!
729!     taking network MPC's into account
730!
731      if(networkmpcs.eq.1) then
732         do i=1,nmpc
733            if(labmpc(i)(1:7).eq.'NETWORK') nteq=nteq+1
734         enddo
735      endif
736!
737!     numbering the active degrees of freedom
738!
739      ntq=0
740      do i=1,ntg
741         node=itg(i)
742         do  j=0,3
743            if (nactdog(j,node).ne.0) then
744               ntq=ntq+1
745               nactdog(j,node)=ntq
746            endif
747         enddo
748      enddo
749c
750c      open(30,file='dummy',status='unknown')
751c      write(30,*) 'nactdog'
752c      do i=1,ntg
753c         write(30,*) itg(i),(nactdog(j,itg(i)),j=0,3)
754c      enddo
755c
756c      write(30,*) ''
757c      write(30,*) 'nacteq'
758c      do i=1,ntg
759c         write(30,*) itg(i),(nacteq(j,itg(i)),j=0,3)
760c      enddo
761c      close(30)
762!
763      if(ntq.ne.nteq) then
764         write(*,*) '*ERROR in envtemp:'
765         write(*,*) '*****number of network equations is not equal to'
766         write(*,*) ' number of active degrees of freedom*****'
767         write(*,*) ' # of network equations = ',nteq
768         write(*,*) ' # of active degrees of freedom= ',ntq
769         call exit(201)
770      endif
771!
772!     for isothermal gas pipes the energy equation in the
773!     topologically downstream node is replaced by an equation
774!     expressing the equality of the static temperature at both
775!     ends of the pipe. To this end these downstream nodes are
776!     referring in nacteq(3,*) to the topologically upstream node
777!
778!     if the temperature in the downstream node is a boundary
779!     condition (i.e. the energy equation is not built), the
780!     energy equation in the upstream node is replaced. In that
781!     case nacteq(3,upstreamnode) refers to the downstream node.
782!
783!     if both nodes are boundary conditions, nothing is done
784!
785      do i=1,nflow
786         nelem=ieg(i)
787         if((lakon(nelem)(1:4).eq."DGAP")
788     &        .and.(lakon(nelem)(6:6).eq."I")) then
789!
790            index=ipkon(nelem)
791            node1=kon(index+1)
792            node2=kon(index+3)
793            if((node1.eq.0).or.(node2.eq.0)) cycle
794!
795            if(nacteq(0,node2).ne.0) then
796               nacteq(3,node2)=node1
797            elseif(nacteq(0,node1).ne.0) then
798               nacteq(3,node1)=node2
799            endif
800         endif
801      enddo
802!
803      return
804      end
805
806
807