36
37
38
39 USE python_funct_mod
40 USE finter_mixed_mod
41 USE sensor_mod
42
43
44
45#include "implicit_f.inc"
46
47
48
49#include "com04_c.inc"
50#include "com06_c.inc"
51#include "com08_c.inc"
52#include "units_c.inc"
53#include "scr02_c.inc"
54#include "scr07_c.inc"
55#include "scr18_c.inc"
56#include "task_c.inc"
57#include "param_c.inc"
58
59
60
61 INTEGER ,INTENT(IN) :: NSENSOR
62 INTEGER NPC(*),IVOLU(*),NJET,IBAGJET(NIBJET,*),
63 . NVENT,IBAGHOL(NIBHOL,*),ICBAG(NICBAG,*),PMAIN
65 . tf(*),rvolu(*),rbagjet(nrbjet,*),
66 . rbaghol(nrbhol,*),rcbag(nrcbag,*),rvoluv(nrvolu,*),vol
67 TYPE (SENSOR_STR_) ,DIMENSION(NSENSOR) ,INTENT(IN) :: SENSOR_TAB
68 DOUBLE PRECISION, INTENT(INOUT) :: WFEXT
69 TYPE(PYTHON_), intent(inout) :: PYTHON
70
71
72
73 INTEGER I, IMASS, ITEMP, ISENS, IFLU, IEQUI,
74 . IDEF, NAV, II, IS, IDTPDEF,IDPDEF,
75 . IV, IJ,ITTF
77 . pdef, pext, pvois, dtpdefc,
78 . gama, amu, tstart, tsg,
79 . amtot, p, pold, dv,
80 . rot, ftemp, fmass, dt,
area,
81 . cv, cp, cvg, cpg, cpa, cpb, cpc, cvi, cpi, cpai, cpbi, cpci,
82 . rmwi, rmwg,
83 . rnm_old, rnm, rnmg,
84 . gmtot, gmtot_old, gmass, gmass_old, dgmass, dgmout,
85 . gmi, gmi_old,
86 . efac, dgeout, right, left,
87 . vold, veps, vmin, amtot_old, tbag_old, tbag,scalt,
88 . qold,q,cx,qx,temp
89
90 pext =rvolu(3)
91 pold =rvolu(12)
92 scalt =rvolu(26)
93 efac = zero
94
95 DO iv=1,nvent
96 idef = ibaghol(1,iv)
97 pdef = rbaghol(1,iv)
98 dtpdefc= rbaghol(5,iv)
99 idtpdef= ibaghol(11,iv)
100 idpdef = ibaghol(12,iv)
101
102 IF (idtpdef == 0) THEN
103 IF(idef == 0 .AND. pold > pdef+pext) THEN
104 dtpdefc = dtpdefc+dt1
105 END IF
106 ELSE IF (idtpdef == 1) THEN
107 IF (pold > pdef+pext) THEN
108 idpdef = 1
109 END IF
110 IF (idpdef == 1) THEN
111 dtpdefc = dtpdefc+dt1
112 END IF
113 END IF
114 ibaghol(12,iv) = idpdef
115 rbaghol(5,iv) = dtpdefc
116 ENDDO
117 IF(ispmd+1 == pmain) THEN
118 DO iv=1,nvent
119 rbaghol(19,iv)= rbaghol(19,iv)+rbaghol(21,iv)*dt1
120 rbaghol(20,iv)= rbaghol(20,iv)+rbaghol(22,iv)*dt1
121 ENDDO
122 END IF
123
124 nav = ivolu(3)
125 DO i=1,nav
126 ii = icbag(1,i)
127 idef = icbag(3,i)
128 pdef = rcbag(1,i)
129 dtpdefc= rcbag(5,i)
130 pvois = rvoluv(31,ii)
131 IF(idef == 0 .and .pold > pdef+pvois)rcbag(5,i)=dtpdefc+dt1
132 ENDDO
133
134 amu =rvolu(2)
135
136 iequi=ivolu(15)
137
138 is = 0
139 ittf=ivolu(17)
140 DO ij=1,njet
141 isens=ibagjet(4,ij)
142 IF(isens == 0)THEN
143 tstart=zero
144 ELSE
145 tstart=sensor_tab(isens)%TSTART
146 ENDIF
147 IF(tt >= tstart)is=1
148 IF(iequi == 0)THEN
149 dgmout=rbagjet(9,ij)
150 IF(dgmout /= zero)is=1
151 END IF
152 IF (is == 1 .AND. (ittf == 1 .OR. ittf == 2 .OR. ittf == 3)) THEN
153 ittf=ittf+10
154
155 rvolu(60)=tstart
156 ivolu(17)=ittf
157 END IF
158 ENDDO
159 IF(iequi == 0)THEN
160 tbag_old =rvolu(13)
161 vold =rvolu(16)
162 veps =rvolu(17)
164 rot =rvolu(21)
165 qold =rvolu(23)
166 vol = vol + veps
167
168 dgmout=rvolu(24)
169 IF(dgmout/=zero)is=1
170
171 dv = vol-vold
172 IF(dv /= zero)is=1
173 ELSE
174 IF(is == 1)THEN
175 tbag_old =rvolu(13)
177 IF(iequi==1)THEN
178
179 vmin =em4*
area**three_half
180 veps =
max(zero,vmin-vol)
181 rvolu(17)=veps
182
183
184 rmwi = rvolu(10)
185 gmi = pext*(vol+veps)/(rmwi*tbag_old)
186 rvolu(11)= gmi
187 rvolu(14)= rmwi*gmi
188 rvolu(20)= gmi
189 IF(ispmd+1==pmain) THEN
190 WRITE(iout,*)' *** MONITORED VOLUME : INITIAL EQUILIBRIUM IS SET ***'
191 WRITE(iout,'(A,I10,A,G20.13,A)') ' *** MONITORED VOLUME ',ivolu(1),' VOLUME ',vol,' ***'
192 WRITE(istdo,*)' *** MONITORED VOLUME : INITIAL EQUILIBRIUM IS SET ***'
193 ENDIF
194 ivolu(15)= -1
195 vol =vol + veps
196 vold =vol
197 qold =zero
198 ELSE
199 vold =rvolu(16)
200 veps =rvolu(17)
201 qold =rvolu(23)
202 vol = vol + veps
203 END IF
204 rot =rvolu(21)
205 dv =vol-vold
206 END IF
207 END IF
208
209
210 IF (is /= 0)THEN
211
212
213
214
215
216 rnm_old =rvolu(14)
217 amtot_old=rvolu(20)
218 left =zero
219 right=zero
220 amtot=amtot_old
221
222
223
224 DO ij=1,njet
225 fmass=rbagjet(5,ij)
226 ftemp=rbagjet(6,ij)
227 imass=ibagjet(1,ij)
228 iflu =ibagjet(2,ij)
229 itemp=ibagjet(3,ij)
230 isens=ibagjet(4,ij)
231 IF(isens==0)THEN
232 tstart=zero
233 ELSE
234 tstart=sensor_tab(isens)%TSTART
235 ENDIF
236 rmwg=rbagjet(1,ij)
237 cpa =rbagjet(2,ij)
238 cpb =rbagjet(3,ij)
239 cpc =rbagjet(4,ij)
240 IF(tt >= tstart)THEN
241 tsg = (tt-tstart)*rvolu(26)
242 IF(itemp > 0) THEN
243 temp=ftemp*finter_mixed(python,nfunct,itemp,tsg,npc,tf)
244 ELSE
245 temp=ftemp
246 ENDIF
247 efac= temp*(cpa+half*cpb*temp+third*cpc*temp*temp)
248 ELSE
249 ENDIF
250
251 gmtot_old=rbagjet(8,ij)
252 gmass_old=rbagjet(7,ij)
253 IF(tt >= tstart)THEN
254 IF(imass > 0) THEN
255 gmass=fmass*finter_mixed(python,nfunct,imass,tsg,npc,tf)
256 IF(iflu == 1)gmass = gmass*rvolu(26)*dt1 + gmass_old
257 ELSE
258 gmass=fmass
259 ENDIF
260 dgmass=
max(zero,gmass-gmass_old)
261 right =right+dgmass*efac
262 ELSE
263 dgmass=zero
264 gmass =zero
265 ENDIF
266 dgmout=rbagjet( 9,ij)*dt1
267 gmtot =gmtot_old+dgmass-dgmout
268 dgeout=rbagjet(10,ij)*dt1
269 right= right-dgeout
270 right= right+(dgmout-dgmass)*tbag_old*(cpa+half*cpb*tbag_old+third*cpc*tbag_old*tbag_old-rmwg)
271 cvg = cpa+cpb*tbag_old+cpc*tbag_old*tbag_old-rmwg
272 left = left +gmtot*cvg
273 right= right+gmtot*cvg*tbag_old
274 rbagjet(8,ij)=gmtot
275 rbagjet(7,ij)=gmass
276 amtot=amtot+gmtot-gmtot_old
277 ENDDO
278
279
280
281
282 dgmout=rvolu(24)*dt1
283 cpai =rvolu(7)
284 cpbi =rvolu(8)
285 cpci =rvolu(9)
286 rmwi =rvolu(10)
287 gmi_old=rvolu(11)
288 cvi = cpai+cpbi*tbag_old+cpci*tbag_old*tbag_old-rmwi
289 gmi=gmi_old-dgmout
290 dgeout=rvolu(22)*dt1
291
292 right=right+dgmout*tbag_old*(cpai+half*cpbi*tbag_old+third*cpci*tbag_old*tbag_old-rmwi)-dgeout
293 left = left+gmi*cvi
294 right = right+gmi*cvi*tbag_old
295 amtot=amtot-dgmout
296 rvolu(11)=gmi
297
298
299
300 gmi =rvolu(11)
301 cpai=rvolu(7)
302 cpbi=rvolu(8)
303 cpci=rvolu(9)
304 rmwi=rvolu(10)
305 rnm =gmi*rmwi
306 DO ij=1,njet
307 rmwg =rbagjet(1,ij)
308 gmtot=rbagjet(8,ij)
309 rnmg =gmtot*rmwg
310 rnm =rnm+rnmg
311 ENDDO
312
313 left = left + half*rnm*dv/vol
314 right = right- half*rnm_old*tbag_old*dv/vold
315
316 tbag = right/left
317 tbag =
max(tbag,zero)
318 p=rnm*tbag/vol
319
320
321
322 cpi=cpai+cpbi*tbag+cpci*tbag*tbag
323 cvi=cpi-rmwi
324 cp =gmi*cpi
325 cv =gmi*cvi
326 DO ij=1,njet
327 rmwg =rbagjet(1,ij)
328 cpa =rbagjet(2,ij)
329 cpb =rbagjet(3,ij)
330 cpc =rbagjet(4,ij)
331 gmtot=rbagjet(8,ij)
332 cpg =cpa+cpb*tbag+cpc*tbag*tbag
333 cvg =cpg-rmwg
334 cp =cp+gmtot*cpg
335 cv =cv+gmtot*cvg
336 ENDDO
337 gama=cp/cv
338
339 rvolu(1) =gama
340 rvolu(14)=rnm
341
342 IF(dt1==zero.OR.dv>zero)THEN
343 q=zero
344 ELSE
345 q=-amu*sqrt(p*
area*rot/vol)*dv/
area/dt1
346 ENDIF
347 IF (ispmd+1==pmain) THEN
348 wfext=wfext+(half*(p+pold+q+qold)-pext)*dv
349 ENDIF
350
351 rvolu(12)=p
352 rvolu(13)=tbag
353 rvolu(16)=vol
354 rvolu(20)=amtot
355 rvolu(23)=q
356
357 rvolu(22)=zero
358 rvolu(24)=zero
359 DO ij=1,njet
360 rbagjet( 9,ij)=zero
361 rbagjet(10,ij)=zero
362 ENDDO
363
364
365 cx = sqrt(two*gama*p*vol/(gama-one)/(amtot+
area*rot))
366 qx = amu*cx
367 dt = dtfac1(9)*vol/
area/
max(em20,qx+sqrt(qx*qx+cx*cx))
368 IF(dt<dt2)THEN
369 dt2=dt
370 nelts =ivolu(1)
371 itypts=9
372 ENDIF
373 IF(idtmin(9) == 1 .AND. dt < dtmin1(9)) THEN
374 tstop = tt
375 IF (ispmd+1==pmain) THEN
376 WRITE(iout,*) '-- MINIMUM MONITORED VOLUME TIME STEP '
377 WRITE(istdo,*) '-- MINIMUM MONITORED VOLUME TIME STEP '
378 ENDIF
379 ELSEIF(idtmin(9) == 5 .AND. dt < dtmin1(9)) THEN
380 mstop = 2
381 IF (ispmd+1==pmain) THEN
382 WRITE(iout,*) '-- MINIMUM MONITORED VOLUME TIME STEP '
383 WRITE(istdo,*) '-- MINIMUM MONITORED VOLUME TIME STEP '
384 ENDIF
385 ENDIF
386 ENDIF
387
388 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)