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 initialconditionss(inpc,textpart,set,istartset,iendset,
20     &     ialset,nset,t0,t1,prestr,iprestr,ithermal,veold,inoelfree,
21     &     nk_,
22     &     mi,istep,istat,n,iline,ipol,inl,ipoinp,inp,lakon,kon,co,ne,
23     &     ipkon,vold,ipoinpc,xstate,nstate_,nk,t0g,t1g,iaxial,ielprop,
24     &     prop,ier,nuel_)
25!
26!     reading the input deck: *INITIAL CONDITIONS
27!
28      implicit none
29!
30      logical user
31!
32      character*1 inpc(*)
33      character*8 lakon(*)
34      character*80 rebarn
35      character*81 set(*),noset
36      character*132 textpart(16)
37!
38      integer istartset(*),iendset(*),ialset(*),nset,iprestr,
39     &     ithermal(*),m,id,nuel_,
40     &     istep,istat,n,i,j,k,l,ii,key,idir,ipos,inoelfree,nk_,mi(*),
41     &     iline,ipol,inl,ipoinp(2,*),inp(3,*),ij,jj,ntens,ncrds,layer,
42     &     kspt,lrebar,iflag,i1,mint3d,nope,kon(*),konl(20),indexe,
43     &     ipkon(*),ne,ipoinpc(0:*),nstate_,nk,jmax,ntot,numberoflines,
44     &     iaxial,null,ielprop(*),ier
45!
46      real*8 t0(*),t1(*),beta(8),prestr(6,mi(1),*),veold(0:mi(2),*),
47     &     temperature,velocity,tempgrad1,tempgrad2,pgauss(3),
48     &     shp(4,20),xsj,xl(3,20),xi,et,ze,weight,co(3,*),pressure,
49     &     vold(0:mi(2),*),xstate(nstate_,mi(1),*),dispvelo,totpres,
50     &     xmassflow,t0g(2,*),t1g(2,*),prop(*),turbini(0:mi(2))
51!
52      include "gauss.f"
53!
54      null=0
55!
56      if(istep.gt.0) then
57        write(*,*)
58     &       '*ERROR reading *INITIAL CONDITIONS: *INITIAL CONDITIONS'
59        write(*,*) '  should be placed before all step definitions'
60        ier=1
61        return
62      endif
63!
64      do ij=2,n
65        if(textpart(ij)(1:16).eq.'TYPE=TEMPERATURE') then
66!
67          ithermal(1)=1
68          do
69            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
70     &           ipoinp,inp,ipoinpc)
71            if((istat.lt.0).or.(key.eq.1)) return
72            read(textpart(2)(1:20),'(f20.0)',iostat=istat)
73     &           temperature
74            if(istat.gt.0) then
75              call inputerror(inpc,ipoinpc,iline,
76     &             "*INITIAL CONDITIONS%",ier)
77              return
78            endif
79            temperature=1.d-6*int(1.d6*temperature+0.5d0)
80!
81            if((inoelfree.ne.0).or.(nuel_.gt.0)) then
82              tempgrad1=0.d0
83              tempgrad2=0.d0
84              if(n.gt.2) then
85                read(textpart(3)(1:20),'(f20.0)',iostat=istat)
86     &               tempgrad1
87                if(istat.gt.0) then
88                  call inputerror(inpc,ipoinpc,iline,
89     &                 "*INITIAL CONDITIONS%",ier)
90                  return
91                endif
92              endif
93              if(n.gt.3) then
94                read(textpart(4)(1:20),'(f20.0)',iostat=istat)
95     &               tempgrad2
96                if(istat.gt.0) then
97                  call inputerror(inpc,ipoinpc,iline,
98     &                 "*INITIAL CONDITIONS%",ier)
99                  return
100                endif
101              endif
102            endif
103!
104            read(textpart(1)(1:10),'(i10)',iostat=istat) l
105            if(istat.eq.0) then
106              if(l.gt.nk) then
107                write(*,*)
108     &               '*WARNING reading *INITIAL CONDITIONS: node ',l
109                write(*,*)'          exceeds the largest defined ',
110     &               'node number'
111                cycle
112              endif
113              t0(l)=temperature
114              t1(l)=temperature
115              vold(0,l)=temperature
116              if((inoelfree.ne.0).or.(nuel_.gt.0)) then
117                t0g(1,l)=tempgrad1
118                t0g(2,l)=tempgrad2
119                t1g(1,l)=tempgrad1
120                t1g(2,l)=tempgrad2
121              endif
122            else
123              read(textpart(1)(1:80),'(a80)',iostat=istat) noset
124              noset(81:81)=' '
125              ipos=index(noset,' ')
126              noset(ipos:ipos)='N'
127c              do ii=1,nset
128c                if(set(ii).eq.noset) exit
129c              enddo
130              call cident81(set,noset,nset,id)
131              ii=nset+1
132              if(id.gt.0) then
133                if(noset.eq.set(id)) then
134                  ii=id
135                endif
136              endif
137              if(ii.gt.nset) then
138                noset(ipos:ipos)=' '
139                write(*,*)
140     &               '*ERROR reading *INITIAL CONDITIONS: node set '
141     &               ,noset
142                write(*,*)'  has not yet been defined. '
143                call inputerror(inpc,ipoinpc,iline,
144     &               "*INITIAL CONDITIONS%",ier)
145                return
146              endif
147              do j=istartset(ii),iendset(ii)
148                if(ialset(j).gt.0) then
149                  t0(ialset(j))=temperature
150                  t1(ialset(j))=temperature
151                  vold(0,ialset(j))=temperature
152                  if((inoelfree.ne.0).or.(nuel_.gt.0)) then
153                    t0g(1,ialset(j))=tempgrad1
154                    t0g(2,ialset(j))=tempgrad2
155                    t1g(1,ialset(j))=tempgrad1
156                    t1g(2,ialset(j))=tempgrad2
157                  endif
158                else
159                  k=ialset(j-2)
160                  do
161                    k=k-ialset(j)
162                    if(k.ge.ialset(j-1)) exit
163                    t0(k)=temperature
164                    t1(k)=temperature
165                    vold(0,k)=temperature
166                    if((inoelfree.ne.0).or.(nuel_.gt.0)) then
167                      t0g(1,k)=tempgrad1
168                      t0g(2,k)=tempgrad2
169                      t1g(1,k)=tempgrad1
170                      t1g(2,k)=tempgrad2
171                    endif
172                  enddo
173                endif
174              enddo
175            endif
176          enddo
177          return
178        elseif(textpart(ij)(1:11).eq.'TYPE=STRESS') then
179!
180          iprestr=1
181          do jj=1,n
182            if(textpart(jj)(1:4).eq.'USER') then
183!
184!     residual stresses are defined by user subroutine
185!     sigini
186!
187              iflag=1
188              ntens=6
189              ncrds=3
190              lrebar=0
191              do i=1,ne
192                indexe=ipkon(i)
193                if(lakon(i)(4:4).eq.'2') then
194                  nope=20
195                elseif(lakon(i)(4:4).eq.'8') then
196                  nope=8
197                elseif(lakon(i)(4:5).eq.'10') then
198                  nope=10
199                elseif(lakon(i)(4:4).eq.'4') then
200                  nope=4
201                elseif(lakon(i)(4:5).eq.'15') then
202                  nope=15
203                elseif(lakon(i)(4:4).eq.'6') then
204                  nope=6
205                else
206                  cycle
207                endif
208!
209                if(lakon(i)(4:5).eq.'8R') then
210                  mint3d=1
211                elseif(lakon(i)(4:7).eq.'20RB') then
212                  if((lakon(i)(8:8).eq.'R').or.
213     &                 (lakon(i)(8:8).eq.'C')) then
214                    mint3d=50
215                  else
216                    call beamintscheme(lakon(i),mint3d,
217     &                   ielprop(i),prop,
218     &                   null,xi,et,ze,weight)
219                  endif
220                elseif((lakon(i)(4:4).eq.'8').or.
221     &                 (lakon(i)(4:6).eq.'20R')) then
222                  mint3d=8
223                elseif(lakon(i)(4:4).eq.'2') then
224                  mint3d=27
225                elseif(lakon(i)(4:5).eq.'10') then
226                  mint3d=4
227                elseif(lakon(i)(4:4).eq.'4') then
228                  mint3d=1
229                elseif(lakon(i)(4:5).eq.'15') then
230                  mint3d=9
231                elseif(lakon(i)(4:4).eq.'6') then
232                  mint3d=2
233                endif
234!
235                do j=1,nope
236                  konl(j)=kon(indexe+j)
237                  do k=1,3
238                    xl(k,j)=co(k,konl(j))
239                  enddo
240                enddo
241!
242                do j=1,mint3d
243                  if(lakon(i)(4:5).eq.'8R') then
244                    xi=gauss3d1(1,j)
245                    et=gauss3d1(2,j)
246                    ze=gauss3d1(3,j)
247                    weight=weight3d1(j)
248                  elseif(lakon(i)(4:7).eq.'20RB') then
249                    if((lakon(i)(8:8).eq.'R').or.
250     &                   (lakon(i)(8:8).eq.'C')) then
251                      xi=gauss3d13(1,j)
252                      et=gauss3d13(2,j)
253                      ze=gauss3d13(3,j)
254                      weight=weight3d13(j)
255                    else
256                      call beamintscheme(lakon(i),mint3d,
257     &                     ielprop(i),prop,
258     &                     j,xi,et,ze,weight)
259                    endif
260                  elseif((lakon(i)(4:4).eq.'8').or.
261     &                   (lakon(i)(4:6).eq.'20R'))
262     &                   then
263                    xi=gauss3d2(1,j)
264                    et=gauss3d2(2,j)
265                    ze=gauss3d2(3,j)
266                    weight=weight3d2(j)
267                  elseif(lakon(i)(4:4).eq.'2') then
268                    xi=gauss3d3(1,j)
269                    et=gauss3d3(2,j)
270                    ze=gauss3d3(3,j)
271                    weight=weight3d3(j)
272                  elseif(lakon(i)(4:5).eq.'10') then
273                    xi=gauss3d5(1,j)
274                    et=gauss3d5(2,j)
275                    ze=gauss3d5(3,j)
276                    weight=weight3d5(j)
277                  elseif(lakon(i)(4:4).eq.'4') then
278                    xi=gauss3d4(1,j)
279                    et=gauss3d4(2,j)
280                    ze=gauss3d4(3,j)
281                    weight=weight3d4(j)
282                  elseif(lakon(i)(4:5).eq.'15') then
283                    xi=gauss3d8(1,j)
284                    et=gauss3d8(2,j)
285                    ze=gauss3d8(3,j)
286                    weight=weight3d8(j)
287                  elseif(lakon(i)(4:4).eq.'6') then
288                    xi=gauss3d7(1,j)
289                    et=gauss3d7(2,j)
290                    ze=gauss3d7(3,j)
291                    weight=weight3d7(j)
292                  endif
293!
294                  if(nope.eq.20) then
295                    call shape20h(xi,et,ze,xl,xsj,shp,iflag)
296                  elseif(nope.eq.8) then
297                    call shape8h(xi,et,ze,xl,xsj,shp,iflag)
298                  elseif(nope.eq.10) then
299                    call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
300                  elseif(nope.eq.4) then
301                    call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
302                  elseif(nope.eq.15) then
303                    call shape15w(xi,et,ze,xl,xsj,shp,iflag)
304                  else
305                    call shape6w(xi,et,ze,xl,xsj,shp,iflag)
306                  endif
307!
308                  do k=1,3
309                    pgauss(k)=0.d0
310                    do i1=1,nope
311                      pgauss(k)=pgauss(k)+
312     &                     shp(4,i1)*co(k,konl(i1))
313                    enddo
314                  enddo
315!
316                  call sigini(prestr(1,j,i),pgauss,ntens,ncrds,
317     &                 i,j,layer,kspt,lrebar,rebarn)
318!
319                enddo
320              enddo
321              call getnewline(inpc,textpart,istat,n,key,iline,ipol,
322     &             inl,ipoinp,inp,ipoinpc)
323              return
324            endif
325          enddo
326!
327!     residual stresses are written explicitly in the input deck
328!
329          do
330            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
331     &           ipoinp,inp,ipoinpc)
332            if((istat.lt.0).or.(key.eq.1)) return
333            do j=1,6
334              read(textpart(j+2)(1:20),'(f20.0)',iostat=istat)
335     &             beta(j)
336              if(istat.gt.0) then
337                call inputerror(inpc,ipoinpc,iline,
338     &               "*INITIAL CONDITIONS%",ier)
339                return
340              endif
341            enddo
342            read(textpart(1)(1:10),'(i10)',iostat=istat) l
343            if(istat.ne.0) then
344              call inputerror(inpc,ipoinpc,iline,
345     &             "*INITIAL CONDITIONS%",ier)
346              return
347            endif
348            if(l.gt.ne) then
349              write(*,*)
350     &             '*WARNING reading *INITIAL CONDITIONS: element ',l
351              write(*,*)'          exceeds the largest defined ',
352     &             'element number'
353              cycle
354            endif
355            read(textpart(2)(1:10),'(i10)',iostat=istat) k
356            if(istat.eq.0) then
357              do j=1,6
358                prestr(j,k,l)=beta(j)
359              enddo
360            else
361              call inputerror(inpc,ipoinpc,iline,
362     &             "*INITIAL CONDITIONS%",ier)
363              return
364            endif
365          enddo
366          return
367        elseif(textpart(ij)(1:18).eq.'TYPE=PLASTICSTRAIN') then
368!
369          iprestr=2
370          do
371            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
372     &           ipoinp,inp,ipoinpc)
373            if((istat.lt.0).or.(key.eq.1)) return
374            do j=1,6
375              read(textpart(j+2)(1:20),'(f20.0)',iostat=istat)
376     &             beta(j)
377              if(istat.gt.0) then
378                call inputerror(inpc,ipoinpc,iline,
379     &               "*INITIAL CONDITIONS%",ier)
380                return
381              endif
382            enddo
383            read(textpart(1)(1:10),'(i10)',iostat=istat) l
384            if(istat.ne.0) then
385              call inputerror(inpc,ipoinpc,iline,
386     &             "*INITIAL CONDITIONS%",ier)
387              return
388            endif
389            if(l.gt.ne) then
390              write(*,*)
391     &             '*WARNING reading *INITIAL CONDITIONS: element ',l
392              write(*,*)'          exceeds the largest defined ',
393     &             'element number'
394              cycle
395            endif
396            read(textpart(2)(1:10),'(i10)',iostat=istat) k
397            if(istat.eq.0) then
398              do j=1,6
399                prestr(j,k,l)=beta(j)
400              enddo
401            else
402              call inputerror(inpc,ipoinpc,iline,
403     &             "*INITIAL CONDITIONS%",ier)
404              return
405            endif
406          enddo
407          return
408        elseif((textpart(ij)(1:17).eq.'TYPE=DISPLACEMENT').or.
409     &         (textpart(ij)(1:18).eq.'TYPE=FLUIDVELOCITY')) then
410!
411          do
412            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
413     &           ipoinp,inp,ipoinpc)
414            if((istat.lt.0).or.(key.eq.1)) return
415            read(textpart(2)(1:10),'(i10)',iostat=istat) idir
416            if(istat.gt.0) then
417              call inputerror(inpc,ipoinpc,iline,
418     &             "*INITIAL CONDITIONS%",ier)
419              return
420            endif
421            read(textpart(3)(1:20),'(f20.0)',iostat=istat) dispvelo
422            if(istat.gt.0) then
423              call inputerror(inpc,ipoinpc,iline,
424     &             "*INITIAL CONDITIONS%",ier)
425              return
426            endif
427            read(textpart(1)(1:10),'(i10)',iostat=istat) l
428            if(istat.eq.0) then
429              if(l.gt.nk) then
430                write(*,*)
431     &               '*WARNING reading *INITIAL CONDITIONS: node ',l
432                write(*,*)'          exceeds the largest defined ',
433     &               'node number'
434                cycle
435              endif
436              vold(idir,l)=dispvelo
437            else
438              read(textpart(1)(1:80),'(a80)',iostat=istat) noset
439              noset(81:81)=' '
440              ipos=index(noset,' ')
441              noset(ipos:ipos)='N'
442c              do ii=1,nset
443c                if(set(ii).eq.noset) exit
444c              enddo
445              call cident81(set,noset,nset,id)
446              ii=nset+1
447              if(id.gt.0) then
448                if(noset.eq.set(id)) then
449                  ii=id
450                endif
451              endif
452              if(ii.gt.nset) then
453                noset(ipos:ipos)=' '
454                write(*,*)
455     &               '*ERROR reading *INITIAL CONDITIONS: node set '
456     &               ,noset
457                write(*,*)'  has not yet been defined. '
458                call inputerror(inpc,ipoinpc,iline,
459     &               "*INITIAL CONDITIONS%",ier)
460                return
461              endif
462              do j=istartset(ii),iendset(ii)
463                if(ialset(j).gt.0) then
464                  vold(idir,ialset(j))=dispvelo
465                else
466                  k=ialset(j-2)
467                  do
468                    k=k-ialset(j)
469                    if(k.ge.ialset(j-1)) exit
470                    vold(idir,k)=dispvelo
471                  enddo
472                endif
473              enddo
474            endif
475          enddo
476          return
477        elseif(textpart(ij)(1:13).eq.'TYPE=VELOCITY') then
478!
479          do
480            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
481     &           ipoinp,inp,ipoinpc)
482            if((istat.lt.0).or.(key.eq.1)) return
483            read(textpart(2)(1:10),'(i10)',iostat=istat) idir
484            if(istat.gt.0) then
485              call inputerror(inpc,ipoinpc,iline,
486     &             "*INITIAL CONDITIONS%",ier)
487              return
488            endif
489            read(textpart(3)(1:20),'(f20.0)',iostat=istat) velocity
490            if(istat.gt.0) then
491              call inputerror(inpc,ipoinpc,iline,
492     &             "*INITIAL CONDITIONS%",ier)
493              return
494            endif
495            read(textpart(1)(1:10),'(i10)',iostat=istat) l
496            if(istat.eq.0) then
497              if(l.gt.nk) then
498                write(*,*)
499     &               '*WARNING reading *INITIAL CONDITIONS: node ',l
500                write(*,*)'          exceeds the largest defined ',
501     &               'node number'
502                cycle
503              endif
504              veold(idir,l)=velocity
505            else
506              read(textpart(1)(1:80),'(a80)',iostat=istat) noset
507              noset(81:81)=' '
508              ipos=index(noset,' ')
509              noset(ipos:ipos)='N'
510c              do ii=1,nset
511c                if(set(ii).eq.noset) exit
512c              enddo
513              call cident81(set,noset,nset,id)
514              ii=nset+1
515              if(id.gt.0) then
516                if(noset.eq.set(id)) then
517                  ii=id
518                endif
519              endif
520              if(ii.gt.nset) then
521                noset(ipos:ipos)=' '
522                write(*,*)
523     &               '*ERROR reading *INITIAL CONDITIONS: node set '
524     &               ,noset
525                write(*,*)'  has not yet been defined. '
526                call inputerror(inpc,ipoinpc,iline,
527     &               "*INITIAL CONDITIONS%",ier)
528                return
529              endif
530              do j=istartset(ii),iendset(ii)
531                if(ialset(j).gt.0) then
532                  veold(idir,ialset(j))=velocity
533                else
534                  k=ialset(j-2)
535                  do
536                    k=k-ialset(j)
537                    if(k.ge.ialset(j-1)) exit
538                    veold(idir,k)=velocity
539                  enddo
540                endif
541              enddo
542            endif
543          enddo
544          return
545        elseif(textpart(ij)(1:13).eq.'TYPE=PRESSURE') then
546!
547          do
548            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
549     &           ipoinp,inp,ipoinpc)
550            if((istat.lt.0).or.(key.eq.1)) return
551            read(textpart(2)(1:20),'(f20.0)',iostat=istat) pressure
552            if(istat.gt.0) then
553              call inputerror(inpc,ipoinpc,iline,
554     &             "*INITIAL CONDITIONS%",ier)
555              return
556            endif
557            read(textpart(1)(1:10),'(i10)',iostat=istat) l
558            if(istat.eq.0) then
559              if(l.gt.nk) then
560                write(*,*)
561     &               '*WARNING reading *INITIAL CONDITIONS: node ',l
562                write(*,*)'          exceeds the largest defined ',
563     &               'node number'
564                cycle
565              endif
566              vold(4,l)=pressure
567            else
568              read(textpart(1)(1:80),'(a80)',iostat=istat) noset
569              noset(81:81)=' '
570              ipos=index(noset,' ')
571              noset(ipos:ipos)='N'
572c              do ii=1,nset
573c                if(set(ii).eq.noset) exit
574c              enddo
575              call cident81(set,noset,nset,id)
576              ii=nset+1
577              if(id.gt.0) then
578                if(noset.eq.set(id)) then
579                  ii=id
580                endif
581              endif
582              if(ii.gt.nset) then
583                noset(ipos:ipos)=' '
584                write(*,*)
585     &               '*ERROR reading *INITIAL CONDITIONS: node set '
586     &               ,noset
587                write(*,*)'  has not yet been defined. '
588                call inputerror(inpc,ipoinpc,iline,
589     &               "*INITIAL CONDITIONS%",ier)
590                return
591              endif
592              do j=istartset(ii),iendset(ii)
593                if(ialset(j).gt.0) then
594                  vold(4,ialset(j))=pressure
595                else
596                  k=ialset(j-2)
597                  do
598                    k=k-ialset(j)
599                    if(k.ge.ialset(j-1)) exit
600                    vold(4,k)=pressure
601                  enddo
602                endif
603              enddo
604            endif
605          enddo
606          return
607        elseif(textpart(ij)(1:15).eq.'TYPE=TURBULENCE') then
608!
609          do
610            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
611     &           ipoinp,inp,ipoinpc)
612            if((istat.lt.0).or.(key.eq.1)) return
613            do m=5,mi(2)
614              read(textpart(m-3)(1:20),'(f20.0)',iostat=istat)
615     &             turbini(m)
616            enddo
617            if(istat.gt.0) then
618              call inputerror(inpc,ipoinpc,iline,
619     &             "*INITIAL CONDITIONS%",ier)
620              return
621            endif
622            read(textpart(1)(1:10),'(i10)',iostat=istat) l
623            if(istat.eq.0) then
624              if(l.gt.nk) then
625                write(*,*)
626     &               '*WARNING reading *INITIAL CONDITIONS: node ',l
627                write(*,*)'          exceeds the largest defined ',
628     &               'node number'
629                cycle
630              endif
631              do m=5,mi(2)
632                vold(m,l)=turbini(m)
633              enddo
634            else
635              read(textpart(1)(1:80),'(a80)',iostat=istat) noset
636              noset(81:81)=' '
637              ipos=index(noset,' ')
638              noset(ipos:ipos)='N'
639c              do ii=1,nset
640c                if(set(ii).eq.noset) exit
641c              enddo
642              call cident81(set,noset,nset,id)
643              ii=nset+1
644              if(id.gt.0) then
645                if(noset.eq.set(id)) then
646                  ii=id
647                endif
648              endif
649              if(ii.gt.nset) then
650                noset(ipos:ipos)=' '
651                write(*,*)
652     &               '*ERROR reading *INITIAL CONDITIONS: node set '
653     &               ,noset
654                write(*,*)'  has not yet been defined. '
655                call inputerror(inpc,ipoinpc,iline,
656     &               "*INITIAL CONDITIONS%",ier)
657                return
658              endif
659              do j=istartset(ii),iendset(ii)
660                if(ialset(j).gt.0) then
661                  do m=5,mi(2)
662                    vold(m,ialset(j))=turbini(m)
663                  enddo
664                else
665                  k=ialset(j-2)
666                  do
667                    k=k-ialset(j)
668                    if(k.ge.ialset(j-1)) exit
669                    do m=5,mi(2)
670                      vold(m,k)=turbini(m)
671                    enddo
672                  enddo
673                endif
674              enddo
675            endif
676          enddo
677          return
678        elseif(textpart(ij)(1:18).eq.'TYPE=TOTALPRESSURE') then
679!
680          do
681            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
682     &           ipoinp,inp,ipoinpc)
683            if((istat.lt.0).or.(key.eq.1)) return
684            read(textpart(2)(1:20),'(f20.0)',iostat=istat) totpres
685            if(istat.gt.0) then
686              call inputerror(inpc,ipoinpc,iline,
687     &             "*INITIAL CONDITIONS%",ier)
688              return
689            endif
690            read(textpart(1)(1:10),'(i10)',iostat=istat) l
691            if(istat.eq.0) then
692              if(l.gt.nk) then
693                write(*,*)
694     &               '*WARNING reading *INITIAL CONDITIONS: node ',l
695                write(*,*)'          exceeds the largest defined ',
696     &               'node number'
697                cycle
698              endif
699              vold(2,l)=totpres
700            else
701              read(textpart(1)(1:80),'(a80)',iostat=istat) noset
702              noset(81:81)=' '
703              ipos=index(noset,' ')
704              noset(ipos:ipos)='N'
705c              do ii=1,nset
706c                if(set(ii).eq.noset) exit
707c              enddo
708              call cident81(set,noset,nset,id)
709              ii=nset+1
710              if(id.gt.0) then
711                if(noset.eq.set(id)) then
712                  ii=id
713                endif
714              endif
715              if(ii.gt.nset) then
716                noset(ipos:ipos)=' '
717                write(*,*)
718     &               '*ERROR reading *INITIAL CONDITIONS: node set '
719     &               ,noset
720                write(*,*)'  has not yet been defined. '
721                call inputerror(inpc,ipoinpc,iline,
722     &               "*INITIAL CONDITIONS%",ier)
723                return
724              endif
725              do j=istartset(ii),iendset(ii)
726                if(ialset(j).gt.0) then
727                  vold(2,ialset(j))=totpres
728                else
729                  k=ialset(j-2)
730                  do
731                    k=k-ialset(j)
732                    if(k.ge.ialset(j-1)) exit
733                    vold(2,k)=totpres
734                  enddo
735                endif
736              enddo
737            endif
738          enddo
739          return
740        elseif(textpart(ij)(1:13).eq.'TYPE=MASSFLOW') then
741!
742          do
743            call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
744     &           ipoinp,inp,ipoinpc)
745            if((istat.lt.0).or.(key.eq.1)) return
746            read(textpart(2)(1:20),'(f20.0)',iostat=istat) xmassflow
747            if(iaxial.eq.180) xmassflow=xmassflow/iaxial
748            if(istat.gt.0) then
749              call inputerror(inpc,ipoinpc,iline,
750     &             "*INITIAL CONDITIONS%",ier)
751              return
752            endif
753            read(textpart(1)(1:10),'(i10)',iostat=istat) l
754            if(istat.eq.0) then
755              if(l.gt.nk) then
756                write(*,*)
757     &               '*WARNING reading *INITIAL CONDITIONS: node ',l
758                write(*,*)'          exceeds the largest defined ',
759     &               'node number'
760                cycle
761              endif
762              vold(1,l)=xmassflow
763            else
764              read(textpart(1)(1:80),'(a80)',iostat=istat) noset
765              noset(81:81)=' '
766              ipos=index(noset,' ')
767              noset(ipos:ipos)='N'
768c              do ii=1,nset
769c                if(set(ii).eq.noset) exit
770c              enddo
771              call cident81(set,noset,nset,id)
772              ii=nset+1
773              if(id.gt.0) then
774                if(noset.eq.set(id)) then
775                  ii=id
776                endif
777              endif
778              if(ii.gt.nset) then
779                noset(ipos:ipos)=' '
780                write(*,*)
781     &               '*ERROR reading *INITIAL CONDITIONS: node set '
782     &               ,noset
783                write(*,*)'  has not yet been defined. '
784                call inputerror(inpc,ipoinpc,iline,
785     &               "*INITIAL CONDITIONS%",ier)
786                return
787              endif
788              do j=istartset(ii),iendset(ii)
789                if(ialset(j).gt.0) then
790                  vold(1,ialset(j))=xmassflow
791                else
792                  k=ialset(j-2)
793                  do
794                    k=k-ialset(j)
795                    if(k.ge.ialset(j-1)) exit
796                    vold(1,k)=xmassflow
797                  enddo
798                endif
799              enddo
800            endif
801          enddo
802          return
803!
804        elseif(textpart(ij)(1:13).eq.'TYPE=SOLUTION') then
805          user=.false.
806          do j=2,n
807            if(textpart(j)(1:4).eq.'USER') user=.true.
808          enddo
809c     if(.not.user) then
810c     write(*,*)
811c     &            '*ERROR reading *INITIAL CONDITIONS: TYPE=SOLUTION'
812c     write(*,*) '       can only be used in combination with'
813c     write(*,*) '       USER'
814c     ier=1
815c     return
816c     endif
817          if(user) then
818!
819!     internal state variables are read in file sdvini.f
820!
821            iflag=1
822            ncrds=3
823            do i=1,ne
824              indexe=ipkon(i)
825              if(lakon(i)(4:4).eq.'2') then
826                nope=20
827              elseif(lakon(i)(4:4).eq.'8') then
828                nope=8
829              elseif(lakon(i)(4:5).eq.'10') then
830                nope=10
831              elseif(lakon(i)(4:4).eq.'4') then
832                nope=4
833              elseif(lakon(i)(4:5).eq.'15') then
834                nope=15
835              elseif(lakon(i)(4:4).eq.'6') then
836                nope=6
837              else
838                cycle
839              endif
840!
841              if(lakon(i)(4:5).eq.'8R') then
842                mint3d=1
843              elseif((lakon(i)(4:4).eq.'8').or.
844     &               (lakon(i)(4:6).eq.'20R')) then
845                mint3d=8
846              elseif(lakon(i)(4:4).eq.'2') then
847                mint3d=27
848              elseif(lakon(i)(4:5).eq.'10') then
849                mint3d=4
850              elseif(lakon(i)(4:4).eq.'4') then
851                mint3d=1
852              elseif(lakon(i)(4:5).eq.'15') then
853                mint3d=9
854              elseif(lakon(i)(4:4).eq.'6') then
855                mint3d=2
856              endif
857!
858              do j=1,nope
859                konl(j)=kon(indexe+j)
860                do k=1,3
861                  xl(k,j)=co(k,konl(j))
862                enddo
863              enddo
864!
865              do j=1,mint3d
866                if(lakon(i)(4:5).eq.'8R') then
867                  xi=gauss3d1(1,j)
868                  et=gauss3d1(2,j)
869                  ze=gauss3d1(3,j)
870                  weight=weight3d1(j)
871                elseif((lakon(i)(4:4).eq.'8').or.
872     &                 (lakon(i)(4:6).eq.'20R'))
873     &                 then
874                  xi=gauss3d2(1,j)
875                  et=gauss3d2(2,j)
876                  ze=gauss3d2(3,j)
877                  weight=weight3d2(j)
878                elseif(lakon(i)(4:4).eq.'2') then
879                  xi=gauss3d3(1,j)
880                  et=gauss3d3(2,j)
881                  ze=gauss3d3(3,j)
882                  weight=weight3d3(j)
883                elseif(lakon(i)(4:5).eq.'10') then
884                  xi=gauss3d5(1,j)
885                  et=gauss3d5(2,j)
886                  ze=gauss3d5(3,j)
887                  weight=weight3d5(j)
888                elseif(lakon(i)(4:4).eq.'4') then
889                  xi=gauss3d4(1,j)
890                  et=gauss3d4(2,j)
891                  ze=gauss3d4(3,j)
892                  weight=weight3d4(j)
893                elseif(lakon(i)(4:5).eq.'15') then
894                  xi=gauss3d8(1,j)
895                  et=gauss3d8(2,j)
896                  ze=gauss3d8(3,j)
897                  weight=weight3d8(j)
898                elseif(lakon(i)(4:4).eq.'6') then
899                  xi=gauss3d7(1,j)
900                  et=gauss3d7(2,j)
901                  ze=gauss3d7(3,j)
902                  weight=weight3d7(j)
903                endif
904!
905                if(nope.eq.20) then
906                  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
907                elseif(nope.eq.8) then
908                  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
909                elseif(nope.eq.10) then
910                  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
911                elseif(nope.eq.4) then
912                  call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
913                elseif(nope.eq.15) then
914                  call shape15w(xi,et,ze,xl,xsj,shp,iflag)
915                else
916                  call shape6w(xi,et,ze,xl,xsj,shp,iflag)
917                endif
918!
919                do k=1,3
920                  pgauss(k)=0.d0
921                  do i1=1,nope
922                    pgauss(k)=pgauss(k)+
923     &                   shp(4,i1)*co(k,konl(i1))
924                  enddo
925                enddo
926!
927                call sdvini(xstate(1,j,i),pgauss,nstate_,ncrds,
928     &               i,j,layer,kspt)
929!
930              enddo
931            enddo
932            call getnewline(inpc,textpart,istat,n,key,iline,ipol,
933     &           inl,ipoinp,inp,ipoinpc)
934            return
935          else
936!
937!     internal variables are written explicitly in the input deck
938!
939            do
940              call getnewline(inpc,textpart,istat,n,key,iline,ipol,
941     &             inl,ipoinp,inp,ipoinpc)
942              if((istat.lt.0).or.(key.eq.1)) return
943!
944              if(nstate_.lt.6) then
945                ntot=nstate_
946              else
947                ntot=6
948              endif
949!
950              do j=1,ntot
951                read(textpart(j+2)(1:20),'(f20.0)',iostat=istat)
952     &               beta(j)
953                if(istat.gt.0) then
954                  call inputerror(inpc,ipoinpc,iline,
955     &                 "*INITIAL CONDITIONS%",ier)
956                  return
957                endif
958              enddo
959              read(textpart(1)(1:10),'(i10)',iostat=istat) l
960              if(istat.ne.0) then
961                call inputerror(inpc,ipoinpc,iline,
962     &               "*INITIAL CONDITIONS%",ier)
963                return
964              endif
965              if(l.gt.ne) then
966                write(*,*)
967     &               '*WARNING reading *INITIAL CONDITIONS: element ',l
968                write(*,*)'          exceeds the largest defined ',
969     &               'element number'
970                cycle
971              endif
972              read(textpart(2)(1:10),'(i10)',iostat=istat) k
973              if(istat.eq.0) then
974                do j=1,ntot
975                  xstate(j,k,l)=beta(j)
976                enddo
977              else
978                call inputerror(inpc,ipoinpc,iline,
979     &               "*INITIAL CONDITIONS%",ier)
980                return
981              endif
982!
983              if(nstate_.gt.6) then
984                numberoflines=(nstate_-7)/8+1
985                do ii=1,numberoflines
986                  if(ii.lt.numberoflines) then
987                    jmax=8
988                  else
989                    jmax=nstate_-ntot
990                  endif
991                  call getnewline(inpc,textpart,istat,n,key,iline,
992     &                 ipol,inl,ipoinp,inp,ipoinpc)
993                  if((istat.lt.0).or.(key.eq.1)) return
994                  do j=1,jmax
995                    read(textpart(j+2)(1:20),'(f20.0)',
996     &                   iostat=istat) beta(j)
997                    if(istat.gt.0)
998     &                   call inputerror(inpc,ipoinpc,iline,
999     &                   "*INITIAL CONDITIONS%",ier)
1000                    return
1001                    xstate(ntot+j,k,l)=beta(j)
1002                  enddo
1003                  ntot=ntot+jmax
1004                enddo
1005              endif
1006!
1007            enddo
1008            return
1009          endif
1010!
1011        else
1012          write(*,*)
1013     &'*WARNING reading *INITIAL CONDITIONS: parameter not recognized:'
1014          write(*,*) '         ',
1015     &         textpart(ij)(1:index(textpart(ij),' ')-1)
1016          call inputwarning(inpc,ipoinpc,iline,
1017     &         "*INITIAL CONDITIONS%")
1018        endif
1019      enddo
1020!
1021      write(*,*) '*ERROR reading *INITIAL CONDITIONS: unknown type'
1022      write(*,*) '  '
1023      call inputerror(inpc,ipoinpc,iline,
1024     &     "*INITIAL CONDITIONS%",ier)
1025!
1026      return
1027      end
1028
1029