OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigpoly.F File Reference
#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine polystress2 (nel, matb, c10, c01, c20, c11, c02, c30, c21, c12, c03, d1, d2, d3, sig, bi1, bi2, jdet, flag_mul, nvarf, coefr, betaf, coefm, uvarf, rbulk, iform)
subroutine polystrest2 (nel, matb, c10, c01, c20, c11, c02, c30, c21, c12, c03, d1, d2, d3, sig, bi1, bi2, jdet, flag_mul, nvarf, coefr, betaf, coefm, uvarf, cii, rbulk, iform)
subroutine polystress (nel, ev, nuparam, uparam, t1, t2, t3, stiff, bi1, bi2, jdet, evb)

Function/Subroutine Documentation

◆ polystress()

subroutine polystress ( integer nel,
ev,
integer nuparam,
uparam,
t1,
t2,
t3,
stiff,
bi1,
bi2,
jdet,
evb )

Definition at line 384 of file sigpoly.F.

387C-----------------------------------------------
388C I M P L I C I T T Y P E S
389C-----------------------------------------------
390#include "implicit_f.inc"
391C----------------------------------------------------------------
392C I N P U T A R G U M E N T S
393C----------------------------------------------------------------
394 INTEGER NEL , NUPARAM
395 my_real
396 . ev(nel,3),uparam(nuparam)
397C----------------------------------------------------------------
398C O U T P U T A R G U M E N T S
399C----------------------------------------------------------------
400 my_real
401 . t1(nel), t2(nel),t3(nel)
402C----------------------------------------------------------------
403C I N P U T O U T P U T A R G U M E N T S
404C----------------------------------------------------------------
405C----------------------------------------------------------------
406C L O C A L V A R I B L E S
407C----------------------------------------------------------------
408 INTEGER I
409
410 my_real
411 . lpchain(nel), trace(nel),traceb(nel),jdet(nel),stiff(nel),
412 . sb1(nel), sb2(nel),sb3(nel),tbnorm(nel),dgamma(nel),
413 . evpn(nel,3),evel(nel,3),evb(nel,3),
414 . aa,bb,cc,c10,c01,c20,c11,c02,c30,c21,c12,c03,d1,d2,d3,
415 . i2(nel),jthird(nel),j4third(nel),bi1(nel),
416 . bi2(nel),dphidi1(nel) ,dphidi2(nel) , dphidj(nel) ,inv2j(nel)
417C----------------------------------------------------------------
418 c10 = uparam(1)
419 c01 = uparam(2)
420 c20 = uparam(3)
421 c11 = uparam(4)
422 c02 = uparam(5)
423 c30 = uparam(6)
424 c21 = uparam(7)
425 c12 = uparam(8)
426 c03 = uparam(9)
427 d1 = uparam(11) !ONE/D1
428 !IF(D1 /= ZERO) D1 = ONE/D1
429 d2 = uparam(12) !ONE/D2
430 d3 = uparam(13) !ONE/D3
431
432 DO i = 1,nel
433 !J = RHO0/RHO = RELATIVE VOLUME = DET F (F = GRADIENT OF DEFORMATION)
434 jdet(i) = ev(i,1)*ev(i,2)*ev(i,3)
435 !SECOND INVARIANT
436 i2(i)= ev(i,1)*ev(i,1)*ev(i,2)*ev(i,2)
437 . +ev(i,3)*ev(i,3)*ev(i,2)*ev(i,2)
438 . +ev(i,3)*ev(i,3)*ev(i,1)*ev(i,1)
439 IF(jdet(i)>zero) THEN
440 jthird(i) = exp((-third )*log(jdet(i)))
441 j4third(i) = exp((-four_over_3)*log(jdet(i)))
442 ELSE
443 jthird(i) = zero
444 j4third(i)= zero
445 ENDIF
446 ENDDO
447 !-------------------------
448 !chain A: compute stress :
449 !-------------------------
450 DO i = 1,nel
451 !LAMBDA BAR
452 evb(i,1)=ev(i,1)*jthird(i)
453 evb(i,2)=ev(i,2)*jthird(i)
454 evb(i,3)=ev(i,3)*jthird(i)
455 !first invariant deviator BI1:
456 bi1(i) = evb(i,1)*evb(i,1)+evb(i,2)*evb(i,2)+evb(i,3)*evb(i,3)
457 !SECOND INVARIANT BAR
458 bi2(i)= evb(i,1)*evb(i,1)*evb(i,2)*evb(i,2)
459 . + evb(i,3)*evb(i,3)*evb(i,2)*evb(i,2)
460 . + evb(i,3)*evb(i,3)*evb(i,1)*evb(i,1)
461 ENDDO
462 DO i = 1,nel
463
464 dphidi1(i) = c10 + two *c20 *(bi1(i)-three)+ three*c30 *(bi1(i)-three)**2
465 . + c11 *(bi2(i)-three)+ c12 *(bi2(i)-three)**2
466 . + two *c21 *(bi1(i)-three)*(bi2(i)-three)
467
468 dphidi2(i) = c01 + two*c02*(bi2(i)-three) +three*c03*(bi2(i)-three)**2
469 . + c11*(bi1(i)-three) + c21*(bi1(i)-three)**2
470 . + two*c12*(bi1(i)-three)*(bi2(i)-three)
471 dphidj(i) = two*d1* (jdet(i)-one) + four * d2 * (jdet(i)-one)**3 + six*d3*(jdet(i)-one)**5
472
473 inv2j(i)=two/max(em20,jdet(i))
474 ENDDO
475
476 DO i=1,nel
477 !CAUCHY STRESS
478 aa = (dphidi1(i) + dphidi2(i) * bi1(i))*inv2j(i)
479 bb = dphidi2(i) *inv2j(i)
480 cc = third* inv2j(i)* ( bi1(i)* dphidi1(i)+two* bi2(i)*dphidi2(i))
481 t1(i) = aa*evb(i,1)*evb(i,1)
482 . - bb*evb(i,1)**4
483 . - cc
484 . + dphidj(i)
485 t2(i) = aa*evb(i,2)*evb(i,2)
486 . - bb*evb(i,2)**4
487 . - cc
488 . + dphidj(i)
489 t3(i) = aa*evb(i,3)*evb(i,3)
490 . - bb*evb(i,3)**4
491 . - cc
492 . + dphidj(i)
493
494! STIFF(I) =( TWO *C20 + THREE*C30 *(BI1(I)-THREE) ) !IF YEOH
495! . * ((TWO * EVB(I,1)**2 - EVB(I,2)**2 - EVB(I,3)**2)**2 )*TWO/NINE
496! . + (c10 + two *c20 *(bi1(i)-three)+ three*c30 *(bi1(i)-three)**2)
497! . * (FOUR * EVB(I,1)**2 + EVB(I,2)**2 + EVB(I,3)**2 )*FOUR/NINE
498! . + TWO*D1*( JDET(I)*(JDET(I)-ONE) + JDET(I)**2 )
499
500 ENDDO
501C
502 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21

