51
52
53
54 use glob_therm_mod
55 USE prop_param_mod , only : n_var_igeo
56
57
58
59#include "implicit_f.inc"
60
61
62
63#include "mvsiz_p.inc"
64
65
66
67#include "com08_c.inc"
68#include "param_c.inc"
69#include "scr05_c.inc"
70
71
72
73 INTEGER, INTENT(IN) :: JSMS
74 INTEGER, INTENT(IN) :: JTUR
75 INTEGER, INTENT(IN) :: ITY
76 INTEGER, INTENT(IN) :: JTHE
77 INTEGER, INTENT(IN) :: ISMSTR
78 INTEGER, INTENT(IN) :: JSPH,NPG
79 INTEGER,INTENT(IN) :: NUMGEO
80
81 INTEGER NEL
83 . pm(npropm,*), sig(nel,6), rho(mvsiz),
84 . d1(*), d2(*), d3(*), d4(*),d5(*), d6(*),
85 . bpreld(nel,*), off(*)
86 INTEGER NELTST,ITYPTST,PID(*),G_DT
87 INTEGER MAT(*),NGL(*),IPM(NPROPMI,*)
89 . dt2t
90
92 . eint(*), qold(*),
93 . vol0(*),stifn(*), offg(*),geo(npropg,*),mumax(*)
95 . vnew(*), vd2(*), deltax(*), ssp(*), aire(*), vis(*),
96 . psh(*), pnew(*),qnew(*) ,ssp_eq(*), dvol(*),
97 . sold1(*), sold2(*), sold3(*), sold4(*), sold5(*), sold6(*),
98 . mssa(*), dmels(*),conde(*),amu(*),vol_avg(*),
dtel(*),
99 . rhoref(*) ,rhosp(*)
100 DOUBLE PRECISION VOL0DP(*)
101 type (glob_therm_) ,intent(inout) :: glob_therm
102 integer,dimension(n_var_igeo,numgeo),intent(in) :: igeo
103
104
105
106 INTEGER I, MX
108 . reduc ,reduc1 , t1, t2 , ee, nu, ggt, treps, p, tp
110 . c1(mvsiz), lmbd(mvsiz), gg(mvsiz),
111 . rho0(mvsiz), rho0_1,c1_1
112
114 . facq0,
115 . e1, e2, e3, e4, e5, e6, einc,
116 . bid1, bid2, bid3, dta
117
118 facq0 = one
119 reduc1 = em04
120
121 DO i=1,nel
122 t1 = bpreld(i,1)+zep4*(bpreld(i,2)-bpreld(i,1))
123 t2 = bpreld(i,1)+zep7*(bpreld(i,2)-bpreld(i,1))
124 reduc = reduc1
125 tp = tt-t1
126 IF(tp > zero ) THEN
127 reduc =
min(reduc1*(one-tp/(t2-t1))+tp/(t2-t1), one)
128 ENDIF
129
130 mx = mat(i)
131 rho0_1 = pm(1, mx)
132 nu = pm(21, mx)
133 ee = pm(20, mx)*reduc
134 gg(i) = pm(22, mx)*reduc
135 c1_1 = pm(32, mx)
136 lmbd(i) = dt1*ee*nu/((one+nu)*(one-two*nu))
137 c1(i) = c1_1*reduc
138 rho0(i) =rho0_1
139
140
141 ssp(i)=sqrt((onep333*pm(22, mx)+c1_1)/rho0(i))
142 ENDDO
143
144 IF (jsph==0)THEN
146 1 pm, off, rho, bid1,
147 2 bid2, ssp, bid3, stifn,
148 3 dt2t, neltst, ityptst, aire,
149 4 offg, geo, pid, vnew,
150 5 vd2, deltax, vis, d1,
151 6 d2, d3, pnew, psh,
152 7 mat, ngl, qnew, ssp_eq,
153 8 vol0, mssa, dmels, igeo,
154 9 facq0, conde,
dtel, g_dt,
155 a ipm, rhoref, rhosp, nel,
156 b ity, ismstr, jtur, jthe,
157 c jsms, npg , glob_therm)
158 ELSE
160 1 pm, off, rho, bid1,
161 2 bid2, bid3, stifn, dt2t,
162 3 neltst, ityptst, offg, geo,
163 4 pid, mumax, ssp, vnew,
164 5 vd2, deltax, vis, d1,
165 6 d2, d3, pnew, psh,
166 7 mat, ngl, qnew, ssp_eq,
167 8 g_dt,
dtel, nel, ity,
168 9 jtur, jthe)
169 ENDIF
170
171 dta =half*dt1
172
173 DO i=1,nel
174 pnew(i) = c1(i)*amu(i)
175 treps = d1(i)+d2(i)+d3(i)
176 ggt = dt1*gg(i)
177 sig(i,1) = sig(i,1)+two*ggt*d1(i)+lmbd(i)*treps
178 sig(i,2) = sig(i,2)+two*ggt*d2(i)+lmbd(i)*treps
179 sig(i,3) = sig(i,3)+two*ggt*d3(i)+lmbd(i)*treps
180 sig(i,4) = sig(i,4)+ggt*d4(i)
181 sig(i,5) = sig(i,5)+ggt*d5(i)
182 sig(i,6) = sig(i,6)+ggt*d6(i)
183 qold(i) = qnew(i)
184 ENDDO
185
186
187
188
189
190
191
192 IF(ismstr == 1)THEN
193 IF(iresp==0)THEN
194 DO i=1,nel
195 IF(tt >= (bpreld(i,2)-two*dt1)) THEN
196 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
197 rho(i)=rho0_1*(one+p/c1_1)
198
199
200 ENDIF
201 ENDDO
202 ELSE
203 DO i=1,nel
204 IF(tt >= (bpreld(i,2)-two*dt1)) THEN
205 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
206 rho(i)=rho0_1*(one+p/c1_1)
207 vol0dp(i) = vol0(i)*(one+p/c1_1)
208
209 ENDIF
210 ENDDO
211 END IF
212 ELSEIF(ismstr==2)THEN
213 IF(iresp==0)THEN
214 DO i=1,nel
215 IF(offg(i) > one)THEN
216 IF(tt >= (bpreld(i,2)-two*dt1)) THEN
217 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
218 rho(i)=rho0_1*(one+p/c1_1)
219
220
221 ENDIF
222 ELSE
223 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
224 vol0(i) = vnew(i)*(one+p/c1_1)
225 rho(i) = rho0_1*vol0(i)/vnew(i)
226
227
228 END IF
229 ENDDO
230 ELSE
231 DO i=1,nel
232 IF(offg(i) > one)THEN
233 IF(tt >= (bpreld(i,2)-two*dt1)) THEN
234 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
235 rho(i) = rho0_1*(one+p/c1_1)
236 vol0dp(i) = vol0(i)*(one+p/c1_1)
237
238 ENDIF
239 ELSE
240 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
241 vol0(i) = vnew(i)*(one+p/c1_1)
242 vol0dp(i) = vol0(i)
243
244 rho(i) = rho0_1*vol0(i)/vnew(i)
245 END IF
246 ENDDO
247 END IF
248 ELSEIF(ismstr<=4)THEN
249 IF(iresp==0)THEN
250 DO i=1,nel
251 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
252 vol0(i)=vnew(i)*(one+p/c1_1)
253 rho(i)= rho0_1*vol0(i)/vnew(i)
254
255 ENDDO
256 ELSE
257 DO i=1,nel
258 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
259 vol0(i) = vnew(i)*(one+p/c1_1)
260 rho(i) = rho0_1*vol0(i)/vnew(i)
261 vol0dp(i) = vol0(i)
262
263 ENDDO
264 END IF
265 END IF
266
267 DO i=1,nel
268 e1=d1(i)*(sold1(i)+sig(i,1))
269 e2=d2(i)*(sold2(i)+sig(i,2))
270 e3=d3(i)*(sold3(i)+sig(i,3))
271 e4=d4(i)*(sold4(i)+sig(i,4))
272 e5=d5(i)*(sold5(i)+sig(i,5))
273 e6=d6(i)*(sold6(i)+sig(i,6))
274 einc= zero
275 eint(i)=(eint(i)+einc*off(i)) /
max(em15,vol0(i))
276 ENDDO
277
278 RETURN
subroutine dtel(ssp, pm, geo, pid, mat, rho0, vis, deltax, aire, vol, dtx)
subroutine mdtsph(pm, off, rho, rk, t, re, sti, dt2t, neltst, ityptst, offg, geo, pid, mumax, ssp, vol, vd2, deltax, vis, d1, d2, d3, pnew, psh, mat, ngl, qvis, ssp_eq, g_dt, dtsph, nel, ity, jtur, jthe)
subroutine mqviscb(pm, off, rho, rk, temp, ssp, re, sti, dt2t, neltst, ityptst, aire, offg, geo, pid, vol, vd2, deltax, vis, d1, d2, d3, pnew, psh, mat, ngl, qvis, ssp_eq, vol0, mssa, dmels, igeo, facq0, conde, dtel, g_dt, ipm, rhoref, rhosp, nel, ity, ismstr, jtur, jthe, jsms, npg, glob_therm)