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 con2phys(vold,vcon,nk,ntmat_,shcon,nshcon,rhcon,nrhcon,
20     &     physcon,ithermal,compressible,iturbulent,inomat,mi,
21     &     ierr,ifreesurface,dgravity,depth)
22!
23!     calculates the physical variable from the conservative
24!     variables:
25!     for gases: temperature,velocity and pressure
26!     for thermal liquids: temperature, velocity and density
27!     for athermal liquids: velocity and density
28!
29      implicit none
30!
31      integer compressible,ifreesurface,nrhcon(*),ntmat_,iturbulent,
32     &     mi(*),nshcon(*),nk,ithermal(*),i,j,k,imat,ierr,inomat(*)
33!
34      real*8 vold(0:mi(2),*),vcon(nk,0:mi(2)),rhcon(0:1,ntmat_,*),rho,
35     &     c1,cp,r,temp,temp0,c2,shcon(0:3,ntmat_,*),physcon(*),
36     &     dgravity,depth(*)
37!
38      if(ithermal(1).gt.1) then
39!
40        do i=1,nk
41          imat=inomat(i)
42          temp=vold(0,i)
43!
44          if(compressible.eq.1) then
45!
46!     gas : rho*epsilon_t, rho and rho*v known
47!
48!     calculation of the static temperature
49!
50            rho=vcon(i,4)
51            r=shcon(3,1,imat)
52            c1=(vcon(i,0)-(vcon(i,1)**2+vcon(i,2)**2+
53     &           vcon(i,3)**2)/(2.d0*rho))/rho
54!
55!     the energy density per volume unit cannot be
56!     negative
57!
58            if(c1.lt.1.d-10) c1=1.d-10
59!
60            temp0=temp
61            j=0
62            do
63              call materialdata_cp_sec(imat,ntmat_,temp,shcon,
64     &             nshcon,cp,physcon)
65              temp=c1/(cp-r)+physcon(1)
66              j=j+1
67              if((dabs(temp-temp0).lt.1.d-4*dabs(temp)).or.
68     &             (dabs(temp-temp0).lt.1.d-10)) then
69                vold(0,i)=temp
70                exit
71              endif
72              if(j.gt.100) then
73                write(*,*)
74     &               '*ERROR in con2phys: too many iterations'
75                write(*,*) '       for node',i
76                write(*,*) '       increment is recalculated'
77                write(*,*) '       with a higher shock smoothing'
78                ierr=1
79                return
80              endif
81              temp0=temp
82            enddo
83!
84            if(ifreesurface.eq.0) then
85!
86!     determining the pressure (gas equation)
87!
88              vold(4,i)=rho*r*(temp-physcon(1))
89            else
90!
91!     determining the pressure (shallow water)
92!
93              vold(4,i)=dgravity*(vcon(i,4)**2-depth(i)**2)/2.d0
94            endif
95!
96!     calculating the velocity from the conservative variables
97!
98            do k=1,3
99              vold(k,i)=vcon(i,k)/rho
100            enddo
101!
102!     calculating the turbulent quantities from the conservative
103!     variables
104!
105            if(iturbulent.ne.0) then
106              do k=5,6
107                vold(k,i)=vcon(i,k)/rho
108              enddo
109            endif
110!
111            cycle
112          else
113!
114!     thermal liquid
115!
116            c1=vcon(i,0)
117            c2=(vcon(i,1)**2+vcon(i,2)**2+vcon(i,3)**2)/2.d0
118            temp0=temp
119            j=0
120!
121!     iterating to find the temperature
122!
123            do
124              call materialdata_cp_sec(imat,ntmat_,temp,shcon,
125     &             nshcon,cp,physcon)
126              call materialdata_rho(rhcon,nrhcon,imat,rho,
127     &             temp,ntmat_,ithermal)
128              temp=(c1-c2/rho)/(rho*cp)+physcon(1)
129              j=j+1
130              if((dabs(temp-temp0).lt.1.d-4*dabs(temp)).or.
131     &             (dabs(temp-temp0).lt.1.d-10)) then
132                vold(0,i)=temp
133                exit
134              endif
135              if(j.gt.100) then
136                write(*,*)
137     &               '*ERROR in con2phys: too many iterations'
138                write(*,*) '       for node',i
139                write(*,*) '       actual temperature ',temp,' K'
140                stop
141              endif
142              temp0=temp
143            enddo
144!
145!     storing the density
146!
147            vcon(i,4)=rho
148!
149!     calculating the velocity
150!
151            do k=1,3
152              vold(k,i)=vcon(i,k)/rho
153            enddo
154!
155!     calculating the turbulent quantities
156!
157            if(iturbulent.ne.0) then
158              do k=5,6
159                vold(k,i)=vcon(i,k)/rho
160              enddo
161            endif
162!
163          endif
164        enddo
165      else
166!
167!     athermal calculations
168!
169        do i=1,nk
170!
171!     calculating the fictitious pressure for shallow water
172!     calculations
173!
174          if(ifreesurface.eq.1) then
175            vold(4,i)=dgravity*(vcon(i,4)**2-depth(i)**2)/2.d0
176          endif
177!
178!     calculating the velocity
179!
180          rho=vcon(i,4)
181          do k=1,3
182            vold(k,i)=vcon(i,k)/rho
183          enddo
184!
185!     calculating the turbulent quantities
186!
187          if(iturbulent.ne.0) then
188            do k=5,6
189              vold(k,i)=vcon(i,k)/rho
190            enddo
191          endif
192!
193        enddo
194      endif
195!
196      return
197      end
198