36
37
38
39#include "implicit_f.inc"
40
41
42
43#include "param_c.inc"
44
45
46
47 INTEGER NEL
48 INTEGER NGL(NEL)
49 my_real,
DIMENSION(NEL,6) :: sig
50 my_real,
DIMENSION(NPROPM) :: pm
51 my_real,
DIMENSION(NEL),
INTENT(IN) :: s01,s02,s03,s04,s05,s06,
52 . scal1,scal2,scal3,vk0,rob,off,seq
53 my_real,
DIMENSION(NEL),
INTENT(OUT) :: vk,s1,s2,s3,s4,s5,s6,
54 . scle1,scle2,scle3,sm,dsm
55
56
57
58 INTEGER I,N,NIT,IBUG,NINDEX,NINDEX2,ICRIT_OUTP
59 INTEGER INDEX()
60 my_real h1, h2, h3, h4, h5, h6, ro0,rok0, tolf,
61 . fc,rt,rc,rct1,rct2,aa,ac,bc,bt,tol
62 my_real ds1(nel),ds2(nel),ds3(nel),ds4(nel),ds5(nel),ds6(nel),
63 . fa(nel), xn(nel),fn(nel),sn1(nel),sn2(nel),
64 . sn3(nel),sn4(nel),sn5(nel),sn6(nel),sn7(nel),rok(nel)
65
66 DATA tolf/0.005/
67
68
69
70
71 DO i=1,nel
72 s1(i) = s01(i) * scal1(i)
73 s2(i) = s02(i) * scal2(i)
74 s3(i) = s03(i) * scal3(i)
75 s4(i) = s04(i) * scal1(i)*scal2(i)
76 s5(i) = s05(i) * scal2(i)*scal3(i)
77 s6(i) = s06(i) * scal3(i)*scal1(i)
78 sm(i) = third * (s1(i) + s2(i) + s3(i))
79
80 ds1(i) = (s1(i) - sig(i,1)) * scal1(i)
81 ds2(i) = (s2(i) - sig(i,2)) * scal2(i)
82 ds3(i) = (s3(i) - sig(i,3)) * scal3(i)
83 ds4(i) = (s4(i) - sig(i,4)) * scal1(i)*scal2(i)
84 ds5(i) = (s5(i) - sig(i,5)) * scal2(i)*scal3(i)
85 ds6(i) = (s6(i) - sig(i,6)) * scal3(i)*scal1(i)
86 dsm(i) = third * (ds1(i)+ds2(i)+ds3(i))
87 ENDDO
88
89 DO i=1,nel
90 s1(i) = s1(i) -sm(i)
91 s2(i) = s2(i) -sm(i)
92 s3(i) = s3(i) -sm(i)
93 ds1(i) = ds1(i)-dsm(i)
94 ds2(i) = ds2(i)-dsm(i)
95 ds3(i) = ds3(i)-dsm(i)
96 ENDDO
97
98 fc = pm(33)
99 rt = pm(34)
100 rc = pm(35)
101 rct1 = pm(36)
102 rct2 = pm(37)
103 aa = pm(38)
104 ac = pm(41)
105 bc = pm(39)
106 bt = pm(40)
107 rok0 = pm(29)
108 ro0 = pm(30)
109 tol = (rt-rc)/twenty
110 nindex = 0
111 DO i = 1,nel
112 scle1(i)=one
113 scle2(i)=zero
114 scle3(i)=-one
115 IF (off(i) >= one) THEN
116 nindex = nindex + 1
117 index(nindex) = i
118 rok(i) = rok0+rob(i)-ro0
119 ENDIF
120 ENDDO
121
122
123 icrit_outp = 0
124 IF (nindex > 0) THEN
125 ibug = nint(pm(59))
126 CALL frv(s1,s2,s3,s4,s5,s6,
127 . sm,vk0,vk,rob,fc,rt,rc,
128 . rct1,rct2,aa,ac,bc,bt,
129 . rok,tol,fa,nindex,index,ibug,
130 . nel,seq,icrit_outp)
131 ENDIF
132
133 nindex2 = 0
134
135#include "vectorize.inc"
136 DO n = 1, nindex
137 i = index(n)
138 IF (fa(i) < zero) THEN
139 scle3(i)=-one
140 ELSEIF(abs(fa(i)) < em10) THEN
141 scle3(i)=one
142 ELSE
143 scle3(i)=one
144 nindex2 = nindex2 + 1
145 index(nindex2) = i
146 xn(i) = one
147 ENDIF
148 ENDDO
149
150
151
152 DO nit = 1,10
153
154#include "vectorize.inc"
155 DO n = 1,nindex2
156 i = index(n)
157 IF (i > 0) THEN
158 sn1(i) = s1(i)-xn(i)*ds1(i)
159 sn2(i) = s2(i)-xn(i)*ds2(i)
160 sn3(i) = s3(i)-xn(i)*ds3(i)
161 sn4(i) = s4(i)-xn(i)*ds4(i)
162 sn5(i) = s5(i)-xn(i)*ds5(i)
163 sn6(i) = s6(i)-xn(i)*ds6(i)
164 sn7(i) = sm(i)-xn(i)*dsm(i)
165 ENDIF
166 ENDDO
167
168 CALL frv(sn1,sn2,sn3,sn4,sn5,sn6,
169 . sn7,vk0,vk,rob,fc,rt,rc,
170 . rct1,rct2,aa,ac,bc,bt,
171 . rok,tol,fn,nindex2,index,ibug,
172 . nel,seq,icrit_outp)
173
174
175#include "vectorize.inc"
176 DO n = 1,nindex2
177 i = index(n)
178 IF (i > 0) THEN
179 IF (nit==1 .AND. fn(i) > -tolf) THEN
180 scle2(i) = one
181 index(n) = 0
182 ELSE
183 scle2(i)=xn(i)/(one-fn(i)/fa(i))
184 IF (abs(fn(i)) < tolf) THEN
185 index(n) = 0
186 scle2(i) =
min(one ,scle2(i))
187 scle2(i) =
max(zero,scle2(i))
188 ELSE
189 xn(i) = scle2(i)
190 ENDIF
191 ENDIF
192 ENDIF
193 ENDDO
194
195 ENDDO
196
197
198
199
200#include "vectorize.inc"
201 DO n = 1,nindex2
202 i = index(n)
203 IF (i/=0) THEN
204 scle2(i) =
min(one,scle2(i))
205 scle2(i) =
max(zero,scle2(i))
206 ENDIF
207 ENDDO
208
209 DO i=1,nel
210 scle1(i) = one-scle2(i)
211 s1(i) = s1(i) - scle2(i)*ds1(i)
212 s2(i) = s2(i) - scle2(i)*ds2(i)
213 s3(i) = s3(i) - scle2(i)*ds3(i)
214 s4(i) = s4(i) - scle2(i)*ds4(i)
215 s5(i) = s5(i) - scle2(i)*ds5(i)
216 s6(i) = s6(i) - scle2(i)*ds6(i)
217 sm(i) = sm(i) - scle2(i)*dsm(i)
218 ENDDO
219
220 RETURN
subroutine frv(s1, s2, s3, s4, s5, s6, sm, vk0, vk, rob, fc, rt, rc, rct1, rct2, aa, ac, bc, bt, rok, tol, fa, nindex, index, ibug, nel, seq, icrit_outp)