1 subroutine ElmGMR (u, y, ac, x, 2 & shp, shgl, iBC, 3 & BC, shpb, shglb, 4 & res, iper, ilwork, 5 & rowp, colm, lhsK, 6 & lhsP, rerr, GradV) 7c 8c---------------------------------------------------------------------- 9c 10c This routine computes the LHS mass matrix, the RHS residual 11c vector, and the preconditioning matrix, for use with the GMRES 12c solver. 13c 14c Zdenek Johan, Winter 1991. (Fortran 90) 15c Chris Whiting, Winter 1998. (Matrix EBE-GMRES) 16c Alberto Figueroa, Winter 2004. CMM-FSI 17c Irene Vignon, Spring 2004. 18c---------------------------------------------------------------------- 19c 20 use pvsQbi ! brings in NABI 21 use stats ! 22 use pointer_data ! brings in the pointers for the blocked arrays 23 use local_mass 24 use timedata 25c 26 include "common.h" 27c 28 dimension y(nshg,ndof), ac(nshg,ndof), 29 & u(nshg,nsd), 30 & x(numnp,nsd), 31 & iBC(nshg), 32 & BC(nshg,ndofBC), 33 & res(nshg,nflow), 34 & iper(nshg) 35c 36 dimension shp(MAXTOP,maxsh,MAXQPT), 37 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 38 & shpb(MAXTOP,maxsh,MAXQPT), 39 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 40c 41 dimension qres(nshg,idflx), rmass(nshg) 42 dimension GradV(nshg,nsdsq) 43c 44 dimension ilwork(nlwork) 45 46 integer rowp(nshg*nnz), colm(nshg+1) 47 48 real*8 lhsK(9,nnz_tot), lhsP(4,nnz_tot) 49 50 real*8, allocatable, dimension(:,:,:,:) :: xKebe, xGoC 51 52 real*8 rerr(nshg,10) 53 54 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 55 real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:) 56 57 real*8 spmasstot(20), ebres(nshg) 58c 59c.... set up the timer 60c 61 62CAD call timer ('Elm_Form') 63c 64c.... --------------------> diffusive flux <-------------------- 65c 66c.... set up parameters 67c 68 ires = 1 69 70 if (idiff==1 .or. idiff==3 .or. isurf==1) then ! global reconstruction 71 ! of qdiff 72c 73c loop over element blocks for the global reconstruction 74c of the diffusive flux vector, q, and lumped mass matrix, rmass 75c 76 qres = zero 77 if(iter == nitr .and. icomputevort == 1 ) then 78 !write(*,*) 'iter:',iter,' - nitr:',nitr 79 !write(*,*) 'Setting GradV to zero' 80 GradV = zero 81 endif 82 rmass = zero 83 84 do iblk = 1, nelblk 85 iel = lcblk(1,iblk) 86 lelCat = lcblk(2,iblk) 87 lcsyst = lcblk(3,iblk) 88 iorder = lcblk(4,iblk) 89 nenl = lcblk(5,iblk) ! no. of vertices per element 90 nshl = lcblk(10,iblk) 91 mattyp = lcblk(7,iblk) 92 ndofl = lcblk(8,iblk) 93 nsymdl = lcblk(9,iblk) 94 npro = lcblk(1,iblk+1) - iel 95 ngauss = nint(lcsyst) 96c 97c.... compute and assemble diffusive flux vector residual, qres, 98c and lumped mass matrix, rmass 99 100 if(iter == nitr .and. icomputevort == 1 ) then 101 !write(*,*) 'Calling AsIqGradV' 102 call AsIqGradV (y, x, 103 & shp(lcsyst,1:nshl,:), 104 & shgl(lcsyst,:,1:nshl,:), 105 & mien(iblk)%p, 106 & GradV) 107 endif 108 call AsIq (y, x, 109 & shp(lcsyst,1:nshl,:), 110 & shgl(lcsyst,:,1:nshl,:), 111 & mien(iblk)%p, mxmudmi(iblk)%p, 112 & qres, rmass ) 113 enddo 114 115c 116c.... form the diffusive flux approximation 117c 118 call qpbc( rmass, qres, iBC, iper, ilwork ) 119 if(iter == nitr .and. icomputevort == 1 ) then 120 !write(*,*) 'Calling solveGradV' 121 call solveGradV( rmass, GradV, iBC, iper, ilwork ) 122 endif 123c 124 endif 125c 126c.... --------------------> interior elements <-------------------- 127c 128 res = zero 129 if (stsResFlg .ne. 1) then 130 flxID = zero 131 endif 132 133 if (lhs .eq. 1) then 134 lhsp = zero 135 lhsk = zero 136 endif 137c 138c.... loop over the element-blocks 139c 140 do iblk = 1, nelblk 141 iblock = iblk ! used in local mass inverse (p>2) 142 iblkts = iblk ! used in timeseries 143 iel = lcblk(1,iblk) 144 lelCat = lcblk(2,iblk) 145 lcsyst = lcblk(3,iblk) 146 iorder = lcblk(4,iblk) 147 nenl = lcblk(5,iblk) ! no. of vertices per element 148 nshl = lcblk(10,iblk) 149 mattyp = lcblk(7,iblk) 150 ndofl = lcblk(8,iblk) 151 nsymdl = lcblk(9,iblk) 152 npro = lcblk(1,iblk+1) - iel 153 inum = iel + npro - 1 154 ngauss = nint(lcsyst) 155c 156c.... allocate the element matrices 157c 158 allocate ( xKebe(npro,9,nshl,nshl) ) 159 allocate ( xGoC (npro,4,nshl,nshl) ) 160c 161c.... compute and assemble the residual and tangent matrix 162c 163 allocate (tmpshp(nshl,MAXQPT)) 164 allocate (tmpshgl(nsd,nshl,MAXQPT)) 165 166 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 167 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 168 169 call AsIGMR (y, ac, 170 & x, mxmudmi(iblk)%p, 171 & tmpshp, 172 & tmpshgl, 173 & mien(iblk)%p, 174 & res, 175 & qres, xKebe, 176 & xGoC, rerr) 177c 178c.... satisfy the BC's on the implicit LHS 179c 180 if (impl(1) .ne. 9 .and. lhs .eq. 1) then 181 if(ipord.eq.1) 182 & call bc3lhs (iBC, BC,mien(iblk)%p, xKebe) 183 call fillsparseI (mien(iblk)%p, 184 & xKebe, lhsK, 185 & xGoC, lhsP, 186 & rowp, colm) 187 endif 188 189 deallocate ( xKebe ) 190 deallocate ( xGoC ) 191 deallocate ( tmpshp ) 192 deallocate ( tmpshgl ) 193c 194c.... end of interior element loop 195c 196 enddo 197c$$$ if(ibksiz.eq.20 .and. iwrote.ne.789) then 198c$$$ do i=1,nshg 199c$$$ write(789,*) 'eqn block ',i 200c$$$ do j=colm(i),colm(i+1)-1 201c$$$ write(789,*) 'var block',rowp(j) 202c$$$ 203c$$$ do ii=1,3 204c$$$ write(789,111) (lhsK((ii-1)*3+jj,j),jj=1,3) 205c$$$ enddo 206c$$$ enddo 207c$$$ enddo 208c$$$ close(789) 209c$$$ iwrote=789 210c$$$ endif 211c$$$ 111 format(3(e14.7,2x)) 212c$$$c 213c.... add in lumped mass contributions if needed 214c 215 if((flmpr.ne.0).or.(flmpl.ne.0)) then 216 call lmassadd(ac,res,rowp,colm,lhsK,gmass) 217 endif 218 219 have_local_mass = 1 220c 221c.... time average statistics 222c 223 if ( stsResFlg .eq. 1 ) then 224 225 if (numpe > 1) then 226 call commu (stsVec, ilwork, nResDims , 'in ') 227 endif 228 do j = 1,nshg 229 if (btest(iBC(j),10)) then 230 i = iper(j) 231 stsVec(i,:) = stsVec(i,:) + stsVec(j,:) 232 endif 233 enddo 234c 235 do i = 1,nshg 236 stsVec(i,:) = stsVec(iper(i),:) 237 enddo 238 239 if (numpe > 1) then 240 call commu (stsVec, ilwork, nResDims , 'out') 241 endif 242 return 243 244 endif 245c 246c.... --------------------> boundary elements <-------------------- 247c 248c.... loop over the boundary elements 249c 250 do iblk = 1, nelblb 251c 252c.... set up the parameters 253c 254 iel = lcblkb(1,iblk) 255 lelCat = lcblkb(2,iblk) 256 lcsyst = lcblkb(3,iblk) 257 iorder = lcblkb(4,iblk) 258 nenl = lcblkb(5,iblk) ! no. of vertices per element 259 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 260 nshl = lcblkb(9,iblk) 261 nshlb = lcblkb(10,iblk) 262 mattyp = lcblkb(7,iblk) 263 ndofl = lcblkb(8,iblk) 264 npro = lcblkb(1,iblk+1) - iel 265 266 267 if(lcsyst.eq.3) lcsyst=nenbl 268c 269 if(lcsyst.eq.3 .or. lcsyst.eq.4) then 270 ngaussb = nintb(lcsyst) 271 else 272 ngaussb = nintb(lcsyst) 273 endif 274c 275c.... allocate the element matrices 276c 277 allocate ( xKebe(npro,9,nshl,nshl) ) 278 allocate ( xGoC (npro,4,nshl,nshl) ) 279 280c 281c.... compute and assemble the residuals corresponding to the 282c boundary integral 283c 284 allocate (tmpshpb(nshl,MAXQPT)) 285 allocate (tmpshglb(nsd,nshl,MAXQPT)) 286 287 tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:) 288 tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:) 289 290 call AsBMFG (u, y, 291 & ac, x, 292 & tmpshpb, 293 & tmpshglb, 294 & mienb(iblk)%p, mmatb(iblk)%p, 295 & miBCB(iblk)%p, mBCB(iblk)%p, 296 & res, xKebe) 297 298c 299c.... satisfy (again, for the vessel wall contributions) the BC's on the implicit LHS 300c 301c.... first, we need to make xGoC zero, since it doesn't have contributions from the 302c.... vessel wall elements 303 304 xGoC = zero 305 306 if (impl(1) .ne. 9 .and. lhs .eq. 1) then 307 if(ipord.eq.1) 308 & call bc3lhs (iBC, BC,mienb(iblk)%p, xKebe) 309 call fillsparseI (mienb(iblk)%p, 310 & xKebe, lhsK, 311 & xGoC, lhsP, 312 & rowp, colm) 313 endif 314 315 deallocate ( xKebe ) 316 deallocate ( xGoC ) 317 deallocate (tmpshpb) 318 deallocate (tmpshglb) 319c 320c.... end of boundary element loop 321c 322 enddo 323 324 if(ipvsq.ge.1) then 325c 326c.... pressure vs. resistance boundary condition sets pressure at 327c outflow to linearly increase as flow through that face increases 328c (routine is at bottom of this file) 329c 330 call ElmpvsQ (res,y,-1.0d0) 331 endif 332 333c 334c before the commu we need to rotate the residual vector for axisymmetric 335c boundary conditions (so that off processor periodicity is a dof add instead 336c of a dof combination). Take care of all nodes now so periodicity, like 337c commu is a simple dof add. 338c 339 if(iabc==1) !are there any axisym bc's 340 & call rotabc(res, iBC, 'in ') 341c 342c 343c.... --------------------> communications <------------------------- 344c 345 346 if (numpe > 1) then 347 call commu (res , ilwork, nflow , 'in ') 348 endif 349 350c 351c.... ----------------------> post processing <---------------------- 352c 353c.... satisfy the BCs on the residual 354c 355 call bc3Res (iBC, BC, res, iper, ilwork) 356c 357c.... return 358c 359c call timer ('Back ') 360 return 361 end 362 363 364!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 365!******************************************************************** 366!-------------------------------------------------------------------- 367 368 subroutine ElmGMRSclr (y, ac, x, 369 & shp, shgl, iBC, 370 & BC, shpb, shglb, 371 & res, iper, ilwork, 372 & rowp, colm, lhsS ) 373c 374c---------------------------------------------------------------------- 375c 376c This routine computes the LHS mass matrix, the RHS residual 377c vector, and the preconditioning matrix, for use with the GMRES 378c solver. 379c 380c---------------------------------------------------------------------- 381c 382 use pointer_data 383 use local_mass 384c 385 include "common.h" 386 include "mpif.h" 387c 388 dimension y(nshg,ndof), ac(nshg,ndof), 389 & x(numnp,nsd), iBC(nshg), 390 & BC(nshg,ndofBC), res(nshg), 391 & iper(nshg) 392c 393 dimension shp(MAXTOP,maxsh,MAXQPT), 394 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 395 & shpb(MAXTOP,maxsh,MAXQPT), 396 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 397c 398 dimension qres(nshg,nsd), rmass(nshg) 399c 400 integer ilwork(nlwork), rowp(nshg*nnz), colm(nshg+1) 401 402 real*8 lhsS(nnz_tot) 403 404 real*8, allocatable, dimension(:,:,:) :: xSebe 405c 406c.... set up the timer 407c 408 409CAD call timer ('Elm_Form') 410c 411c.... --------------------> diffusive flux <-------------------- 412c 413 ires = 1 414 415 if (idiff==1 .or. idiff==3) then ! global reconstruction of qdiff 416c 417c loop over element blocks for the global reconstruction 418c of the diffusive flux vector, q, and lumped mass matrix, rmass 419c 420 qres = zero 421 rmass = zero 422 423 do iblk = 1, nelblk 424 iel = lcblk(1,iblk) 425 lcsyst = lcblk(3,iblk) 426 nenl = lcblk(5,iblk) ! no. of vertices per element 427 nshl = lcblk(10,iblk) 428 mattyp = lcblk(7,iblk) 429 ndofl = lcblk(8,iblk) 430 npro = lcblk(1,iblk+1) - iel 431 432 ngauss = nint(lcsyst) 433c 434c.... compute and assemble diffusive flux vector residual, qres, 435c and lumped mass matrix, rmass 436 437 call AsIqSclr (y, x, 438 & shp(lcsyst,1:nshl,:), 439 & shgl(lcsyst,:,1:nshl,:), 440 & mien(iblk)%p, qres, 441 & rmass ) 442 443 enddo 444 445c 446c.... form the diffusive flux approximation 447c 448 call qpbcSclr ( rmass, qres, iBC, iper, ilwork ) 449c 450 endif 451c 452c.... --------------------> interior elements <-------------------- 453c 454 res = zero 455 spmass = zero 456 457 if (lhs .eq. 1) then 458 lhsS = zero 459 endif 460 461 if ((impl(1)/10) .eq. 0) then ! no flow solve so flxID was not zeroed 462 flxID = zero 463 endif 464c 465c.... loop over the element-blocks 466c 467 do iblk = 1, nelblk 468 iblock = iblk ! used in local mass inverse (p>2) 469 iel = lcblk(1,iblk) 470 lcsyst = lcblk(3,iblk) 471 nenl = lcblk(5,iblk) ! no. of vertices per element 472 nshl = lcblk(10,iblk) 473 ndofl = lcblk(8,iblk) 474 npro = lcblk(1,iblk+1) - iel 475 476 ngauss = nint(lcsyst) 477c 478c.... allocate the element matrices 479c 480 allocate ( xSebe(npro,nshl,nshl) ) 481c 482c.... compute and assemble the residual and tangent matrix 483c 484 call AsIGMRSclr(y, ac, 485 & x, 486 & shp(lcsyst,1:nshl,:), 487 & shgl(lcsyst,:,1:nshl,:), 488 & mien(iblk)%p, res, 489 & qres, xSebe, mxmudmi(iblk)%p ) 490c 491c.... satisfy the BC's on the implicit LHS 492c 493 if (impl(1) .ne. 9 .and. lhs .eq. 1) then 494 call fillsparseSclr (mien(iblk)%p, 495 & xSebe, lhsS, 496 & rowp, colm) 497 endif 498 499 deallocate ( xSebe ) 500c 501c.... end of interior element loop 502c 503 enddo 504 505c 506c.... add in lumped mass contributions if needed 507c 508 if((flmpr.ne.0).or.(flmpl.ne.0)) then 509 call lmassaddSclr(ac(:,isclr), res,rowp,colm,lhsS,gmass) 510 endif 511 512 have_local_mass = 1 513c 514c 515c call DtN routine which updates the flux to be consistent with the 516c current solution values. We will put the result in the last slot of 517c BC (we added a space in input.f). That way we can localize this 518c value to the boundary elements. This is important to keep from calling 519c the DtN evaluator more than once per node (it can be very expensive). 520c 521 if(idtn.eq.1) call DtN(iBC,BC,y) 522c 523c.... --------------------> boundary elements <-------------------- 524c 525c 526c.... loop over the boundary elements 527c 528 do iblk = 1, nelblb 529c 530c.... set up the parameters 531c 532 iel = lcblkb(1,iblk) 533 lcsyst = lcblkb(3,iblk) 534 nenl = lcblkb(5,iblk) ! no. of vertices per element 535 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 536 nshl = lcblkb(9,iblk) 537 nshlb = lcblkb(10,iblk) 538 ndofl = lcblkb(8,iblk) 539 npro = lcblkb(1,iblk+1) - iel 540 541 if(lcsyst.eq.3) lcsyst=nenbl 542 if(lcsyst.eq.3 .or. lcsyst.eq.4) then 543 ngaussb = nintb(lcsyst) 544 else 545 ngaussb = nintb(lcsyst) 546 endif 547c 548c localize the dtn boundary condition 549c 550 551 if(idtn.eq.1) call dtnl( iBC, BC, mienb(iblk)%p, 552 & miBCB(iblk)%p, mBCB(iblk)%p) 553 554c 555c.... compute and assemble the residuals corresponding to the 556c boundary integral 557c 558 call AsBSclr (y, x, 559 & shpb(lcsyst,1:nshl,:), 560 & shglb(lcsyst,:,1:nshl,:), 561 & mienb(iblk)%p, mmatb(iblk)%p, 562 & miBCB(iblk)%p, mBCB(iblk)%p, 563 & res) 564c 565c.... end of boundary element loop 566c 567 enddo 568c 569c 570c.... --------------------> communications <------------------------- 571c 572 573 if (numpe > 1) then 574 call commu (res , ilwork, 1 , 'in ') 575 endif 576 577c 578c.... ----------------------> post processing <---------------------- 579c 580c.... satisfy the BCs on the residual 581c 582 call bc3ResSclr (iBC, res, iper, ilwork) 583c 584c.... return 585c 586CAD call timer ('Back ') 587 return 588 end 589 590 591c 592c....routine to compute and return the flow rates for coupled surfaces of a given type 593c 594 subroutine GetFlowQ (qsurf,y,srfIdList,numSrfs) 595 596 use pvsQbi ! brings in NABI 597c 598 include "common.h" 599 include "mpif.h" 600 601 real*8 y(nshg,3) 602 real*8 qsurf(0:MAXSURF),qsurfProc(0:MAXSURF) 603 integer numSrfs, irankCoupled, srfIdList(0:MAXSURF) 604 605c note we only need the first three entries (u) from y 606 607 qsurfProc=zero 608 609 do i = 1,nshg 610 611 if(numSrfs.gt.zero) then 612 do k = 1,numSrfs 613 irankCoupled = 0 614 if (srfIdList(k).eq.ndsurf(i)) then 615 irankCoupled=k 616 do j = 1,3 617 qsurfProc(irankCoupled) = qsurfProc(irankCoupled) 618 & + NABI(i,j)*y(i,j) 619 enddo 620 exit 621 endif 622 enddo 623 endif 624 625 enddo 626c 627c at this point, each qsurf has its "nodes" contributions to Q 628c accumulated into qsurf. Note, because NABI is on processor this 629c will NOT be Q for the surface yet 630c 631c.... reduce integrated Q for each surface, push on qsurf 632c 633 npars=MAXSURF+1 634 if(impistat.eq.1) then 635 iAllR = iAllR+1 636 elseif(impistat.eq.2) then 637 iAllR = iAllR+1 638 endif 639 if(impistat2.eq.1) call MPI_BARRIER (MPI_COMM_WORLD, ierr) 640 if(impistat.gt.0) rmpitmr = TMRC() 641 call MPI_ALLREDUCE (qsurfProc, qsurf(:), npars, 642 & MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr) 643 if(impistat.eq.1) then 644 rAllR = rAllR+TMRC()-rmpitmr 645 elseif(impistat.eq.2) then 646 rAllRScal = rAllRScal+TMRC()-rmpitmr 647 endif 648 649c 650c.... return 651c 652 return 653 end 654 655 656 657c 658c... routine to couple pressure with flow rate for each coupled surface 659c 660 subroutine ElmpvsQ (res,y,sign) 661 662 use pvsQbi ! brings in NABI 663 use convolImpFlow !brings in the current part of convol coef for imp BC 664 use convolRCRFlow !brings in the current part of convol coef for RCR BC 665 666 include "common.h" 667 include "mpif.h" 668 669 real*8 res(nshg,ndof),y(nshg,3) 670 real*8 p(0:MAXSURF) 671 integer irankCoupled 672 673c 674c... get p for the resistance BC 675c 676 if(numResistSrfs.gt.zero) then 677 call GetFlowQ(p,y,nsrflistResist,numResistSrfs) !Q pushed into p but at this point 678 ! p is just the full Q for each surface 679 p(:)=sign*p(:)*ValueListResist(:) ! p=QR now we have the true pressure on each 680 ! outflow surface. Note sign is -1 681 ! for RHS, +1 for LHS 682c 683c.... multiply it by integral NA n_i 684c 685 do i = 1,nshg 686 do k = 1,numResistSrfs 687 irankCoupled = 0 688 if (nsrflistResist(k).eq.ndsurf(i)) then 689 irankCoupled=k 690 res(i,1:3)=res(i,1:3) + p(irankCoupled)*NABI(i,1:3) 691 exit 692 endif 693 enddo 694 enddo 695 696 endif !end of coupling for Resistance BC 697 698 699c 700c... get p for the impedance BC 701c 702 if(numImpSrfs.gt.zero) then 703 call GetFlowQ(p,y,nsrflistImp,numImpSrfs) !Q pushed into p but at this point 704 ! p is just the full Q for each surface 705 do j = 1,numImpSrfs 706 if(sign.lt.zero) then ! RHS so -1 707 p(j)= sign*(poldImp(j) + p(j)*ImpConvCoef(ntimeptpT+2,j)) !pressure p=pold+ Qbeta 708 elseif(sign.gt.zero) then ! LHS so sign is positive 709 p(j)= sign*p(j)*ImpConvCoef(ntimeptpT+2,j) 710 endif 711 enddo 712 713c 714c.... multiply it by integral NA n_i 715c 716 do i = 1,nshg 717 do k = 1,numImpSrfs 718 irankCoupled = 0 719 if (nsrflistImp(k).eq.ndsurf(i)) then 720 irankCoupled=k 721 res(i,1:3)=res(i,1:3) + p(irankCoupled)*NABI(i,1:3) 722 exit 723 endif 724 enddo 725 enddo 726 727 endif !end of coupling for Impedance BC 728c 729c... get p for the RCR BC 730c 731 if(numRCRSrfs.gt.zero) then 732 call GetFlowQ(p,y,nsrflistRCR,numRCRSrfs) !Q pushed into p but at this point 733 ! p is just the full Q for each surface 734 do j = 1,numRCRSrfs 735 if(sign.lt.zero) then ! RHS so -1 736 p(j)= sign*(poldRCR(j) + p(j)*RCRConvCoef(lstep+2,j)) !pressure p=pold+ Qbet 737 p(j)= p(j) - HopRCR(j) ! H operator contribution 738 elseif(sign.gt.zero) then ! LHS so sign is positive 739 p(j)= sign*p(j)*RCRConvCoef(lstep+2,j) 740 endif 741 enddo 742 743c 744c.... multiply it by integral NA n_i 745c 746 do i = 1,nshg 747 do k = 1,numRCRSrfs 748 irankCoupled = 0 749 if (nsrflistRCR(k).eq.ndsurf(i)) then 750 irankCoupled=k 751 res(i,1:3)=res(i,1:3) + p(irankCoupled)*NABI(i,1:3) 752 exit 753 endif 754 enddo 755 enddo 756 757 endif !end of coupling for RCR BC 758 759 return 760 end 761 762 763 764 765 766 767 768