◆ polystress2()

subroutine polystress2 ( integer, intent(in) nel,
dimension(nel,3,3), intent(in) matb,
intent(in) c10,
intent(in) c01,
intent(in) c20,
intent(in) c11,
intent(in) c02,
intent(in) c30,
intent(in) c21,
intent(in) c12,
intent(in) c03,
intent(in) d1,
intent(in) d2,
intent(in) d3,
intent(out) sig,
intent(out) bi1,
intent(out) bi2,
intent(out) jdet,
integer, intent(in) flag_mul,
integer, intent(in) nvarf,
intent(in) coefr,
intent(in) betaf,
intent(in) coefm,
intent(inout) uvarf,
rbulk,
integer, intent(in) iform )

Definition at line 36 of file sigpoly.F.

43C-----------------------------------------------
44C I M P L I C I T T Y P E S
45C-----------------------------------------------
46#include "implicit_f.inc"
47C----------------------------------------------------------------
48C I N P U T A R G U M E N T S
49C----------------------------------------------------------------
50 INTEGER, INTENT(IN) :: NEL,FLAG_MUL,NVARF,IFORM
51 my_real, INTENT(IN) :: c10,c01,c20,c11,c02,c30,c21,c12,c03,
52 . d1,d2,d3,
53 . coefr, betaf,coefm
54
55 my_real, DIMENSION(NEL, 3,3), INTENT(IN) :: matb(nel,3,3)
56C----------------------------------------------------------------
57C O U T P U T A R G U M E N T S
58C----------------------------------------------------------------
59 my_real, DIMENSION(NEL), INTENT(OUT) :: bi1,bi2,jdet
60 my_real, DIMENSION(NEL, 3,3), INTENT(OUT) :: sig
61C----------------------------------------------------------------
62C I N P U T O U T P U T A R G U M E N T S
63C----------------------------------------------------------------
64 my_real, DIMENSION(NEL,NVARF), INTENT(INOUT) :: uvarf
65C----------------------------------------------------------------
66C L O C A L V A R I B L E S
67C----------------------------------------------------------------
68 INTEGER I
69
70 my_real
71 . wmax,lpchain(nel), trace(nel),traceb(nel),j2third(nel),
72 . sb1(nel), sb2(nel),sb3(nel),tbnorm(nel),dgamma(nel),i1(nel),
73 . aa,bb,cc,trb2,trb22,rbulk,
74 . i2(nel),jthird(nel),j4third(nel),eta(nel),ww(nel) ,
75 . dphidi1(nel) ,dphidi2(nel) , dphidj(nel) ,inv2j(nel),
76 . matb2(nel,3,3)
77C----------------------------------------------------------------
78 eta(1:nel) = one !MULLINS DAMAGE FACTOR
79
80 CALL prodmat (matb , matb, matb2, nel) ! MATB2 = MATB**2
81 DO i = 1,nel
82 !J = RHO0/RHO = RELATIVE VOLUME = DET F (F = GRADIENT OF DEFORMATION)
83 jdet(i)=matb(i,1,1)*matb(i,2,2)*matb(i,3,3)-matb(i,1,1)*matb(i,2,3)*matb(i,3,2) -
84 . matb(i,3,3)*matb(i,1,2)*matb(i,2,1) +matb(i,1,2)*matb(i,2,3)*matb(i,3,1) +
85 . matb(i,2,1)*matb(i,3,2)*matb(i,1,3) -matb(i,2,2)*matb(i,3,1)*matb(i,1,3)
86 jdet(i)= sqrt(max(em20, jdet(i)))
87 !SECOND INVARIANT
88 i1(i) = matb(i,1,1)+matb(i,2,2)+matb(i,3,3)
89
90 trb2 = i1(i)**2
91 trb22= matb2(i,1,1) + matb2(i,2,2) +matb2(i,3,3)
92
93 !I2 = ((TrB)**2 - Tr (B**2) )/2
94 i2(i)= (trb2 - trb22)/two
95
96 IF(jdet(i)>zero) THEN
97 jthird(i) = exp((-third )*log(jdet(i)))
98 j2third(i) = jthird(i)**2
99 j4third(i) = jthird(i)**4
100 ELSE
101 jthird(i) = zero
102 j2third(i) = zero
103 j4third(i) = zero
104 ENDIF
105 ENDDO
106 DO i = 1,nel
107 !INVARIANT BAR
108 !first invariant deviator BI1:
109 bi1(i) = i1(i) * j2third(i)
110 !SECOND INVARIANT BAR
111 bi2(i) = i2(i)*j4third(i)
112 ENDDO
113 wmax = zero
114
115 DO i = 1,nel ! helmotz free energy for polynomial
116 ww(i) = c10 *(bi1(i)-three) + c20 *(bi1(i)-three)**2 + c30 *(bi1(i)-three)**3
117 . + c01 *(bi2(i)-three) + c02 *(bi2(i)-three)**2 + c03 *(bi2(i)-three)**3
118 . + c11 *(bi1(i)-three)*(bi2(i)-three)
119 . + c12 *(bi1(i)-three)*(bi2(i)-three)**2
120 . + c21 *(bi2(i)-three)*(bi1(i)-three)**2
121 wmax = max(wmax, ww(i)) ! need deviatoric part for Mullins
122 ENDDO
123
124 !====================================
125 !dammge by mullins effect
126 !====================================
127 IF(flag_mul == 1)THEN
128 CALL mullins_or(
129 1 nel ,nvarf, coefr,betaf ,
130 2 coefm, ww , uvarf,eta )
131 ENDIF
132 !====================================
133
134 DO i = 1,nel
135 eta(i) = max(min(eta(i),one),em20)
136 ! derivatives of energy function
137 dphidi1(i) = eta(i)*(c10 + two *c20 *(bi1(i)-three)+ three*c30 *(bi1(i)-three)**2
138 . + c11 *(bi2(i)-three)+ c12 *(bi2(i)-three)**2
139 . + two *c21 *(bi1(i)-three)*(bi2(i)-three))
140
141 dphidi2(i) = eta(i)*(c01 + two*c02*(bi2(i)-three) +three*c03*(bi2(i)-three)**2
142 . + c11*(bi1(i)-three) + c21*(bi1(i)-three)**2
143 . + two*c12*(bi1(i)-three)*(bi2(i)-three))
144
145 inv2j(i)=two/max(em20,jdet(i))
146 ENDDO
147 IF (iform == 1) THEN
148 DO i = 1,nel
149 dphidj(i) = two*d1* (jdet(i)-one) + four * d2 * (jdet(i)-one)**3 + six*d3*(jdet(i)-one)**5
150 ENDDO
151 ELSEIF (iform == 2) THEN
152 DO i = 1,nel
153 dphidj(i) = rbulk * (one - one / jdet(i))
154 ENDDO
155 ENDIF
156
157
158 DO i=1,nel
159 !CAUCHY STRESS
160 aa = (dphidi1(i) + dphidi2(i) * bi1(i))*inv2j(i)*j2third(i)
161 bb = dphidi2(i) *inv2j(i)*j4third(i)
162 cc = third* inv2j(i)* ( bi1(i)* dphidi1(i)+two* bi2(i)*dphidi2(i))
163
164 sig(i,1,1) = aa*matb(i,1,1)
165 . - bb*matb2(i,1,1)
166 . - cc
167 . + dphidj(i)
168 sig(i,2,2) = aa*matb(i,2,2)
169 . - bb*matb2(i,2,2)
170 . - cc
171 . + dphidj(i)
172 sig(i,3,3) = aa*matb(i,3,3)
173 . - bb*matb2(i,3,3)
174 . - cc
175 . + dphidj(i)
176 sig(i,1,2) = aa*matb(i,1,2)
177 . - bb*matb2(i,1,2)
178 sig(i,2,3) = aa*matb(i,2,3)
179 . - bb*matb2(i,2,3)
180 sig(i,3,1) = aa*matb(i,3,1)
181 . - bb*matb2(i,3,1)
182 sig(i,2,1)=sig(i,1,2)
183 sig(i,3,2)=sig(i,2,3)
184 sig(i,1,3)=sig(i,3,1)
185 ENDDO
186C
187 RETURN
mullins_or
subroutine mullins_or(nel, nuvar, coefr, betaf, coefm, ww, uvar, eta)
Definition fail_mullins_OR_s.F:130
min
#define min(a, b)
Definition macros.h:20
prodmat
subroutine prodmat(a, b, c, nel)
Definition prodmat.F:35

