!************************************************************** !* AceGen 1.1 Windows (13 Aug 07) * !* Co. J. Korelc 2007 13 Aug 07 15:09:52* !************************************************************** ! User : USER ! Evaluation time : 13 s Mode : Optimal ! Number of formulae : 174 Method: Automatic ! Subroutine : SKR2999 size :4863 ! Total size of Mathematica code : 4863 subexpressions ! Total size of Fortran code : 16096 bytes SUBROUTINE SMI2999 & (switch,lesvbl,ISENMR,mswitch,meuvbl,nehist,jfile,nstate, & state,mstate,order,morder) implicit none integer mstate,mgdata,mgpost,mnpost,isenmr parameter(mnpost=1000,mgpost=1000,mgdata=1000) include 'sms.h' integer mswitch,switch(mswitch),meuvbl,lesvbl(meuvbl),nehist integer morder,norder,order(morder),ordert(8),jfile,i integer idata(IData_Length),ngpo integer ngspost,nnpost,nnspost,nstate,ngdata,ngpost character*50 SELEM character*8 gdata(mgdata),gpost(mgpost) character*8 npost(mnpost),state(mstate) logical DEBUG parameter (DEBUG=.false., & SELEM="ExamplesHypersolidELFEN") data (ordert(i),i=1,8) /1,2,3,4,5,6,7,8/ do i = 1,8 order(i) = ordert(i) enddo C switch - Array of the input data for SMS C switch(1) NGAUSS - Number of gauss points switch(1)=8 C switch(2) NSENPA - Number of sensitivity parameters C lesvbl - Array of the element variables vlues C lesvbl(1) NANODE - Number of active nodes per element lesvbl(1)=8 C lesvbl(2) NPNODE - Number of plot nodes per element lesvbl(2)=8 C lesvbl(3) NATYPE - Analysis type lesvbl(3)=1 C lesvbl(4) NDOFND - Default number of degrees of freedom per node lesvbl(4)=3 C lesvbl(5) NELMOD - Element module number lesvbl(5)=1 C lesvbl(6) NEVAB - Number of variables lesvbl(6)=24 C lesvbl(7) NEVABF - Number of variables of fluid element lesvbl(7)=0 C lesvbl(8) NFACE - Number of faces per element lesvbl(8)=6 C lesvbl(9) NFNOD - Maximum number of nodes per face lesvbl(9)=4 C lesvbl(10) NFNODF - Maximum number of nodes per face for fluid element lesvbl(10)=0 C lesvbl(11) NGAUSA - Alternative ( OPT(56) ) number of Gauss points C used to evaluate the stiffness matrix lesvbl(11)=switch(1) C lesvbl(12) NGAUSB - Number of Gauss points used to evaluate the body forces lesvbl(12)=switch(1) C lesvbl(13) NGAUSM - Number of Gauss points used to evaluate the mass matrix lesvbl(13)=switch(1) C lesvbl(14) NGAUSS - Number of Gauss points used to evaluate the stiffness matrix lesvbl(14)=switch(1) C lesvbl(15) NGAUSF - Number of Gauss points used to evaluate H and S matrices lesvbl(15)=switch(1) C lesvbl(16) NGPRP - Number of geometric properties required per element lesvbl(16)=0 C lesvbl(17) NNODE - Number of nodes per element lesvbl(17)=8 C lesvbl(18) NNODEF - Number of nodes per fluid element lesvbl(18)=0 C lesvbl(19) NSFR - Shape parameter lesvbl(19)=8 C lesvbl(20) NSFR - Shape parameter for fluid element lesvbl(20)=0 C lesvbl(21) NSTRE - Number of stresses lesvbl(21)=6 C lesvbl(22) NSTRAN - Number of strains lesvbl(22)=6 C lesvbl(23) NEDNOD - Number of nodes on each edge lesvbl(23)=2 C lesvbl(24) NSTTYP - Stress type number lesvbl(24)=1 C lesvbl(25) NHGSSF - Number of hourglass forces lesvbl(25)=0 C lesvbl(26) NEVABA - Number of variables for coupled problems lesvbl(26)=0 C lesvbl(27) NFTYP - Face type C 1 - Line 2 - Triangle 3 - Serendipity Quadrilateral 4 - Lagrange Quadrilateral lesvbl(27)=3 C lesvbl(28) standard integration code (SMS) lesvbl(28)=7 idata(ID_NoSensParameters)=switch(2) ngpo=switch(1) idata(ID_MaxMessages)=10000 idata(ID_NoMessages)=0 write(jfile,*)"INITIALIZE ELEMENT : ",SELEM c.....Description of the state variables nstate=0 nehist=0 c.....Description of the input data ngdata=2 write(jfile,"(A)")"Description of the element group input data:" gdata(1)="Elastic modulus" gdata(2)="Poisson ratio" write(jfile,"(i10,A3,A)") & (i," = ",gdata(i),i=1,2) c.... description of the postprocessing data if(DEBUG) then write(*,*)"INITIALIZE ELEMENT : ",SELEM endif if(nehist.eq.0) nehist=1 if(nstate.eq.0) nstate=1 RETURN END subroutine SMS2999(gdata,xlu,ul,ul0,sig,sig0,eps,eps0,sta,sta0,h, & sxd,s,p,ielf,nehist,nhsets,mgdata,jfile,ngpo) implicit none integer ielf(10) integer mgpost,mnpost parameter(mgpost=1000,mnpost=1000) include 'sms.h' integer nehist,nhsets,jfile,idata(IData_Length) integer ngpo,ifix,hlen,mgdata,i,j,icode DOUBLE PRECISION xlu(3,8),xl(3,8) DOUBLE PRECISION gdata(mgdata),ul(48),ul0(48),sxd(24) DOUBLE PRECISION sig(6,ngpo),sig0(6,ngpo),eps(6,ngpo),eps0(6,ngpo) DOUBLE PRECISION sta(0,ngpo),sta0(0,ngpo) DOUBLE PRECISION h(nehist,nhsets),s(24,24),p(24) DOUBLE PRECISION gpost(ngpo,mgpost),npost(8,mnpost) DOUBLE PRECISION gp(4,100),v(5305),rdata(RData_Length) character*50 SELEM logical DEBUG parameter (DEBUG=.false., & SELEM="ExamplesHypersolidELFEN") 1235 format(i3,">",3g17.9) 1236 format(" >",5g17.9) 1237 format(i3,"> ",24g11.5) 1238 format(i3,"> ",24g11.5) 1239 format(10x,f15.5," = ",A) idata(ID_NoSensParameters)=ielf(2) idata(ID_SensIndex)=ielf(4) idata(ID_SensType)=ielf(3) idata(ID_SensTypeIndex)=ielf(5) idata(ID_ErrorStatus)=ielf(6) rdata(RD_SubIterationTolerance)=1d-9 if(ielf(1).eq.1)then idata(ID_Task)=TASK_Tangent elseif(ielf(1).eq.2)then idata(ID_Task)=TASK_Residual elseif(ielf(1).eq.3)then idata(ID_Task)=TASK_SensitivityPseudoLoad elseif(ielf(1).eq.4)then idata(ID_Task)=TASK_DependentSensitivity endif idata(ID_OutputFile)=jfile idata(ID_Iteration)=2 idata(ID_GlobalIterationMode)=0 do i=1,3 do j=1,8 xl(i,j)=xlu(i,j)-ul((j-1)*3+i) enddo enddo icode=7 jfile=idata(ID_OutputFile) call SMSIntPoints(icode,ngpo,gp) if(DEBUG) then write(*,*) SELEM ," task=",idata(ID_Task) write(jfile,*)SELEM c.....Description of the input data write(jfile,1239)gdata(1), & "Elastic modulus" write(jfile,1239)gdata(2), & "Poisson ratio" hlen=0 endif idata(ID_PostIteration)=0 select case(idata(ID_Task)) c.....presolution case (TASK_Tangent) idata(ID_SkipTangent)=0 idata(ID_SkipResidual)=1 idata(ID_SubIterationMode)=1; if(DEBUG) then write(jfile,*)"#################################################" write(jfile,*)"TANGENT for ",SELEM write(jfile,*)"NODE X" write(jfile,1235)(i,(xl(j,i),j=1,3),i=1,8) write(jfile,*)"NODE ap" write(jfile,1235)(i,(ul0((i-1)*3+j),j=1,3),i=1,8) write(jfile,*)"NODE at" write(jfile,1235)(i,(ul((i-1)*3+j),j=1,3),i=1,8) write(jfile,*)"ELEMENT hp" write(jfile,1236)(h(i,2),i=1,hlen) do j=1,ngpo write(jfile,*)"INTEGRATION POINT ",j write(jfile,*)"STRESS-p" write(jfile,1236)(sig0(i,j),i=1,6) write(jfile,*)"STRAIN-p" write(jfile,1236)(eps0(i,j),i=1,6) write(jfile,*)"STATE-p" write(jfile,1236)(sta0(i,j),i=1,0) enddo endif call SMSZero(s,576) call SKR2999(v,gdata,rdata,xl,ul,ul0,sig,sig0,eps, & eps0,sta,sta0,s,p,h,gp,idata,nehist,nhsets,ngpo) do i=2,24 do j=1,i-1 s(i,j)=s(j,i) enddo enddo if(DEBUG) then write(jfile,*)"ELEMENT ht" write(jfile,1236)(h(i,1),i=1,hlen) do j=1,ngpo write(jfile,*)"INTEGRATION POINT ",j write(jfile,*)"STRESS-t" write(jfile,1236)(sig(i,j),i=1,6) write(jfile,*)"STRAIN-t" write(jfile,1236)(eps(i,j),i=1,6) write(jfile,*)"STATE-t" write(jfile,1236)(sta(i,j),i=1,0) enddo write(jfile,*)"MatrixLocal=" write(jfile,1237)( i,(s(i,j),j=1,24),i=1,24) endif c.....postsolution case (TASK_Residual) idata(ID_SkipTangent)=1 idata(ID_SkipResidual)=0 idata(ID_SubIterationMode)=0; if(DEBUG) then write(jfile,*)"#################################################" write(jfile,*)"RESIDUAL for ",SELEM write(jfile,*)"NODE X" write(jfile,1235)(i,(xl(j,i),j=1,3),i=1,8) write(jfile,*)"NODE ap" write(jfile,1235)(i,(ul0((i-1)*3+j),j=1,3),i=1,8) write(jfile,*)"NODE at" write(jfile,1235)(i,(ul((i-1)*3+j),j=1,3),i=1,8) write(jfile,*)"ELEMENT hp" write(jfile,1236)(h(i,2),i=1,hlen) do j=1,ngpo write(jfile,*)"INTEGRATION POINT ",j write(jfile,*)"STRESS-p" write(jfile,1236)(sig0(i,j),i=1,6) write(jfile,*)"STRAIN-p" write(jfile,1236)(eps0(i,j),i=1,6) write(jfile,*)"STATE-p" write(jfile,1236)(sta0(i,j),i=1,0) enddo endif call SMSZero(p,24) call SKR2999(v,gdata,rdata,xl,ul,ul0,sig,sig0,eps, & eps0,sta,sta0,s,p,h,gp,idata,nehist,nhsets,ngpo) if(DEBUG) then write(jfile,*)"ELEMENT ht" write(jfile,1236)(h(i,1),i=1,hlen) do j=1,ngpo write(jfile,*)"INTEGRATION POINT ",j write(jfile,*)"STRESS-t" write(jfile,1236)(sig(i,j),i=1,6) write(jfile,*)"STRAIN-t" write(jfile,1236)(eps(i,j),i=1,6) write(jfile,*)"STATE-t" write(jfile,1236)(sta(i,j),i=1,0) enddo write(jfile,*)"VectorLocal=" write(jfile,1237)1,(p(i),i=1,24) endif c.....sensitivity analysis external part case(TASK_SensitivityPseudoLoad) c.....sensitivity analysis internal part case (TASK_DependentSensitivity) case default write(*,*)"WARNING: Unsupported task",idata(ID_Task) idata(ID_ErrorStatus)=ERROR_Fatal end select ielf(6)=idata(ID_ErrorStatus) End !******************* S U B R O U T I N E ********************** SUBROUTINE SKR2999(v,d,rdata,xl,ul,ul0,sig,sig0,eps,eps0,sta &,sta0,s,p,h,gp,idata,nenshs,nssets,ngpo) IMPLICIT NONE include 'sms.h' INTEGER idata(IData_Length),nenshs,nssets,ngpo,i1,i207,i223 &,i238 LOGICAL b206 DOUBLE PRECISION v(5305),d(2),rdata(RData_Length),xl(3,8),ul(48 &),ul0(48),sig(6,ngpo),sig0(6,ngpo),eps(6,ngpo),eps0(6,ngpo),sta &(0,ngpo),sta0(0,ngpo),s(24,24),p(24),h(nenshs,nssets),gp(4,ngpo &) v(271)=2d0*d(2) v(203)=1d0/(1d0+d(2)) v(204)=(d(1)*v(203))/2d0 v(202)=(v(204)*v(271))/(1d0-v(271)) DO i1=1,int(ngpo) v(52)=gp(1,i1) v(63)=1d0-v(52) v(59)=1d0+v(52) v(53)=gp(2,i1) v(69)=1d0+v(53) v(78)=v(69)/8d0 v(86)=-(v(63)*v(78)) v(85)=-(v(59)*v(78)) v(65)=1d0-v(53) v(79)=v(65)/8d0 v(84)=-(v(59)*v(79)) v(82)=-(v(63)*v(79)) v(54)=gp(3,i1) v(70)=1d0+v(54) v(77)=v(70)/8d0 v(88)=-(v(59)*v(77)) v(87)=-(v(63)*v(77)) v(72)=v(69)*v(77) v(67)=v(65)*v(77) v(60)=1d0-v(54) v(80)=v(60)/8d0 v(83)=-(v(59)*v(80)) v(81)=-(v(63)*v(80)) v(62)=v(69)*v(80) v(57)=v(65)*v(80) v(55)=gp(4,i1) v(89)=v(57)*(-xl(1,1)+xl(1,2))+v(62)*(xl(1,3)-xl(1,4))+v(67)*( & -xl(1,5)+xl(1,6))+v(72)*(xl(1,7)-xl(1,8)) v(90)=v(57)*(-xl(2,1)+xl(2,2))+v(62)*(xl(2,3)-xl(2,4))+v(67)*( & -xl(2,5)+xl(2,6))+v(72)*(xl(2,7)-xl(2,8)) v(91)=v(57)*(-xl(3,1)+xl(3,2))+v(62)*(xl(3,3)-xl(3,4))+v(67)*( & -xl(3,5)+xl(3,6))+v(72)*(xl(3,7)-xl(3,8)) v(92)=v(83)*(xl(1,2)-xl(1,3))+v(81)*(xl(1,1)-xl(1,4))+v(88)* & (xl(1,6)-xl(1,7))+v(87)*(xl(1,5)-xl(1,8)) v(93)=v(83)*(xl(2,2)-xl(2,3))+v(81)*(xl(2,1)-xl(2,4))+v(88)* & (xl(2,6)-xl(2,7))+v(87)*(xl(2,5)-xl(2,8)) v(274)=-(v(90)*v(92))+v(89)*v(93) v(94)=v(83)*(xl(3,2)-xl(3,3))+v(81)*(xl(3,1)-xl(3,4))+v(88)* & (xl(3,6)-xl(3,7))+v(87)*(xl(3,5)-xl(3,8)) v(273)=v(91)*v(92)-v(89)*v(94) v(272)=-(v(91)*v(93))+v(90)*v(94) v(95)=v(82)*(xl(1,1)-xl(1,5))+v(84)*(xl(1,2)-xl(1,6))+v(85)* & (xl(1,3)-xl(1,7))+v(86)*(xl(1,4)-xl(1,8)) v(96)=v(82)*(xl(2,1)-xl(2,5))+v(84)*(xl(2,2)-xl(2,6))+v(85)* & (xl(2,3)-xl(2,7))+v(86)*(xl(2,4)-xl(2,8)) v(97)=v(82)*(xl(3,1)-xl(3,5))+v(84)*(xl(3,2)-xl(3,6))+v(85)* & (xl(3,3)-xl(3,7))+v(86)*(xl(3,4)-xl(3,8)) v(98)=v(272)*v(95)+v(273)*v(96)+v(274)*v(97) v(275)=v(55)*v(98) v(102)=(-(v(94)*v(96))+v(93)*v(97))/v(98) v(165)=v(102)*v(72) v(156)=-(v(102)*v(67)) v(147)=v(102)*v(62) v(138)=-(v(102)*v(57)) v(103)=(v(94)*v(95)-v(92)*v(97))/v(98) v(167)=v(103)*v(72) v(158)=-(v(103)*v(67)) v(149)=v(103)*v(62) v(140)=-(v(103)*v(57)) v(104)=(-(v(93)*v(95))+v(92)*v(96))/v(98) v(169)=v(104)*v(72) v(160)=-(v(104)*v(67)) v(151)=v(104)*v(62) v(142)=-(v(104)*v(57)) v(105)=(v(91)*v(96)-v(90)*v(97))/v(98) v(126)=v(105)*v(80) v(114)=v(105)*v(77) v(106)=(-(v(91)*v(95))+v(89)*v(97))/v(98) v(129)=v(106)*v(80) v(116)=v(106)*v(77) v(107)=(v(90)*v(95)-v(89)*v(96))/v(98) v(132)=v(107)*v(80) v(118)=v(107)*v(77) v(108)=v(272)/v(98) v(127)=v(108)*v(79) v(120)=v(108)*v(78) v(109)=v(273)/v(98) v(130)=v(109)*v(79) v(122)=v(109)*v(78) v(110)=v(274)/v(98) v(133)=v(110)*v(79) v(124)=v(110)*v(78) v(111)=v(114)+v(120) v(112)=v(116)+v(122) v(113)=v(118)+v(124) v(115)=-v(114)+v(127) v(117)=-v(116)+v(130) v(119)=-v(118)+v(133) v(121)=-v(120)+v(126) v(123)=-v(122)+v(129) v(125)=-v(124)+v(132) v(128)=-v(126)-v(127) v(131)=-v(129)-v(130) v(134)=-v(132)-v(133) v(135)=v(138)+v(128)*v(63) v(136)=v(140)+v(131)*v(63) v(137)=v(142)+v(134)*v(63) v(139)=-v(138)+v(128)*v(59) v(141)=-v(140)+v(131)*v(59) v(143)=-v(142)+v(134)*v(59) v(144)=v(147)+v(121)*v(59) v(145)=v(149)+v(123)*v(59) v(146)=v(151)+v(125)*v(59) v(148)=-v(147)+v(121)*v(63) v(150)=-v(149)+v(123)*v(63) v(152)=-v(151)+v(125)*v(63) v(153)=v(156)+v(115)*v(63) v(154)=v(158)+v(117)*v(63) v(155)=v(160)+v(119)*v(63) v(157)=-v(156)+v(115)*v(59) v(159)=-v(158)+v(117)*v(59) v(161)=-v(160)+v(119)*v(59) v(162)=v(165)+v(111)*v(59) v(163)=v(167)+v(112)*v(59) v(164)=v(169)+v(113)*v(59) v(166)=-v(165)+v(111)*v(63) v(168)=-v(167)+v(112)*v(63) v(170)=-v(169)+v(113)*v(63) v(172)=ul(2)*v(135)+ul(5)*v(139)+ul(8)*v(144)+ul(11)*v(148)+ul & (14)*v(153)+ul(17)*v(157)+ul(20)*v(162)+ul(23)*v(166) v(173)=ul(3)*v(135)+ul(6)*v(139)+ul(9)*v(144)+ul(12)*v(148)+ul & (15)*v(153)+ul(18)*v(157)+ul(21)*v(162)+ul(24)*v(166) v(174)=ul(1)*v(136)+ul(4)*v(141)+ul(7)*v(145)+ul(10)*v(150)+ul & (13)*v(154)+ul(16)*v(159)+ul(19)*v(163)+ul(22)*v(168) v(176)=ul(3)*v(136)+ul(6)*v(141)+ul(9)*v(145)+ul(12)*v(150)+ul & (15)*v(154)+ul(18)*v(159)+ul(21)*v(163)+ul(24)*v(168) v(177)=ul(1)*v(137)+ul(4)*v(143)+ul(7)*v(146)+ul(10)*v(152)+ul & (13)*v(155)+ul(16)*v(161)+ul(19)*v(164)+ul(22)*v(170) v(178)=ul(2)*v(137)+ul(5)*v(143)+ul(8)*v(146)+ul(11)*v(152)+ul & (14)*v(155)+ul(17)*v(161)+ul(20)*v(164)+ul(23)*v(170) v(192)=1d0+ul(1)*v(135)+ul(4)*v(139)+ul(7)*v(144)+ul(10)*v(148 & )+ul(13)*v(153)+ul(16)*v(157)+ul(19)*v(162)+ul(22)*v(166) v(254)=v(172)*v(177)-v(178)*v(192) v(252)=v(173)*v(174)-v(176)*v(192) v(193)=1d0+ul(2)*v(136)+ul(5)*v(141)+ul(8)*v(145)+ul(11)*v(150 & )+ul(14)*v(154)+ul(17)*v(159)+ul(20)*v(163)+ul(23)*v(168) v(256)=v(174)*v(178)-v(177)*v(193) v(253)=v(172)*v(176)-v(173)*v(193) v(249)=-(v(172)*v(174))+v(192)*v(193) v(194)=1d0+ul(3)*v(137)+ul(6)*v(143)+ul(9)*v(146)+ul(12)*v(152 & )+ul(15)*v(155)+ul(18)*v(161)+ul(21)*v(164)+ul(24)*v(170) v(257)=v(176)*v(177)-v(174)*v(194) v(255)=v(173)*v(178)-v(172)*v(194) v(251)=-(v(176)*v(178))+v(193)*v(194) v(250)=-(v(173)*v(177))+v(192)*v(194) v(201)=v(194)*v(249)+v(178)*v(252)+v(177)*v(253) v(226)=((-1d0)+v(201))*v(202)-v(204)/v(201) v(235)=v(172)*v(204)+v(226)*v(257) v(234)=v(173)*v(204)+v(226)*v(256) v(233)=v(174)*v(204)+v(226)*v(255) v(232)=v(176)*v(204)+v(226)*v(254) v(231)=v(177)*v(204)+v(226)*v(253) v(230)=v(178)*v(204)+v(226)*v(252) v(229)=v(192)*v(204)+v(226)*v(251) v(228)=v(193)*v(204)+v(226)*v(250) v(227)=v(194)*v(204)+v(226)*v(249) IF(int(idata(ID_SkipTangent)).eq.1) THEN v(5032)=v(135)*v(229)+v(137)*v(231)+v(136)*v(233) v(5033)=v(136)*v(228)+v(137)*v(230)+v(135)*v(235) v(5034)=v(137)*v(227)+v(136)*v(232)+v(135)*v(234) v(5035)=v(139)*v(229)+v(143)*v(231)+v(141)*v(233) v(5036)=v(141)*v(228)+v(143)*v(230)+v(139)*v(235) v(5037)=v(143)*v(227)+v(141)*v(232)+v(139)*v(234) v(5038)=v(144)*v(229)+v(146)*v(231)+v(145)*v(233) v(5039)=v(145)*v(228)+v(146)*v(230)+v(144)*v(235) v(5040)=v(146)*v(227)+v(145)*v(232)+v(144)*v(234) v(5041)=v(148)*v(229)+v(152)*v(231)+v(150)*v(233) v(5042)=v(150)*v(228)+v(152)*v(230)+v(148)*v(235) v(5043)=v(152)*v(227)+v(150)*v(232)+v(148)*v(234) v(5044)=v(153)*v(229)+v(155)*v(231)+v(154)*v(233) v(5045)=v(154)*v(228)+v(155)*v(230)+v(153)*v(235) v(5046)=v(155)*v(227)+v(154)*v(232)+v(153)*v(234) v(5047)=v(157)*v(229)+v(161)*v(231)+v(159)*v(233) v(5048)=v(159)*v(228)+v(161)*v(230)+v(157)*v(235) v(5049)=v(161)*v(227)+v(159)*v(232)+v(157)*v(234) v(5050)=v(162)*v(229)+v(164)*v(231)+v(163)*v(233) v(5051)=v(163)*v(228)+v(164)*v(230)+v(162)*v(235) v(5052)=v(164)*v(227)+v(163)*v(232)+v(162)*v(234) v(5053)=v(166)*v(229)+v(170)*v(231)+v(168)*v(233) v(5054)=v(168)*v(228)+v(170)*v(230)+v(166)*v(235) v(5055)=v(170)*v(227)+v(168)*v(232)+v(166)*v(234) DO i207=1,24 p(i207)=p(i207)+v(275)*v(5031+i207) ENDDO ELSE v(5256)=0d0 v(5257)=v(135) v(5258)=0d0 v(5259)=0d0 v(5260)=v(139) v(5261)=0d0 v(5262)=0d0 v(5263)=v(144) v(5264)=0d0 v(5265)=0d0 v(5266)=v(148) v(5267)=0d0 v(5268)=0d0 v(5269)=v(153) v(5270)=0d0 v(5271)=0d0 v(5272)=v(157) v(5273)=0d0 v(5274)=0d0 v(5275)=v(162) v(5276)=0d0 v(5277)=0d0 v(5278)=v(166) v(5279)=0d0 v(5232)=0d0 v(5233)=0d0 v(5234)=v(135) v(5235)=0d0 v(5236)=0d0 v(5237)=v(139) v(5238)=0d0 v(5239)=0d0 v(5240)=v(144) v(5241)=0d0 v(5242)=0d0 v(5243)=v(148) v(5244)=0d0 v(5245)=0d0 v(5246)=v(153) v(5247)=0d0 v(5248)=0d0 v(5249)=v(157) v(5250)=0d0 v(5251)=0d0 v(5252)=v(162) v(5253)=0d0 v(5254)=0d0 v(5255)=v(166) v(5208)=v(136) v(5209)=0d0 v(5210)=0d0 v(5211)=v(141) v(5212)=0d0 v(5213)=0d0 v(5214)=v(145) v(5215)=0d0 v(5216)=0d0 v(5217)=v(150) v(5218)=0d0 v(5219)=0d0 v(5220)=v(154) v(5221)=0d0 v(5222)=0d0 v(5223)=v(159) v(5224)=0d0 v(5225)=0d0 v(5226)=v(163) v(5227)=0d0 v(5228)=0d0 v(5229)=v(168) v(5230)=0d0 v(5231)=0d0 v(5184)=0d0 v(5185)=0d0 v(5186)=v(136) v(5187)=0d0 v(5188)=0d0 v(5189)=v(141) v(5190)=0d0 v(5191)=0d0 v(5192)=v(145) v(5193)=0d0 v(5194)=0d0 v(5195)=v(150) v(5196)=0d0 v(5197)=0d0 v(5198)=v(154) v(5199)=0d0 v(5200)=0d0 v(5201)=v(159) v(5202)=0d0 v(5203)=0d0 v(5204)=v(163) v(5205)=0d0 v(5206)=0d0 v(5207)=v(168) v(5160)=v(137) v(5161)=0d0 v(5162)=0d0 v(5163)=v(143) v(5164)=0d0 v(5165)=0d0 v(5166)=v(146) v(5167)=0d0 v(5168)=0d0 v(5169)=v(152) v(5170)=0d0 v(5171)=0d0 v(5172)=v(155) v(5173)=0d0 v(5174)=0d0 v(5175)=v(161) v(5176)=0d0 v(5177)=0d0 v(5178)=v(164) v(5179)=0d0 v(5180)=0d0 v(5181)=v(170) v(5182)=0d0 v(5183)=0d0 v(5136)=0d0 v(5137)=v(137) v(5138)=0d0 v(5139)=0d0 v(5140)=v(143) v(5141)=0d0 v(5142)=0d0 v(5143)=v(146) v(5144)=0d0 v(5145)=0d0 v(5146)=v(152) v(5147)=0d0 v(5148)=0d0 v(5149)=v(155) v(5150)=0d0 v(5151)=0d0 v(5152)=v(161) v(5153)=0d0 v(5154)=0d0 v(5155)=v(164) v(5156)=0d0 v(5157)=0d0 v(5158)=v(170) v(5159)=0d0 v(5112)=v(135) v(5113)=0d0 v(5114)=0d0 v(5115)=v(139) v(5116)=0d0 v(5117)=0d0 v(5118)=v(144) v(5119)=0d0 v(5120)=0d0 v(5121)=v(148) v(5122)=0d0 v(5123)=0d0 v(5124)=v(153) v(5125)=0d0 v(5126)=0d0 v(5127)=v(157) v(5128)=0d0 v(5129)=0d0 v(5130)=v(162) v(5131)=0d0 v(5132)=0d0 v(5133)=v(166) v(5134)=0d0 v(5135)=0d0 v(5088)=0d0 v(5089)=v(136) v(5090)=0d0 v(5091)=0d0 v(5092)=v(141) v(5093)=0d0 v(5094)=0d0 v(5095)=v(145) v(5096)=0d0 v(5097)=0d0 v(5098)=v(150) v(5099)=0d0 v(5100)=0d0 v(5101)=v(154) v(5102)=0d0 v(5103)=0d0 v(5104)=v(159) v(5105)=0d0 v(5106)=0d0 v(5107)=v(163) v(5108)=0d0 v(5109)=0d0 v(5110)=v(168) v(5111)=0d0 v(5064)=0d0 v(5065)=0d0 v(5066)=v(137) v(5067)=0d0 v(5068)=0d0 v(5069)=v(143) v(5070)=0d0 v(5071)=0d0 v(5072)=v(146) v(5073)=0d0 v(5074)=0d0 v(5075)=v(152) v(5076)=0d0 v(5077)=0d0 v(5078)=v(155) v(5079)=0d0 v(5080)=0d0 v(5081)=v(161) v(5082)=0d0 v(5083)=0d0 v(5084)=v(164) v(5085)=0d0 v(5086)=0d0 v(5087)=v(170) DO i223=1,24 v(240)=v(275)*v(5063+i223) v(241)=v(275)*v(5087+i223) v(242)=v(275)*v(5111+i223) v(243)=v(275)*v(5135+i223) v(244)=v(275)*v(5159+i223) v(245)=v(275)*v(5183+i223) v(246)=v(275)*v(5207+i223) v(247)=v(275)*v(5231+i223) v(248)=v(275)*v(5255+i223) v(258)=(v(202)+v(204)/v(201)**2)*(v(240)*v(249)+v(241)*v(250 & )+v(242)*v(251)+v(243)*v(252)+v(244)*v(253)+v(245)*v(254)+v & (246)*v(255)+v(247)*v(256)+v(248)*v(257)) v(259)=v(204)*v(240)+v(226)*(v(192)*v(241)+v(193)*v(242)-v & (172)*v(246)-v(174)*v(248))+v(249)*v(258) v(260)=v(204)*v(241)+v(226)*(v(192)*v(240)+v(194)*v(242)-v & (173)*v(244)-v(177)*v(247))+v(250)*v(258) v(261)=v(204)*v(242)+v(226)*(v(193)*v(240)+v(194)*v(241)-v & (176)*v(243)-v(178)*v(245))+v(251)*v(258) v(262)=v(204)*v(243)+v(226)*(-(v(176)*v(242))-v(192)*v(245) & +v(173)*v(246)+v(174)*v(247))+v(252)*v(258) v(263)=v(204)*v(244)+v(226)*(-(v(173)*v(241))+v(172)*v(245) & -v(193)*v(247)+v(176)*v(248))+v(253)*v(258) v(264)=v(204)*v(245)+v(226)*(-(v(178)*v(242))-v(192)*v(243) & +v(172)*v(244)+v(177)*v(248))+v(254)*v(258) v(265)=v(204)*v(246)+v(226)*(-(v(172)*v(240))+v(173)*v(243) & +v(178)*v(247)-v(194)*v(248))+v(255)*v(258) v(266)=v(226)*(-(v(177)*v(241))+v(174)*v(243)-v(193)*v(244) & +v(178)*v(246))+v(204)*v(247)+v(256)*v(258) v(267)=v(226)*(-(v(174)*v(240))+v(176)*v(244)+v(177)*v(245) & -v(194)*v(246))+v(204)*v(248)+v(257)*v(258) v(5280)=v(135)*v(261)+v(137)*v(263)+v(136)*v(265) v(5281)=v(136)*v(260)+v(137)*v(262)+v(135)*v(267) v(5282)=v(137)*v(259)+v(136)*v(264)+v(135)*v(266) v(5283)=v(139)*v(261)+v(143)*v(263)+v(141)*v(265) v(5284)=v(141)*v(260)+v(143)*v(262)+v(139)*v(267) v(5285)=v(143)*v(259)+v(141)*v(264)+v(139)*v(266) v(5286)=v(144)*v(261)+v(146)*v(263)+v(145)*v(265) v(5287)=v(145)*v(260)+v(146)*v(262)+v(144)*v(267) v(5288)=v(146)*v(259)+v(145)*v(264)+v(144)*v(266) v(5289)=v(148)*v(261)+v(152)*v(263)+v(150)*v(265) v(5290)=v(150)*v(260)+v(152)*v(262)+v(148)*v(267) v(5291)=v(152)*v(259)+v(150)*v(264)+v(148)*v(266) v(5292)=v(153)*v(261)+v(155)*v(263)+v(154)*v(265) v(5293)=v(154)*v(260)+v(155)*v(262)+v(153)*v(267) v(5294)=v(155)*v(259)+v(154)*v(264)+v(153)*v(266) v(5295)=v(157)*v(261)+v(161)*v(263)+v(159)*v(265) v(5296)=v(159)*v(260)+v(161)*v(262)+v(157)*v(267) v(5297)=v(161)*v(259)+v(159)*v(264)+v(157)*v(266) v(5298)=v(162)*v(261)+v(164)*v(263)+v(163)*v(265) v(5299)=v(163)*v(260)+v(164)*v(262)+v(162)*v(267) v(5300)=v(164)*v(259)+v(163)*v(264)+v(162)*v(266) v(5301)=v(166)*v(261)+v(170)*v(263)+v(168)*v(265) v(5302)=v(168)*v(260)+v(170)*v(262)+v(166)*v(267) v(5303)=v(170)*v(259)+v(168)*v(264)+v(166)*v(266) DO i238=i223,24 s(i223,i238)=s(i223,i238)+v(5279+i238) ENDDO ENDDO ENDIF ENDDO END