1 SUBROUTINE DEF_JOINT_AD_ASS_RES(X1, g1, XP1, gP1, R1_0, W1_0, 2 $ X2, g2, XP2, gP2, R2_0, W2_0, S1, D1, o1, o2, F1, M1, F2, M2) 3 IMPLICIT NONE 4 DOUBLE PRECISION X1(3), g1(3), XP1(3), gP1(3), R1_0(3, 3), W1_0(3) 5 DOUBLE PRECISION X2(3), g2(3), XP2(3), gP2(3), R2_0(3, 3), W2_0(3) 6 DOUBLE PRECISION S1(3, 3), D1(3, 3), o1(3), o2(3) 7 DOUBLE PRECISION R1(3, 3), W1(3) 8 DOUBLE PRECISION R2(3, 3), W2(3) 9 DOUBLE PRECISION R1o1(3), R2o2(3) 10 DOUBLE PRECISION dXR1(3), dXPR1(3), dX(3), dXP(3) 11 DOUBLE PRECISION W1xR1o1(3), W2xR2o2(3) 12 DOUBLE PRECISION F1R1(3), F1(3), M1(3), F2(3), M2(3) 13 INTEGER I, J 14 EXTERNAL UPDATE_ROT 15 16 CALL UPDATE_ROT(g1, gP1, R1_0, W1_0, R1, W1) 17 CALL UPDATE_ROT(g2, gP2, R2_0, W2_0, R2, W2) 18 19 DO I=1, 3 20 R1o1(I) = 0D0 21 R2o2(I) = 0D0 22 DO J=1, 3 23 R1o1(I) = R1o1(I) + R1(I, J) * o1(J) 24 R2o2(I) = R2o2(I) + R2(I, J) * o2(J) 25 END DO 26 END DO 27 28 W1xR1o1(1) = W1(2)*R1o1(3)-R1o1(2)*W1(3) 29 W1xR1o1(2) = R1o1(1)*W1(3)-W1(1)*R1o1(3) 30 W1xR1o1(3) = W1(1)*R1o1(2)-R1o1(1)*W1(2) 31 32 W2xR2o2(1) = W2(2)*R2o2(3)-R2o2(2)*W2(3) 33 W2xR2o2(2) = R2o2(1)*W2(3)-W2(1)*R2o2(3) 34 W2xR2o2(3) = W2(1)*R2o2(2)-R2o2(1)*W2(2) 35 36 DO I=1, 3 37 dXR1(I) = X1(I) + R1o1(I) - X2(I) - R2o2(I) 38 dXPR1(I) = XP1(I) + W1xR1o1(I) - XP2(I) - W2xR2o2(I) 39 END DO 40 41 DO I=1, 3 42 dX(I) = 0D0 43 dXP(I) = 0D0 44 DO J=1, 3 45 dX(I) = dX(I) + R1(J, I) * dXR1(J) 46 dXP(I) = dXP(I) + R1(J, I) * dXPR1(J) 47 END DO 48 END DO 49 50 DO I=1, 3 51 F1R1(I) = 0D0 52 DO J=1, 3 53 F1R1(I) = F1R1(I) + S1(I, J) * dX(J) + D1(I, J) * dXP(J) 54 END DO 55 END DO 56 57 DO I=1, 3 58 F1(I) = 0D0 59 DO J=1, 3 60 F1(I) = F1(I) - R1(I, J) * F1R1(J) 61 END DO 62 F2(I) = -F1(I) 63 END DO 64 65 M1(1) = R1o1(2)*F1(3)-F1(2)*R1o1(3) 66 M1(2) = F1(1)*R1o1(3)-R1o1(1)*F1(3) 67 M1(3) = R1o1(1)*F1(2)-F1(1)*R1o1(2) 68 69 M2(1) = R2o2(2)*F2(3)-F2(2)*R2o2(3) 70 M2(2) = F2(1)*R2o2(3)-R2o2(1)*F2(3) 71 M2(3) = R2o2(1)*F2(2)-F2(1)*R2o2(2) 72 END 73 74 SUBROUTINE UPDATE_ROT(g, gP, R_0, W_0, R, W) 75 IMPLICIT NONE 76 DOUBLE PRECISION g(3), gP(3), R_0(3, 3), W_0(3) 77 DOUBLE PRECISION R(3, 3), W(3) 78 DOUBLE PRECISION RDelta(3, 3), Gr(3, 3), d 79 DOUBLE PRECISION tmp1, tmp2, tmp3, tmp4, tmp5, tmp6 80 DOUBLE PRECISION tmp7, tmp8, tmp9 81 INTEGER I, J, K 82 83 tmp1 = g(3)**2 84 tmp2 = g(2)**2 85 tmp3 = g(1)**2 86 tmp4 = g(1)*g(2)*0.5D0 87 tmp5 = g(2)*g(3)*0.5D0 88 tmp6 = g(1)*g(3)*0.5D0 89 90 d=4D0/(tmp1+tmp2+tmp3+4D0) 91 92 RDelta(1,1) = (-tmp1-tmp2)*d*0.5D0+1D0 93 RDelta(1,2) = (tmp4-g(3))*d 94 RDelta(1,3) = (tmp6+g(2))*d 95 RDelta(2,1) = (g(3)+tmp4)*d 96 RDelta(2,2) = (-tmp1-tmp3)*d*0.5D0+1D0 97 RDelta(2,3) = (tmp5-g(1))*d 98 RDelta(3,1) = (tmp6-g(2))*d 99 RDelta(3,2) = (tmp5+g(1))*d 100 RDelta(3,3) = (-tmp2-tmp3)*d*0.5D0+1D0 101 102 tmp7 = 0.5D0*g(3)*d 103 tmp8 = 0.5D0*g(2)*d 104 tmp9 = 0.5D0*g(1)*d 105 106 Gr(1,1) = d 107 Gr(1,2) = -tmp7 108 Gr(1,3) = tmp8 109 Gr(2,1) = tmp7 110 Gr(2,2) = d 111 Gr(2,3) = -tmp9 112 Gr(3,1) = -tmp8 113 Gr(3,2) = tmp9 114 Gr(3,3) = d 115 116 DO I=1, 3 117 DO J=1, 3 118 R(I, J)=0D0 119 DO K=1, 3 120 R(I, J) = R(I, J) + RDelta(I, K) * R_0(K, J) 121 END DO 122 END DO 123 END DO 124 125 DO I=1, 3 126 W(I) = 0D0 127 DO J=1, 3 128 W(I) = W(I) + Gr(I, J) * gP(J) + RDelta(I, J) * W_0(J) 129 END DO 130 END DO 131 132 END 133 134 135C Generated by TAPENADE (INRIA, Tropics team) 136C Tapenade 3.5 (r3618) - 22 Dec 2010 11:39 137C 138C Differentiation of def_joint_ad_ass_res in forward (tangent) mode: (multi-directional mode) 139C variations of useful results: f1 f2 m1 m2 140C with respect to varying inputs: gp1 gp2 xp1 xp2 g1 g2 x1 x2 141C RW status of diff variables: gp1:in gp2:in f1:out f2:out m1:out 142C m2:out xp1:in xp2:in g1:in g2:in x1:in x2:in 143 SUBROUTINE DEF_JOINT_AD_ASS_RES_DV(x1, x1d, g1, g1d, xp1, xp1d, 144 + gp1, gp1d, r1_0, w1_0, x2, x2d 145 + , g2, g2d, xp2, xp2d, gp2, gp2d 146 + , r2_0, w2_0, s1, d1, o1, o2, 147 + f1, f1d, m1, m1d, f2, f2d, m2, 148 + m2d, nbdirs) 149 IMPLICIT INTEGER (N-N) 150 PARAMETER (nbdirsmax=12) 151C Hint: nbdirsmax should be the maximum number of differentiation directions 152 DOUBLE PRECISION x1(3), g1(3), xp1(3), gp1(3), r1_0(3, 3), w1_0(3) 153 DOUBLE PRECISION x1d(nbdirsmax, 3), g1d(nbdirsmax, 3), xp1d( 154 + nbdirsmax, 3), gp1d(nbdirsmax, 3) 155 DOUBLE PRECISION x2(3), g2(3), xp2(3), gp2(3), r2_0(3, 3), w2_0(3) 156 DOUBLE PRECISION x2d(nbdirsmax, 3), g2d(nbdirsmax, 3), xp2d( 157 + nbdirsmax, 3), gp2d(nbdirsmax, 3) 158 DOUBLE PRECISION s1(3, 3), d1(3, 3), o1(3), o2(3) 159 DOUBLE PRECISION r1(3, 3), w1(3) 160 DOUBLE PRECISION r1d(nbdirsmax, 3, 3), w1d(nbdirsmax, 3) 161 DOUBLE PRECISION r2(3, 3), w2(3) 162 DOUBLE PRECISION r2d(nbdirsmax, 3, 3), w2d(nbdirsmax, 3) 163 DOUBLE PRECISION r1o1(3), r2o2(3) 164 DOUBLE PRECISION r1o1d(nbdirsmax, 3), r2o2d(nbdirsmax, 3) 165 DOUBLE PRECISION dxr1(3), dxpr1(3), dx(3), dxp(3) 166 DOUBLE PRECISION dxr1d(nbdirsmax, 3), dxpr1d(nbdirsmax, 3), dxd( 167 + nbdirsmax, 3), dxpd(nbdirsmax, 3) 168 DOUBLE PRECISION w1xr1o1(3), w2xr2o2(3) 169 DOUBLE PRECISION w1xr1o1d(nbdirsmax, 3), w2xr2o2d(nbdirsmax, 3) 170 DOUBLE PRECISION f1r1(3), f1(3), m1(3), f2(3), m2(3) 171 DOUBLE PRECISION f1r1d(nbdirsmax, 3), f1d(nbdirsmax, 3), m1d( 172 + nbdirsmax, 3), f2d(nbdirsmax, 3), m2d(nbdirsmax, 173 + 3) 174 INTEGER i, j 175 EXTERNAL UPDATE_ROT 176 EXTERNAL UPDATE_ROT_DV 177 INTEGER nd 178 INTEGER nbdirs 179 INTEGER ii1 180 CALL UPDATE_ROT_DV(g1, g1d, gp1, gp1d, r1_0, w1_0, r1, r1d, w1, 181 + w1d, nbdirs) 182 CALL UPDATE_ROT_DV(g2, g2d, gp2, gp2d, r2_0, w2_0, r2, r2d, w2, 183 + w2d, nbdirs) 184 DO nd=1,nbdirs 185 DO ii1=1,3 186 r2o2d(nd, ii1) = 0.D0 187 ENDDO 188 ENDDO 189 DO nd=1,nbdirs 190 DO ii1=1,3 191 r1o1d(nd, ii1) = 0.D0 192 ENDDO 193 ENDDO 194 DO i=1,3 195 DO nd=1,nbdirs 196 r1o1d(nd, i) = 0.D0 197 r2o2d(nd, i) = 0.D0 198 ENDDO 199 r1o1(i) = 0d0 200 r2o2(i) = 0d0 201 DO j=1,3 202 DO nd=1,nbdirs 203 r1o1d(nd, i) = r1o1d(nd, i) + o1(j)*r1d(nd, i, j) 204 r2o2d(nd, i) = r2o2d(nd, i) + o2(j)*r2d(nd, i, j) 205 ENDDO 206 r1o1(i) = r1o1(i) + r1(i, j)*o1(j) 207 r2o2(i) = r2o2(i) + r2(i, j)*o2(j) 208 ENDDO 209 ENDDO 210 DO nd=1,nbdirs 211 DO ii1=1,3 212 w1xr1o1d(nd, ii1) = 0.D0 213 ENDDO 214 w1xr1o1d(nd, 1) = w1d(nd, 2)*r1o1(3) + w1(2)*r1o1d(nd, 3) - 215 + r1o1d(nd, 2)*w1(3) - r1o1(2)*w1d(nd, 3) 216 w1xr1o1d(nd, 2) = r1o1d(nd, 1)*w1(3) + r1o1(1)*w1d(nd, 3) - w1d( 217 + nd, 1)*r1o1(3) - w1(1)*r1o1d(nd, 3) 218 w1xr1o1d(nd, 3) = w1d(nd, 1)*r1o1(2) + w1(1)*r1o1d(nd, 2) - 219 + r1o1d(nd, 1)*w1(2) - r1o1(1)*w1d(nd, 2) 220 DO ii1=1,3 221 w2xr2o2d(nd, ii1) = 0.D0 222 ENDDO 223 w2xr2o2d(nd, 1) = w2d(nd, 2)*r2o2(3) + w2(2)*r2o2d(nd, 3) - 224 + r2o2d(nd, 2)*w2(3) - r2o2(2)*w2d(nd, 3) 225 w2xr2o2d(nd, 2) = r2o2d(nd, 1)*w2(3) + r2o2(1)*w2d(nd, 3) - w2d( 226 + nd, 1)*r2o2(3) - w2(1)*r2o2d(nd, 3) 227 w2xr2o2d(nd, 3) = w2d(nd, 1)*r2o2(2) + w2(1)*r2o2d(nd, 2) - 228 + r2o2d(nd, 1)*w2(2) - r2o2(1)*w2d(nd, 2) 229 ENDDO 230 w1xr1o1(1) = w1(2)*r1o1(3) - r1o1(2)*w1(3) 231 w1xr1o1(2) = r1o1(1)*w1(3) - w1(1)*r1o1(3) 232 w1xr1o1(3) = w1(1)*r1o1(2) - r1o1(1)*w1(2) 233 w2xr2o2(1) = w2(2)*r2o2(3) - r2o2(2)*w2(3) 234 w2xr2o2(2) = r2o2(1)*w2(3) - w2(1)*r2o2(3) 235 w2xr2o2(3) = w2(1)*r2o2(2) - r2o2(1)*w2(2) 236 DO nd=1,nbdirs 237 DO ii1=1,3 238 dxr1d(nd, ii1) = 0.D0 239 ENDDO 240 ENDDO 241 DO nd=1,nbdirs 242 DO ii1=1,3 243 dxpr1d(nd, ii1) = 0.D0 244 ENDDO 245 ENDDO 246 DO i=1,3 247 DO nd=1,nbdirs 248 dxr1d(nd, i) = x1d(nd, i) + r1o1d(nd, i) - x2d(nd, i) - r2o2d( 249 + nd, i) 250 dxpr1d(nd, i) = xp1d(nd, i) + w1xr1o1d(nd, i) - xp2d(nd, i) - 251 + w2xr2o2d(nd, i) 252 ENDDO 253 dxr1(i) = x1(i) + r1o1(i) - x2(i) - r2o2(i) 254 dxpr1(i) = xp1(i) + w1xr1o1(i) - xp2(i) - w2xr2o2(i) 255 ENDDO 256 DO nd=1,nbdirs 257 DO ii1=1,3 258 dxd(nd, ii1) = 0.D0 259 ENDDO 260 ENDDO 261 DO nd=1,nbdirs 262 DO ii1=1,3 263 dxpd(nd, ii1) = 0.D0 264 ENDDO 265 ENDDO 266 DO i=1,3 267 DO nd=1,nbdirs 268 dxd(nd, i) = 0.D0 269 dxpd(nd, i) = 0.D0 270 ENDDO 271 dx(i) = 0d0 272 dxp(i) = 0d0 273 DO j=1,3 274 DO nd=1,nbdirs 275 dxd(nd, i) = dxd(nd, i) + r1d(nd, j, i)*dxr1(j) + r1(j, i)* 276 + dxr1d(nd, j) 277 dxpd(nd, i) = dxpd(nd, i) + r1d(nd, j, i)*dxpr1(j) + r1(j, i 278 + )*dxpr1d(nd, j) 279 ENDDO 280 dx(i) = dx(i) + r1(j, i)*dxr1(j) 281 dxp(i) = dxp(i) + r1(j, i)*dxpr1(j) 282 ENDDO 283 ENDDO 284 DO nd=1,nbdirs 285 DO ii1=1,3 286 f1r1d(nd, ii1) = 0.D0 287 ENDDO 288 ENDDO 289 DO i=1,3 290 DO nd=1,nbdirs 291 f1r1d(nd, i) = 0.D0 292 ENDDO 293 f1r1(i) = 0d0 294 DO j=1,3 295 DO nd=1,nbdirs 296 f1r1d(nd, i) = f1r1d(nd, i) + s1(i, j)*dxd(nd, j) + d1(i, j) 297 + *dxpd(nd, j) 298 ENDDO 299 f1r1(i) = f1r1(i) + s1(i, j)*dx(j) + d1(i, j)*dxp(j) 300 ENDDO 301 ENDDO 302 DO nd=1,nbdirs 303 DO ii1=1,3 304 f1d(nd, ii1) = 0.D0 305 ENDDO 306 ENDDO 307 DO nd=1,nbdirs 308 DO ii1=1,3 309 f2d(nd, ii1) = 0.D0 310 ENDDO 311 ENDDO 312 DO i=1,3 313 DO nd=1,nbdirs 314 f1d(nd, i) = 0.D0 315 ENDDO 316 f1(i) = 0d0 317 DO j=1,3 318 DO nd=1,nbdirs 319 f1d(nd, i) = f1d(nd, i) - r1d(nd, i, j)*f1r1(j) - r1(i, j)* 320 + f1r1d(nd, j) 321 ENDDO 322 f1(i) = f1(i) - r1(i, j)*f1r1(j) 323 ENDDO 324 DO nd=1,nbdirs 325 f2d(nd, i) = -f1d(nd, i) 326 ENDDO 327 f2(i) = -f1(i) 328 ENDDO 329 DO nd=1,nbdirs 330 DO ii1=1,3 331 m1d(nd, ii1) = 0.D0 332 ENDDO 333 m1d(nd, 1) = r1o1d(nd, 2)*f1(3) + r1o1(2)*f1d(nd, 3) - f1d(nd, 2 334 + )*r1o1(3) - f1(2)*r1o1d(nd, 3) 335 m1d(nd, 2) = f1d(nd, 1)*r1o1(3) + f1(1)*r1o1d(nd, 3) - r1o1d(nd 336 + , 1)*f1(3) - r1o1(1)*f1d(nd, 3) 337 m1d(nd, 3) = r1o1d(nd, 1)*f1(2) + r1o1(1)*f1d(nd, 2) - f1d(nd, 1 338 + )*r1o1(2) - f1(1)*r1o1d(nd, 2) 339 DO ii1=1,3 340 m2d(nd, ii1) = 0.D0 341 ENDDO 342 m2d(nd, 1) = r2o2d(nd, 2)*f2(3) + r2o2(2)*f2d(nd, 3) - f2d(nd, 2 343 + )*r2o2(3) - f2(2)*r2o2d(nd, 3) 344 m2d(nd, 2) = f2d(nd, 1)*r2o2(3) + f2(1)*r2o2d(nd, 3) - r2o2d(nd 345 + , 1)*f2(3) - r2o2(1)*f2d(nd, 3) 346 m2d(nd, 3) = r2o2d(nd, 1)*f2(2) + r2o2(1)*f2d(nd, 2) - f2d(nd, 1 347 + )*r2o2(2) - f2(1)*r2o2d(nd, 2) 348 ENDDO 349 m1(1) = r1o1(2)*f1(3) - f1(2)*r1o1(3) 350 m1(2) = f1(1)*r1o1(3) - r1o1(1)*f1(3) 351 m1(3) = r1o1(1)*f1(2) - f1(1)*r1o1(2) 352 m2(1) = r2o2(2)*f2(3) - f2(2)*r2o2(3) 353 m2(2) = f2(1)*r2o2(3) - r2o2(1)*f2(3) 354 m2(3) = r2o2(1)*f2(2) - f2(1)*r2o2(2) 355 END 356 357C Generated by TAPENADE (INRIA, Tropics team) 358C Tapenade 3.5 (r3618) - 22 Dec 2010 11:39 359C 360C Differentiation of update_rot in forward (tangent) mode: (multi-directional mode) 361C variations of useful results: r w 362C with respect to varying inputs: g gp 363 SUBROUTINE UPDATE_ROT_DV(g, gd, gp, gpd, r_0, w_0, r, rd, w, wd, 364 + nbdirs) 365 IMPLICIT INTEGER (N-N) 366 PARAMETER (nbdirsmax=12) 367C Hint: nbdirsmax should be the maximum number of differentiation directions 368 DOUBLE PRECISION g(3), gp(3), r_0(3, 3), w_0(3) 369 DOUBLE PRECISION gd(nbdirsmax, 3), gpd(nbdirsmax, 3) 370 DOUBLE PRECISION r(3, 3), w(3) 371 DOUBLE PRECISION rd(nbdirsmax, 3, 3), wd(nbdirsmax, 3) 372 DOUBLE PRECISION rdelta(3, 3), gr(3, 3), d 373 DOUBLE PRECISION rdeltad(nbdirsmax, 3, 3), grd(nbdirsmax, 3, 3), 374 + dd(nbdirsmax) 375 DOUBLE PRECISION tmp1, tmp2, tmp3, tmp4, tmp5, tmp6 376 DOUBLE PRECISION tmp1d(nbdirsmax), tmp2d(nbdirsmax), tmp3d( 377 + nbdirsmax), tmp4d(nbdirsmax), tmp5d(nbdirsmax), 378 + tmp6d(nbdirsmax) 379 DOUBLE PRECISION tmp7, tmp8, tmp9 380 DOUBLE PRECISION tmp7d(nbdirsmax), tmp8d(nbdirsmax), tmp9d( 381 + nbdirsmax) 382 INTEGER i, j, k 383 INTEGER nd 384 INTEGER nbdirs 385 INTEGER ii2 386 INTEGER ii1 387 tmp1 = g(3)**2 388 tmp2 = g(2)**2 389 tmp3 = g(1)**2 390 tmp4 = g(1)*g(2)*0.5d0 391 tmp5 = g(2)*g(3)*0.5d0 392 tmp6 = g(1)*g(3)*0.5d0 393 d = 4d0/(tmp1+tmp2+tmp3+4d0) 394 DO nd=1,nbdirs 395 tmp1d(nd) = 2*g(3)*gd(nd, 3) 396 tmp2d(nd) = 2*g(2)*gd(nd, 2) 397 tmp3d(nd) = 2*g(1)*gd(nd, 1) 398 tmp4d(nd) = 0.5d0*(gd(nd, 1)*g(2)+g(1)*gd(nd, 2)) 399 tmp5d(nd) = 0.5d0*(gd(nd, 2)*g(3)+g(2)*gd(nd, 3)) 400 tmp6d(nd) = 0.5d0*(gd(nd, 1)*g(3)+g(1)*gd(nd, 3)) 401 dd(nd) = -(4d0*(tmp1d(nd)+tmp2d(nd)+tmp3d(nd))/(tmp1+tmp2+tmp3+ 402 + 4d0)**2) 403 DO ii1=1,3 404 DO ii2=1,3 405 rdeltad(nd, ii2, ii1) = 0.D0 406 ENDDO 407 ENDDO 408 rdeltad(nd, 1, 1) = 0.5d0*((-tmp1d(nd)-tmp2d(nd))*d+(-tmp1-tmp2) 409 + *dd(nd)) 410 rdeltad(nd, 1, 2) = (tmp4d(nd)-gd(nd, 3))*d + (tmp4-g(3))*dd(nd) 411 rdeltad(nd, 1, 3) = (tmp6d(nd)+gd(nd, 2))*d + (tmp6+g(2))*dd(nd) 412 rdeltad(nd, 2, 1) = (gd(nd, 3)+tmp4d(nd))*d + (g(3)+tmp4)*dd(nd) 413 rdeltad(nd, 2, 2) = 0.5d0*((-tmp1d(nd)-tmp3d(nd))*d+(-tmp1-tmp3) 414 + *dd(nd)) 415 rdeltad(nd, 2, 3) = (tmp5d(nd)-gd(nd, 1))*d + (tmp5-g(1))*dd(nd) 416 rdeltad(nd, 3, 1) = (tmp6d(nd)-gd(nd, 2))*d + (tmp6-g(2))*dd(nd) 417 rdeltad(nd, 3, 2) = (tmp5d(nd)+gd(nd, 1))*d + (tmp5+g(1))*dd(nd) 418 rdeltad(nd, 3, 3) = 0.5d0*((-tmp2d(nd)-tmp3d(nd))*d+(-tmp2-tmp3) 419 + *dd(nd)) 420 tmp7d(nd) = 0.5d0*(gd(nd, 3)*d+g(3)*dd(nd)) 421 tmp8d(nd) = 0.5d0*(gd(nd, 2)*d+g(2)*dd(nd)) 422 tmp9d(nd) = 0.5d0*(gd(nd, 1)*d+g(1)*dd(nd)) 423 DO ii1=1,3 424 DO ii2=1,3 425 grd(nd, ii2, ii1) = 0.D0 426 ENDDO 427 ENDDO 428 grd(nd, 1, 1) = dd(nd) 429 grd(nd, 1, 2) = -tmp7d(nd) 430 grd(nd, 1, 3) = tmp8d(nd) 431 grd(nd, 2, 1) = tmp7d(nd) 432 grd(nd, 2, 2) = dd(nd) 433 grd(nd, 2, 3) = -tmp9d(nd) 434 grd(nd, 3, 1) = -tmp8d(nd) 435 grd(nd, 3, 2) = tmp9d(nd) 436 grd(nd, 3, 3) = dd(nd) 437 ENDDO 438 rdelta(1, 1) = (-tmp1-tmp2)*d*0.5d0 + 1d0 439 rdelta(1, 2) = (tmp4-g(3))*d 440 rdelta(1, 3) = (tmp6+g(2))*d 441 rdelta(2, 1) = (g(3)+tmp4)*d 442 rdelta(2, 2) = (-tmp1-tmp3)*d*0.5d0 + 1d0 443 rdelta(2, 3) = (tmp5-g(1))*d 444 rdelta(3, 1) = (tmp6-g(2))*d 445 rdelta(3, 2) = (tmp5+g(1))*d 446 rdelta(3, 3) = (-tmp2-tmp3)*d*0.5d0 + 1d0 447 tmp7 = 0.5d0*g(3)*d 448 tmp8 = 0.5d0*g(2)*d 449 tmp9 = 0.5d0*g(1)*d 450 gr(1, 1) = d 451 gr(1, 2) = -tmp7 452 gr(1, 3) = tmp8 453 gr(2, 1) = tmp7 454 gr(2, 2) = d 455 gr(2, 3) = -tmp9 456 gr(3, 1) = -tmp8 457 gr(3, 2) = tmp9 458 gr(3, 3) = d 459 DO nd=1,nbdirs 460 DO ii1=1,3 461 DO ii2=1,3 462 rd(nd, ii2, ii1) = 0.D0 463 ENDDO 464 ENDDO 465 ENDDO 466 DO i=1,3 467 DO j=1,3 468 DO nd=1,nbdirs 469 rd(nd, i, j) = 0.D0 470 ENDDO 471 r(i, j) = 0d0 472 DO k=1,3 473 DO nd=1,nbdirs 474 rd(nd, i, j) = rd(nd, i, j) + r_0(k, j)*rdeltad(nd, i, k) 475 ENDDO 476 r(i, j) = r(i, j) + rdelta(i, k)*r_0(k, j) 477 ENDDO 478 ENDDO 479 ENDDO 480 DO nd=1,nbdirs 481 DO ii1=1,3 482 wd(nd, ii1) = 0.D0 483 ENDDO 484 ENDDO 485 DO i=1,3 486 DO nd=1,nbdirs 487 wd(nd, i) = 0.D0 488 ENDDO 489 w(i) = 0d0 490 DO j=1,3 491 DO nd=1,nbdirs 492 wd(nd, i) = wd(nd, i) + grd(nd, i, j)*gp(j) + gr(i, j)*gpd( 493 + nd, j) + w_0(j)*rdeltad(nd, i, j) 494 ENDDO 495 w(i) = w(i) + gr(i, j)*gp(j) + rdelta(i, j)*w_0(j) 496 ENDDO 497 ENDDO 498 END 499