1 subroutine e3source(xx, src) 2c----------------------------------------------------------------------- 3c 4c this routine computes the body force term. 5c 6c currently this computes a swirl body with the axis alligned with 7c the z-coordinate 8c 9c----------------------------------------------------------------------- 10 11 include "common.h" 12 13 real*8 xx(npro,nsd), src(npro,nsd) 14 15 real*8 nu 16 17 real*8 r, Stheta, dpdz, rP5 18 if(datmat(2,5,1).eq. 0) then ! calc swirl 19c 20c This is the body force which will drive a swirl in a pipe flow 21c 22 bigR = 0.5d0 23 dpdz = datmat(1,5,1) 24 do iel = 1, npro 25 26 r = sqrt( xx(iel,1)**2 + xx(iel,2)**2) 27 rP5 = (r/bigR)**5 28 29 Stheta = dpdz * sin(0.5*pi*rP5) 30 31 src(iel,1) = -xx(iel,2)/r * Stheta 32 src(iel,2) = xx(iel,1)/r * Stheta 33 src(iel,3) = dpdz 34 enddo 35 else ! contrived test problem 36 37c$$$ nu=datmat(1,2,1) 38c$$$ omag=datmat(3,5,1) ! frame rotation rate 39c$$$ e1 = one/sqrt(3.0d0) ! axis of rotation 40c$$$ e2 = e1 41c$$$ e3 = e1 42c$$$ om1=omag*e1 43c$$$ om2=omag*e2 44c$$$ om3=omag*e3 45c$$$ w=0 46c$$$ 47c$$$ do iel = 1, npro 48c$$$ 49c$$$ x = xx(iel,1) 50c$$$ y = xx(iel,2) 51c$$$ z = xx(iel,3) 52c$$$ 53c$$$c 54c$$$c The following are MAPLE generated forcing functions for 55c$$$c a contrived flow in a rotating reference frame. 56c$$$c 57c$$$ 58c$$$ t1 = x**2 59c$$$ t2 = t1**2 60c$$$ t5 = dexp(14.D0*x) 61c$$$ t6 = t2*x*t5 62c$$$ t7 = y**2 63c$$$ t8 = t7**2 64c$$$ t9 = t8*t7 65c$$$ t12 = t2*t5 66c$$$ t15 = nu*t1 67c$$$ t17 = dexp(7.D0*x) 68c$$$ t18 = t7*y 69c$$$ t19 = t17*t18 70c$$$ t22 = t1*x 71c$$$ t23 = nu*t22 72c$$$ t24 = t17*y 73c$$$ t29 = t17*t7 74c$$$ t34 = nu*x 75c$$$ t43 = nu*t2 76c$$$ t46 = -40.D0*t6*t9-6.D0*t12*t7+92.D0*t15*t19+132.D0*t23*t24+80.D0* 77c$$$ #t6*t18-252.D0*t23*t29+168.D0*t23*t19+96.D0*t34*t29-64.D0*t34*t19+2 78c$$$ #2.D0*t15*t24-138.D0*t15*t29+294.D0*t43*t29 79c$$$ t50 = t2*t22*t5 80c$$$ t57 = nu*t17 81c$$$ t60 = t8*y 82c$$$ t65 = t2**2 83c$$$ t66 = t65*t5 84c$$$ t73 = t22*t5 85c$$$ t77 = t2*t1*t5 86c$$$ t80 = -196.D0*t43*t19-96.D0*t50*t9-122.D0*t43*t24-32.D0*t34*t24-4. 87c$$$ #D0*t57*y+36.D0*t12*t60+192.D0*t50*t18+98.D0*t66*t8+24.D0*t12*t18+2 88c$$$ #8.D0*t66*t9-24.D0*t73*t60-336.D0*t77*t60 89c$$$ t106 = -336.D0*t50*t8-12.D0*t12*t9+8.D0*t73*t9-140.D0*t6*t8+28.D0* 90c$$$ #t73*t8+12.D0*t15*t17+12.D0*t43*t17+288.D0*t50*t60+112.D0*t77*t9-22 91c$$$ #4.D0*t77*t18-56.D0*t66*t18+120.D0*t6*t60 92c$$$ t131 = -8.D0*t57*t18-20.D0*t6*t7+56.D0*t77*t7-42.D0*t12*t8-24.D0*t 93c$$$ #23*t17-48.D0*t50*t7+4.D0*t73*t7-16.D0*t73*t18+392.D0*t77*t8+14.D0* 94c$$$ #t66*t7-84.D0*t66*t60+12.D0*t57*t7 95c$$$ fx = t46+t80+t106+t131 96c$$$ 97c$$$ 98c$$$ t1 = x**2 99c$$$ t2 = t1**2 100c$$$ t3 = nu*t2 101c$$$ t5 = dexp(7.D0*x) 102c$$$ t6 = t5*y 103c$$$ t9 = y**2 104c$$$ t10 = t9**2 105c$$$ t11 = nu*t10 106c$$$ t18 = t1*x 107c$$$ t22 = nu*x 108c$$$ t25 = nu*t18 109c$$$ t32 = dexp(14.D0*x) 110c$$$ t33 = t2*t32 111c$$$ t39 = t2*t1*t32 112c$$$ t40 = t10*t9 113c$$$ t43 = t1*t32 114c$$$ t46 = t9*y 115c$$$ t47 = t10*t46 116c$$$ t50 = -84.D0*t3*t6+343.D0*t11*t5*t2+66.D0*t11*t5*x-98.D0*t11*t5*t1 117c$$$ #8-24.D0*t22*t6+120.D0*t25*t6-287.D0*t11*t5*t1-140.D0*t33*t10+14.D0 118c$$$ #*t3*t5-56.D0*t39*t40-28.D0*t43*t40+8.D0*t43*t47 119c$$$ t52 = t2*x*t32 120c$$$ t55 = t18*t32 121c$$$ t62 = t10*y 122c$$$ t69 = nu*t5 123c$$$ t80 = 120.D0*t52*t10-32.D0*t55*t47+56.D0*t33*t47+30.D0*t11*t5-216. 124c$$$ #D0*t52*t62-144.D0*t55*t62+72.D0*t39*t62-60.D0*t69*t46+30.D0*t69*t9 125c$$$ #-20.D0*t25*t5-20.D0*t43*t10+36.D0*t43*t62 126c$$$ t86 = nu*t1 127c$$$ t89 = t5*t9 128c$$$ t92 = t5*t46 129c$$$ t109 = 112.D0*t55*t40+80.D0*t55*t10-12.D0*t86*t6-275.D0*t86*t89+57 130c$$$ #4.D0*t86*t92-218.D0*t25*t89+196.D0*t25*t92+90.D0*t22*t89-132.D0*t2 131c$$$ #2*t92+427.D0*t3*t89+8.D0*t39*t46+4.D0*t43*t46 132c$$$ t134 = 16.D0*t39*t47+28.D0*t33*t46-48.D0*t52*t47+252.D0*t33*t62+16 133c$$$ #8.D0*t52*t40-24.D0*t52*t46-686.D0*t3*t92+2.D0*t86*t5-40.D0*t39*t10 134c$$$ #-196.D0*t33*t40-16.D0*t55*t46+4.D0*t22*t5 135c$$$ fy = t50+t80+t109+t134 136c$$$ 137c$$$ 138c$$$ t3 = om3*x 139c$$$ t5 = dexp(7.D0*x) 140c$$$ t7 = x**2 141c$$$ t12 = y**2 142c$$$ t15 = (-1.D0+y)**2 143c$$$ fxa = 2.D0*om2*w+2.D0*t3*t5*(2.D0+x-10.D0*t7+7.D0*t7*x)*t12*t15+om 144c$$$ #2*(om1*y-om2*x)-om3*(t3-om1*z) 145c$$$ 146c$$$ t1 = x**2 147c$$$ t4 = (-1.D0+x)**2 148c$$$ t7 = dexp(7.D0*x) 149c$$$ t10 = y**2 150c$$$ fya = 4.D0*om3*t1*t4*t7*y*(1.D0-3.D0*y+2.D0*t10)-2.D0*om1*w+om3*(o 151c$$$ #m2*z-om3*y)-om1*(om1*y-om2*x) 152c$$$ 153c$$$ 154c$$$ t3 = dexp(7.D0*x) 155c$$$ t5 = x**2 156c$$$ t10 = y**2 157c$$$ t13 = (-1.D0+y)**2 158c$$$ t19 = (-1.D0+x)**2 159c$$$ fza = -2.D0*om1*x*t3*(2.D0+x-10.D0*t5+7.D0*t5*x)*t10*t13-4.D0*om2* 160c$$$ #t5*t19*t3*y*(1.D0-3.D0*y+2.D0*t10)+om1*(om3*x-om1*z)-om2*(om2*z-om 161c$$$ #3*y) 162c$$$ 163c$$$ src(iel,1) = fx + fxa 164c$$$ src(iel,2) = fy + fya 165c$$$ src(iel,3) = fza 166c$$$ enddo 167!Analytic lid case 168 169 do iel = 1, npro 170 x = xx(iel,1) 171 y = xx(iel,2) 172 z = xx(iel,3) 173 Re = 1.0/datmat(1,2,1) 174c 175c The following are MAPLE generated forcing functions for 176c a Lid Driven cavity flow with an analytic solution 177c 178 t2 = x**2 179 t3 = t2**2 180 t4 = t3*x 181 t7 = t2*x 182 t8 = 8*t7 183 t13 = y**2 184 t20 = t13**2 185 t21 = t20-t13 186 t26 = t3**2 187 t30 = t3*t2 188 t37 = t13*y 189 t54 = -8/Re*(24.E0/5.E0*t4-12*t3+t8+2*(4*t7-6*t2+2*x)*(12 190 & *t13-2)+(24*x-12)*t21)-64*(t26/2-2*t3*t7+3*t30-2*t4+t3 191 & /2)*(-24*t20*y+8*t37-4*y)+64*t21*(4*t37-2*y)*(-4*t30+12 192 & *t4-14*t3+t8-2*t2) 193 194 src(iel,1) = 0.0 195 src(iel,2) = t54 196 src(iel,3) = 0.0 197 enddo 198 199 endif 200 return 201 end 202 203 204!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 205!****************************************************************** 206!------------------------------------------------------------------- 207 208 209 210 subroutine e3sourceSclr ( Sclr, Sdot, gradS, 211 & dwl, shape_funct, shg, 212 & yl, dxidx, rmu, 213 & u1, u2, u3, xl, 214 & srcR, srcL, uMod, 215 & srcRat) 216 217 218c----------------------------------------------------------------------- 219c 220c calculate the coefficients for the Spalart-Allmaras 221c turbulence model and its jacobian 222c 223c input: 224c Sclr : The scalar that is being transported in this pde. 225c gradS : spatial gradient of Sclr 226c rmu : kinematic viscosity 227c 228c output: 229c rmu : diffusion for eddy viscosity equation 230c srcR : residual terms for 231c (srcR * turb) 232c srcL : jacobian of srcR 233c 234c----------------------------------------------------------------------- 235 use turbSA 236 use turbKE 237 include "common.h" 238c coming in 239 real*8 Sclr(npro), Sdot(npro), 240 & gradS(npro,nsd), dwl(npro,nenl), 241 & shape_funct(npro,nshl), shg(npro,nshl,nsd), 242 & yl(npro,nenl,ndof), dxidx(npro,nsd,nsd), 243 & rmu(npro), u1(npro), 244 & u2(npro), u3(npro), 245 & xl(npro,nenl,nsd) 246c going out 247 real*8 srcR(npro), srcL(npro), 248 & uMod(npro,nsd) 249c used locally 250 251 real*8 gradv(npro,nsd,nsd), absVort(npro), 252 & dwall(npro) 253 real*8 chi, chiP3, fv1, fv2, st, r, 254 & g, fw, s, viscInv, k2d2Inv, 255 & dP2Inv, sixth, tmp, tmp1, p, dp, 256 & fv1p, fv2p, stp, gp, fwp, rp, 257 & chiP2, mytmp(npro), sign_levelset(npro), 258 & sclr_ls(npro) 259c Kay-Epsilon 260 real*8 Ry(npro), Rt(npro), RtP2(npro), f2(npro), f1(npro), 261 & kay(npro), epsilon(npro), fmu(npro), fmui(npro), 262 & srcjac(npro,4), srcRat(npro) 263 integer e 264c---------------------------------------------------------------------- 265c 266c -- -- -- -- 267c | cb2 | 1 | | 268c phi, + |u - -----phi, | phi, - -----|(nu+phi)phi, | 269c t | i sigma i| i sigma| i| 270c -- -- -- --, 271c i 272c 273c ~ phi^2 274c - cb1*S*phi + cw1*fw*----- 275c d^2 276c 277c---------------------------------------------------------------------- 278c 279 if(iRANS.eq.-1) then ! spalart almaras model 280c 281c.... compute the global gradient of u 282c 283 gradV = zero 284 do n = 1, nshl 285c 286c du_i/dx_j 287c 288c i j indices match array where V is the velocity (u in our notes) 289 gradV(:,1,1) = gradV(:,1,1) + shg(:,n,1) * yl(:,n,2) 290 gradV(:,2,1) = gradV(:,2,1) + shg(:,n,1) * yl(:,n,3) 291 gradV(:,3,1) = gradV(:,3,1) + shg(:,n,1) * yl(:,n,4) 292c 293 gradV(:,1,2) = gradV(:,1,2) + shg(:,n,2) * yl(:,n,2) 294 gradV(:,2,2) = gradV(:,2,2) + shg(:,n,2) * yl(:,n,3) 295 gradV(:,3,2) = gradV(:,3,2) + shg(:,n,2) * yl(:,n,4) 296c 297 gradV(:,1,3) = gradV(:,1,3) + shg(:,n,3) * yl(:,n,2) 298 gradV(:,2,3) = gradV(:,2,3) + shg(:,n,3) * yl(:,n,3) 299 gradV(:,3,3) = gradV(:,3,3) + shg(:,n,3) * yl(:,n,4) 300c a j u a i 301c from our notes where we had N_{a,j} = dN_a/dx_j note that i is off by one because p was first in yl vector 302c 303 enddo 304c 305c.... magnitude of vorticity 306c 307 absVort = sqrt( (gradV(:,2,3) - gradV(:,3,2)) ** 2 308 & + (gradV(:,3,1) - gradV(:,1,3)) ** 2 309 & + (gradV(:,1,2) - gradV(:,2,1)) ** 2 ) 310 dwall = zero 311 do n = 1, nenl 312 dwall = dwall + shape_funct(:,n) * dwl(:,n) 313 enddo 314 315 sixth = 1.0/6.0 316c 317c.... compute source and its jacobian 318c 319 do e = 1, npro 320 SclrNN= max(Sclr(e),zero) 321 322 dpod=zero ! if not DDES 323 324 if(iles.lt.0) then ! for DES 325 dx=maxval(xl(e,:,1))-minval(xl(e,:,1)) 326 dy=maxval(xl(e,:,2))-minval(xl(e,:,2)) 327 dz=maxval(xl(e,:,3))-minval(xl(e,:,3)) 328 dmax=max(dx,max(dy,dz)) 329 dmax=0.65*dmax 330c 331c Next 5 lines are DDES 332c 333 qfac = sqrt(gradV(e,1,1)**2+gradV(e,1,2)**2+gradV(e,1,3)**2 334 & +gradV(e,2,1)**2+gradV(e,2,2)**2+gradV(e,2,3)**2 335 & +gradV(e,3,1)**2+gradV(e,3,2)**2+gradV(e,3,3)**2) 336 qfac = max(qfac,1.0e-10) 337 rd=SclrNN*saKappaP2Inv/(dwall(e)**2*qfac) 338 tanharg=tanh(8.0000000000000000d0*rd)**3 339 fd=one-tanharg 340c 341c LHS 342c 343 if (dwall(e) .ne. 0) 344 & rdp = saKappaP2Inv /(qfac * dwall(e)**2) 345 346 tanhp = (one-tanharg**2) 347 fdp = -512*rd**2*rdp*tanhp 348 dmaxf = max(zero,dwall(e)-dmax) 349 dwall(e)=dwall(e)-fd*dmaxf 350 dpod = -fdp*dmaxf/dwall(e) 351c DES97 dwall(e)=min(dmax,dwall(e)) 352 endif 353 354 if (dwall(e) .ne. 0) dP2Inv = 1.0 / dwall(e)**2 355 356 k2d2Inv = dP2Inv * saKappaP2Inv 357 358 viscInv = 1.0/datmat(1,2,1) 359 chi = SclrNN * viscInv !skipping the kiy kie for now 360 chiP2 = chi * chi 361 chiP3 = chi * chiP2 362 363 ep=zero 364 if(Sclr(e).gt.zero) ep=one 365 366 tmp = 1.0 / ( chiP3 + saCv1P3 ) 367 fv1 = chiP3 * tmp 368 fv1p = 3.0 * viscInv * saCv1P3 * chiP2 * tmp**2 369 370 tmp = 1.0 / (1.0 + chi * fv1) 371 fv2 = 1.0 - chi * tmp 372 fv2p = (chiP2 * fv1p - viscInv) * tmp**2 373 374 s = absVort(e) !prd(e,2) 375 tmp = SclrNN * k2d2Inv !eOkdP2 376 st = s + fv2 * tmp 377 stp = ep*k2d2Inv * fv2 *(one-two*SclrNN*dpod) + tmp*fv2p 378 379 r = zero 380 rp = zero 381 if (st .gt. epsM ) then 382 r = tmp / st 383 r = min( max(r, -8.0d0), 8.0d0) 384 rp = k2d2Inv / st**2 * (ep*st - SclrNN*(stp+st*two*dpod) 385 endif 386 rP5 = r * (r * r) ** 2 387 tmp = 1.0 + saCw2 * (rP5 - 1.0) 388 g = r * tmp 389 gp = rp * (tmp + 5.0 * saCw2 * rP5) 390 391 gP6 = (g * g * g) ** 2 392 tmp = 1.0 / (gP6 + saCw3P6) 393 tmp1 = ( (1.0 + saCw3P6) * tmp ) ** sixth 394 fw = g * tmp1 395 fwp = gp * tmp1 * saCw3P6 * tmp 396 if(abs(r).gt.7.0d0) fwp=zero 397 398 tmp1 = saCw1 * dP2Inv*SclrNN 399 prat = ep*saCb1*st !prodRat 400 srcRat(e)= prat - ep*fw * tmp1 401 402 tmp = zero 403 if(prat .lt. zero) tmp=prat 404 if(stp .lt. zero) tmp=tmp+SaCb1*stp*SclrNN 405 if(fw .gt. zero) tmp=tmp-two*tmp1*fw 406 if(fwp .gt. zero) tmp=tmp-tmp1*SclrNN*fwp 407 408 srcL(e) = -1.0*tmp 409 srcR(e) = srcRat(e)*SclrNN 410 enddo 411c 412c.... One source term has the form (beta_i)*(phi,_i). Since 413c the convective term has (u_i)*(phi,_i), it is useful to treat 414c beta_i as a "correction" to the velocity. In calculating the 415c stabilization terms, the new "modified" velocity (u_i-beta_i) is 416c then used in place of the pure velocity for stabilization terms, 417c and the source term sneaks into the RHS and LHS. Here, the term 418c is given by beta_i=(cb_2)*(phi,_i)/(sigma) 419c 420 421 tmp = saCb2 * saSigmaInv 422 uMod(:,1) = u1(:) - tmp * gradS(:,1) 423 uMod(:,2) = u2(:) - tmp * gradS(:,2) 424 uMod(:,3) = u3(:) - tmp * gradS(:,3) 425 426 elseif(iRANS.eq.-2) then ! K-epsilon 427c.... compute the global gradient of u 428c 429 gradV = zero 430 do n = 1, nshl 431c 432c du_i/dx_j 433c 434c i j indices match array where V is the velocity (u in our notes) 435 gradV(:,1,1) = gradV(:,1,1) + shg(:,n,1) * yl(:,n,2) 436 gradV(:,2,1) = gradV(:,2,1) + shg(:,n,1) * yl(:,n,3) 437 gradV(:,3,1) = gradV(:,3,1) + shg(:,n,1) * yl(:,n,4) 438c 439 gradV(:,1,2) = gradV(:,1,2) + shg(:,n,2) * yl(:,n,2) 440 gradV(:,2,2) = gradV(:,2,2) + shg(:,n,2) * yl(:,n,3) 441 gradV(:,3,2) = gradV(:,3,2) + shg(:,n,2) * yl(:,n,4) 442c 443 gradV(:,1,3) = gradV(:,1,3) + shg(:,n,3) * yl(:,n,2) 444 gradV(:,2,3) = gradV(:,2,3) + shg(:,n,3) * yl(:,n,3) 445 gradV(:,3,3) = gradV(:,3,3) + shg(:,n,3) * yl(:,n,4) 446c a j u a i 447c from our notes where we had N_{a,j} = dN_a/dx_j note that i is off by one because p was first in yl vector 448c 449 enddo 450 451 dwall = zero 452 do n = 1, nenl 453 dwall = dwall + shape_funct(:,n) * dwl(:,n) 454 enddo 455 456 kay(:)=zero 457 epsilon(:)=zero 458 do ii=1,npro 459 do jj = 1, nshl 460 kay(ii) = kay(ii) + shape_funct(ii,jj) * yl(ii,jj,6) 461 epsilon(ii) = epsilon(ii) 462 & + shape_funct(ii,jj) * yl(ii,jj,7) 463 enddo 464 enddo 465c no source term in the LHS yet 466 srcL(:)=zero 467 468 call elm3keps (kay, epsilon, 469 & dwall, gradV, 470 & srcRat, srcR, srcJac) 471 if(isclr.eq.1) srcL = srcJac(:,1) 472 if(isclr.eq.2) srcL = srcJac(:,4) 473 iadvdiff=0 ! scalar advection-diffusion flag 474 if(iadvdiff.eq.1)then 475 srcL(:)=zero 476 srcR(:)=zero 477 srcRat(:)=zero 478 endif 479c 480c.... No source terms with the form (beta_i)*(phi,_i) for K or E 481c 482 uMod(:,1) = u1(:) - zero 483 uMod(:,2) = u2(:) - zero 484 uMod(:,3) = u3(:) - zero 485 486 487 elseif (iLSet.ne.0) then 488 if (isclr.eq.1) then 489 490 srcR = zero 491 srcL = zero 492 493 elseif (isclr.eq.2) then !we are redistancing level-sets 494 495CAD If Sclr(:,1).gt.zero, result of sign_term function 1 496CAD If Sclr(:,1).eq.zero, result of sign_term function 0 497CAD If Sclr(:,1).lt.zero, result of sign_term function -1 498 499 sclr_ls = zero !zero out temp variable 500 501 do ii=1,npro 502 503 do jj = 1, nshl ! first find the value of levelset at point ii 504 505 sclr_ls(ii) = sclr_ls(ii) 506 & + shape_funct(ii,jj) * yl(ii,jj,6) 507 508 enddo 509 510 if (sclr_ls(ii) .gt. epsilon_ls)then 511 512 sign_levelset(ii) = one 513 514 elseif (abs(sclr_ls(ii)) .le. epsilon_ls)then 515 516 sign_levelset(ii) = sclr_ls(ii)/epsilon_ls 517 & + sin(pi*sclr_ls(ii)/epsilon_ls)/pi 518 519 elseif (sclr_ls(ii) .lt. epsilon_ls) then 520 521 sign_levelset(ii) = - one 522 523 endif 524 525 srcR(ii) = sign_levelset(ii) 526 527 enddo 528 529 srcL = zero 530 531cad The redistancing equation can be written in the following form 532cad 533cad d_{,t} + sign(phi)*( d_{,i}/|d_{,i}| )* d_{,i} = sign(phi) 534cad 535cad This is rewritten in the form 536cad 537cad d_{,t} + u * d_{,i} = sign(phi) 538cad 539 540CAD For the redistancing equation the "velocity" term is calculated as 541CAD follows 542 543 544 545 mytmp = sign_levelset / sqrt( gradS(:,1)*gradS(:,1) 546 & + gradS(:,2)*gradS(:,2) 547 & + gradS(:,3)*gradS(:,3)) 548 549 uMod(:,1) = mytmp * gradS(:,1) ! These are for the LHS 550 uMod(:,2) = mytmp * gradS(:,2) 551 uMod(:,3) = mytmp * gradS(:,3) 552 553 u1 = umod(:,1) ! These are for the RHS 554 u2 = umod(:,2) 555 u3 = umod(:,3) 556 557 558 endif ! close of scalar 2 of level set 559 560 else ! NOT turbulence and NOT level set so this is a simple 561 ! scalar. If your scalar equation has a source term 562 ! then add your own like the above but leave an unforced case 563 ! as an option like you see here 564 565 srcR = zero 566 srcL = zero 567 endif 568 569 return 570 end 571 572 573 574