◆ polystrest2()

subroutine polystrest2 ( integer, intent(in) nel,
dimension(nel,3,3), intent(in) matb,
intent(in) c10,
intent(in) c01,
intent(in) c20,
intent(in) c11,
intent(in) c02,
intent(in) c30,
intent(in) c21,
intent(in) c12,
intent(in) c03,
intent(in) d1,
intent(in) d2,
intent(in) d3,
intent(out) sig,
intent(out) bi1,
intent(out) bi2,
intent(out) jdet,
integer, intent(in) flag_mul,
integer, intent(in) nvarf,
intent(in) coefr,
intent(in) betaf,
intent(in) coefm,
intent(inout) uvarf,
intent(out) cii,
intent(in) rbulk,
integer, intent(in) iform )

Definition at line 197 of file sigpoly.F.

204C-----------------------------------------------
205C I M P L I C I T T Y P E S
206C-----------------------------------------------
207#include "implicit_f.inc"
208C----------------------------------------------------------------
209C I N P U T A R G U M E N T S
210C----------------------------------------------------------------
211 INTEGER, INTENT(IN) :: NEL,FLAG_MUL,NVARF,IFORM
212 my_real, INTENT(IN) :: c10,c01,c20,c11,c02,c30,c21,c12,c03,
213 . d1,d2,d3,
214 . coefr, betaf,coefm ,rbulk
215
216 my_real, DIMENSION(NEL, 3,3), INTENT(IN) :: matb(nel,3,3)
217C----------------------------------------------------------------
218C O U T P U T A R G U M E N T S
219C----------------------------------------------------------------
220 my_real, DIMENSION(NEL), INTENT(OUT) :: bi1,bi2,jdet
221 my_real, DIMENSION(NEL, 3,3), INTENT(OUT) :: sig
222 my_real, DIMENSION(NEL, 3), INTENT(OUT) :: cii
223C----------------------------------------------------------------
224C I N P U T O U T P U T A R G U M E N T S
225C----------------------------------------------------------------
226 my_real, DIMENSION(NEL,NVARF), INTENT(INOUT) :: uvarf
227C----------------------------------------------------------------
228C L O C A L V A R I B L E S
229C----------------------------------------------------------------
230 INTEGER I
231
232 my_real
233 . wmax,lpchain(nel), trace(nel),traceb(nel),j2third(nel),matb2(nel,3,3),
234 . sb1(nel), sb2(nel),sb3(nel),tbnorm(nel),dgamma(nel),i1(nel),
235 . aa,bb,cc,trb2,trb22,
236 . i2(nel),jthird(nel),j4third(nel),eta(nel),ww(nel) ,
237 . dphidi1(nel) ,dphidi2(nel) , dphidj(nel) ,inv2j(nel),
238 . dphi2di1(nel) ,dphi2di2(nel) , dphi2dj(nel),lam_b(3),lam_b_1(3),
239 . bi1_3,bi2_3
240C----------------------------------------------------------------
241 eta(1:nel) = one !MULLINS DAMAGE FACTOR
242
243 CALL prodmat (matb , matb, matb2, nel) ! MATB2 = MATB**2
244
245 DO i = 1,nel
246 !J = RHO0/RHO = RELATIVE VOLUME = DET F (F = GRADIENT OF DEFORMATION)
247 jdet(i)=matb(i,1,1)*matb(i,2,2)*matb(i,3,3)-matb(i,1,1)*matb(i,2,3)*matb(i,3,2) -
248 . matb(i,3,3)*matb(i,1,2)*matb(i,2,1) +matb(i,1,2)*matb(i,2,3)*matb(i,3,1) +
249 . matb(i,2,1)*matb(i,3,2)*matb(i,1,3) -matb(i,2,2)*matb(i,3,1)*matb(i,1,3)
250 jdet(i)= sqrt(max(em20, jdet(i)))
251
252 !SECOND INVARIANT
253 i1(i) = matb(i,1,1)+matb(i,2,2)+matb(i,3,3)
254
255 trb2 = i1(i)**2
256 trb22= matb2(i,1,1) + matb2(i,2,2) +matb2(i,3,3)
257 !I2 = ((TrB)**2 - Tr (B**2) )/2
258 i2(i)= (trb2 - trb22)/two
259
260 IF(jdet(i)>zero) THEN
261 jthird(i) = exp((-third )*log(jdet(i)))
262 j2third(i) = jthird(i)**2
263 j4third(i) = jthird(i)**4
264 ELSE
265 jthird(i) = zero
266 j2third(i) = zero
267 j4third(i) = zero
268 ENDIF
269 ENDDO
270 DO i = 1,nel
271 !INVARIANT BAR
272 !first invariant deviator BI1:
273 bi1(i) = i1(i) * j2third(i)
274 !SECOND INVARIANT BAR
275 bi2(i) = i2(i)*j4third(i)
276 ENDDO
277 wmax = zero
278
279 DO i = 1,nel ! helmotz free energy for polynomial
280 ww(i) = c10 *(bi1(i)-three) + c20 *(bi1(i)-three)**2 + c30 *(bi1(i)-three)**3
281 . + c01 *(bi2(i)-three) + c02 *(bi2(i)-three)**2 + c03 *(bi2(i)-three)**3
282 . + c11 *(bi1(i)-three)*(bi2(i)-three)
283 . + c12 *(bi1(i)-three)*(bi2(i)-three)**2
284 . + c21 *(bi2(i)-three)*(bi1(i)-three)**2
285 wmax = max(wmax, ww(i)) ! need deviatoric part for Mullins
286 ENDDO
287
288 !====================================
289 !dammge by mullins effect
290 !====================================
291 IF(flag_mul == 1)THEN
292 CALL mullins_or(
293 1 nel ,nvarf, coefr,betaf ,
294 2 coefm, ww , uvarf,eta )
295 ENDIF
296 !====================================
297
298 DO i = 1,nel
299 eta(i) = max(min(eta(i),one),em20)
300 ! derivatives of energy function
301 dphidi1(i) = eta(i)*(c10 + two *c20 *(bi1(i)-three)+ three*c30 *(bi1(i)-three)**2
302 . + c11 *(bi2(i)-three)+ c12 *(bi2(i)-three)**2
303 . + two *c21 *(bi1(i)-three)*(bi2(i)-three))
304
305 dphidi2(i) = eta(i)*(c01 + two*c02*(bi2(i)-three) +three*c03*(bi2(i)-three)**2
306 . + c11*(bi1(i)-three) + c21*(bi1(i)-three)**2
307 . + two*c12*(bi1(i)-three)*(bi2(i)-three))
308 inv2j(i)=two/max(em20,jdet(i))
309 ENDDO
310 IF (iform == 1) THEN
311 DO i = 1,nel
312 dphidj(i) = two*d1* (jdet(i)-one) + four * d2 * (jdet(i)-one)**3 + six*d3*(jdet(i)-one)**5
313 ENDDO
314 ELSEIF (iform == 2) THEN
315 DO i = 1,nel
316 dphidj(i) = rbulk * (one - one / jdet(i))
317 ENDDO
318 ENDIF
319
320 DO i=1,nel
321 !CAUCHY STRESS
322 aa = (dphidi1(i) + dphidi2(i) * bi1(i))*inv2j(i)*j2third(i)
323 bb = dphidi2(i) *inv2j(i)*j4third(i)
324 cc = third* inv2j(i)* ( bi1(i)* dphidi1(i)+two* bi2(i)*dphidi2(i))
325
326 sig(i,1,1) = aa*matb(i,1,1)
327 . - bb*matb2(i,1,1)
328 . - cc
329 . + dphidj(i)
330 sig(i,2,2) = aa*matb(i,2,2)
331 . - bb*matb2(i,2,2)
332 . - cc
333 . + dphidj(i)
334 sig(i,3,3) = aa*matb(i,3,3)
335 . - bb*matb2(i,3,3)
336 . - cc
337 . + dphidj(i)
338 sig(i,1,2) = aa*matb(i,1,2)
339 . - bb*matb2(i,1,2)
340 sig(i,2,3) = aa*matb(i,2,3)
341 . - bb*matb2(i,2,3)
342 sig(i,3,1) = aa*matb(i,3,1)
343 . - bb*matb2(i,3,1)
344 sig(i,2,1)=sig(i,1,2)
345 sig(i,3,2)=sig(i,2,3)
346 sig(i,1,3)=sig(i,3,1)
347
348 ENDDO
349
350 IF (iform == 1) THEN
351 DO i = 1,nel
352 dphi2dj(i) = two*d1 + twelve * d2 * (jdet(i)-one)**2 + thirty*d3*(jdet(i)-one)**4
353 ENDDO
354 ELSEIF (iform == 2) THEN
355 DO i = 1,nel
356 dphi2dj(i) = rbulk/(jdet(i)**2)
357 ENDDO
358 ENDIF
359 DO i = 1,nel
360 dphi2di1(i) = two*eta(i)*(c20 + three*c30*(bi1(i)-three)+ c21*(bi2(i)-three))
361
362 dphi2di2(i) = two*eta(i)*(c02 +three*c03*(bi2(i)-three) + c12*(bi1(i)-three))
363
364 lam_b(1) = matb(i,1,1)*j2third(i)
365 lam_b(2) = matb(i,2,2)*j2third(i)
366 lam_b(3) = matb(i,3,3)*j2third(i)
367 lam_b_1(1:3) = one/lam_b(1:3)
368 bi1_3 = third*bi1(i)
369 bi2_3 = third*bi2(i)
370 cii(i,1:3) = two*(two_third*dphidi1(i)*(lam_b(1:3)+bi1_3)+dphi2di1(i)*(lam_b(1:3)-bi1_3)+
371 . two_third*dphidi2(i)*(lam_b_1(1:3)+bi2_3)+dphi2di2(i)*(lam_b_1(1:3)-bi2_3))
372 . + dphi2dj(i)
373 ENDDO
374C
375 RETURN