32
33
34
35#include "implicit_f.inc"
36
37
38
39#include "com04_c.inc"
40#include "com08_c.inc"
41#include "param_c.inc"
42
43
44
45 INTEGER MAT(*),NGEO(*)
46 INTEGER, INTENT(IN) :: LFT, LLT, JPOR, JCLOSE
47
48 my_real pm(npropm,nummat), off(*), rho(*), rk(*), re(*),geo(npropg,numgeo)
49 my_real einc(*), dvol(*), vnew(*), tmu(*), vis(*), vd2(*)
50 my_real,
intent(inout) :: pturb(1:llt)
51
52
53
54 INTEGER I
55 my_real ei, xfac, ark, c1, c2, c3, sgsl, cmu, se, fac, resgs
56
57
58
59
60 IF(jpor /= 2)THEN
61
62 DO i=lft,llt
63 xfac=tmu(i)/vis(i)
64 pturb(i)=two*rk(i)/(vnew(i)*three)
65 ei =xfac*einc(i)-dvol(i)*pturb(i)
67 rk(i) =rk(i)+(ei-re(i)*dt1)*off(i)
69 c1 =pm(82,mat(i))/ark
70 c2 =pm(83,mat(i))/ark
71 c3 =pm(84,mat(i))
72 re(i)=re(i)*(one +(c1*ei-c2*re(i)*dt1+c3*dvol(i)/vnew(i))*off(i))
73 rk(i)=
max(em15,rk(i)/vnew(i))
74 re(i)=
max(em15,re(i)/vnew(i))
75 ENDDO
76 IF(jclose == 0)THEN
77 DO i=lft,llt
78 sgsl=pm(46,mat(i))
79 cmu=pm(81,mat(i))
80 c1 =pm(82,mat(i))
81 c2 =pm(83,mat(i))
82 se =pm(86,mat(i))
83 fac=sqrt(cmu/(se*(c2-c1)))/sgsl
84 resgs=rho(i)*fac*(rk(i)/rho(i))**three_half
85 re(i)=
max(re(i),resgs)
86 ENDDO
87 ELSE
88 DO i=lft,llt
89 sgsl=
max(2.*geo(130,ngeo(i)),pm(46,mat(i)))
90 cmu=pm(81,mat(i))
91 c1 =pm(82,mat(i))
92 c2 =pm(83,mat(i))
93 se =pm(86,mat(i))
94 fac=sqrt(cmu/(se*(c2-c1)))/
max(sgsl,em15)
95 resgs=rho(i)*fac*(rk(i)/rho(i))**three_half
96 re(i)=
max(re(i),resgs)
97 ENDDO
98 ENDIF
99 ELSE
100
101 DO i=lft,llt
102
103 rk(i)=rho(i)*geo(22,ngeo(i))**2*vd2(i)
104 re(i)=(rk(i)**3/rho(i))**half/geo(23,ngeo(i))
105 ENDDO
106
107 endif
108
109 RETURN