1 subroutine itrdrv (y, ac, uold, x, 2 & iBC, BC, 3 & iper, ilwork, shp, 4 & shgl, shpb, shglb, 5 & ifath, velbar, nsons ) 6c 7c---------------------------------------------------------------------- 8c 9c This iterative driver is the semi-discrete, predictor multi-corrector 10c algorithm. It contains the Hulbert Generalized Alpha method which 11c is 2nd order accurate for Rho_inf from 0 to 1. The method can be 12c made first-order accurate by setting Rho_inf=-1. It uses a 13c GMRES iterative solver. 14c 15c working arrays: 16c y (nshg,ndof) : Y variables 17c x (nshg,nsd) : node coordinates 18c iBC (nshg) : BC codes 19c BC (nshg,ndofBC) : BC constraint parameters 20c iper (nshg) : periodicity table 21c 22c shape functions: 23c shp (nshape,ngauss) : interior element shape functions 24c shgl (nsd,nshape,ngauss) : local shape function gradients 25c shpb (nshapeb,ngaussb) : boundary element shape functions 26c shglb (nsd,nshapeb,ngaussb) : bdry. elt. shape gradients 27c 28c Zdenek Johan, Winter 1991. (Fortran 90) 29c---------------------------------------------------------------------- 30c 31 use pvsQbi !gives us splag (the spmass at the end of this run 32 use specialBC !gives us itvn 33 use timedata !allows collection of time series 34 use turbSA 35 36 include "common.h" 37 include "mpif.h" 38 include "auxmpi.h" 39 40c 41 dimension y(nshg,ndof), ac(nshg,ndof), 42 & yold(nshg,ndof), acold(nshg,ndof), 43 & x(numnp,nsd), iBC(nshg), 44 & BC(nshg,ndofBC), ilwork(nlwork), 45 & iper(nshg), uold(nshg,nsd) 46c 47 dimension res(nshg,nflow), BDiag(nshg,nflow,nflow), 48 & rest(nshg), solinc(nshg,ndof) 49c 50 dimension shp(MAXTOP,maxsh,MAXQPT), 51 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 52 & shpb(MAXTOP,maxsh,MAXQPT), 53 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 54 real*8 almit, alfit, gamit 55 dimension ifath(numnp), velbar(nfath,ndof), nsons(nfath) 56 real*8 rerr(nshg,10),ybar(nshg,ndof+8) ! 8 is for avg. of square as uu, vv, ww, pp, TT, uv, uw, and vw 57 integer, allocatable, dimension(:) :: ivarts 58 integer, allocatable, dimension(:) :: ivartsg 59 real*8, allocatable, dimension(:) :: vartssoln 60 real*8, allocatable, dimension(:) :: vartssolng 61 real*8, allocatable, dimension(:,:,:) :: vartsbuff 62! assuming three profiles to control (inlet, bottom FC and top FC) 63! store velocity profile set via BC 64 real*8 vbc_prof(nshg,3) 65 character*20 fname1, fmt1, fname2, fmt2 66 character*4 fname4c ! 4 characters 67 character*60 fvarts 68 character*5 cname 69 character*10 cname2 70 integer ifuncs(6), iarray(10) 71 72 real*8 elDw(numel) ! element average of DES d variable 73 74 real*8, allocatable, dimension(:,:) :: HBrg 75 real*8, allocatable, dimension(:) :: eBrg 76 real*8, allocatable, dimension(:) :: yBrg 77 real*8, allocatable, dimension(:) :: Rcos, Rsin 78c 79c Here are the data structures for sparse matrix GMRES 80c 81 integer, allocatable, dimension(:,:) :: rowp 82 integer, allocatable, dimension(:) :: colm 83 real*8, allocatable, dimension(:,:) :: lhsK 84 real*8, allocatable, dimension(:,:) :: EGmass 85 real*8, allocatable, dimension(:,:) :: EGmasst 86 87 inquire(file='xyzts.dat',exist=exts) 88 lskeep=lstep 89 if(exts) then 90 91 open(unit=626,file='xyzts.dat',status='old') 92 read(626,*) ntspts, freq, tolpt, iterat, varcod 93 call sTD ! sets data structures 94 do jj=1,ntspts ! read coordinate data where solution desired 95 read(626,*) ptts(jj,1),ptts(jj,2),ptts(jj,3) 96 enddo 97 close(626) 98 99 statptts(:,:) = 0 100 parptts(:,:) = zero 101 varts(:,:) = zero 102 103 allocate (ivarts(ntspts*ndof)) 104 allocate (ivartsg(ntspts*ndof)) 105 allocate (vartssoln(ntspts*ndof)) 106 allocate (vartssolng(ntspts*ndof)) 107 108 nbuff=ntout 109 allocate (vartsbuff(ntspts,ndof,nbuff)) 110 111 if (myrank .eq. master) then 112 do jj=1,ntspts 113 fvarts='varts/varts' 114 fvarts=trim(fvarts)//trim(cname2(jj)) 115 fvarts=trim(fvarts)//trim(cname2(lstep)) 116 fvarts=trim(fvarts)//'.dat' 117 fvarts=trim(fvarts) 118 open(unit=1000+jj, file=fvarts, status='unknown') 119 enddo 120 endif 121 122 endif 123 124c 125c 126c.... open history and aerodynamic forces files 127c 128 if (myrank .eq. master) then 129 open (unit=ihist, file=fhist, status='unknown') 130 open (unit=iforce, file=fforce, status='unknown') 131 endif 132c 133c 134c.... initialize 135c 136 ifuncs = 0 ! func. evaluation counter 137 istep = 0 138 ntotGM = 0 ! number of GMRES iterations 139 time = 0 140 yold = y 141 acold = ac 142 143 allocate(HBrg(Kspace+1,Kspace)) 144 allocate(eBrg(Kspace+1)) 145 allocate(yBrg(Kspace)) 146 allocate(Rcos(Kspace)) 147 allocate(Rsin(Kspace)) 148 149 if (mod(impl(1),100)/10 .eq. 1) then 150c 151c generate the sparse data fill vectors 152c 153 allocate (rowp(nshg,nnz)) 154 allocate (colm(nshg+1)) 155 call genadj(colm, rowp, icnt ) ! preprocess the adjacency list 156 157 nnz_tot=icnt ! this is exactly the number of non-zero 158 ! blocks on this proc 159 allocate (lhsK(nflow*nflow,nnz_tot)) 160 endif 161 if (mod(impl(1),100)/10 .eq. 3) then 162c 163c generate the ebe data fill vectors 164c 165 nedof=nflow*nshape 166 allocate (EGmass(numel,nedof*nedof)) 167 endif 168cc 169 170c.......................................... 171 rerr = zero 172 ybar(:,1:ndof) = y(:,1:ndof) 173 do idof=1,5 174 ybar(:,ndof+idof) = y(:,idof)*y(:,idof) 175 enddo 176 ybar(:,ndof+6) = y(:,1)*y(:,2) 177 ybar(:,ndof+7) = y(:,1)*y(:,3) 178 ybar(:,ndof+8) = y(:,2)*y(:,3) 179c......................................... 180 181 182 vbc_prof(:,1:3) = BC(:,3:5) 183 184c 185c.... loop through the time sequences 186c 187 do 3000 itsq = 1, ntseq 188 itseq = itsq 189c 190c.... set up the current parameters 191c 192 nstp = nstep(itseq) 193 nitr = niter(itseq) 194 LCtime = loctim(itseq) 195c 196 call itrSetup ( y, acold) 197 isclr=0 198 199 niter(itseq)=0 ! count number of flow solves in a step 200 ! (# of iterations) 201 do i=1,seqsize 202 if(stepseq(i).eq.0) niter(itseq)=niter(itseq)+1 203 enddo 204 nitr = niter(itseq) 205c 206c.... determine how many scalar equations we are going to need to solve 207c 208 nsclrsol=nsclr ! total number of scalars solved. At 209 ! some point we probably want to create 210 ! a map, considering stepseq(), to find 211 ! what is actually solved and only 212 ! dimension EGmasst to the appropriate 213 ! size. 214 215 if(nsclrsol.gt.0)allocate (EGmasst(numel*nshape*nshape 216 & ,nsclrsol)) 217 218c 219c.... loop through the time steps 220c 221 ttim(1) = REAL(secs(0.0)) / 100. 222 ttim(2) = secs(0.0) 223 224c tcorecp1 = REAL(secs(0.0)) / 100. 225c tcorewc1 = secs(0.0) 226 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 227 if(myrank.eq.0) then 228 tcorecp1 = TMRC() 229 endif 230 231 rmub=datmat(1,2,1) 232 if(rmutarget.gt.0) then 233 rmue=rmutarget 234 xmulfact=(rmue/rmub)**(1.0/nstp) 235 if(myrank.eq.0) then 236 write(*,*) 'viscosity will by multiplied by ', xmulfact 237 write(*,*) 'to bring it from ', rmub,' down to ', rmue 238 endif 239 datmat(1,2,1)=datmat(1,2,1)/xmulfact ! make first step right 240 else 241 rmue=datmat(1,2,1) ! keep constant 242 xmulfact=one 243 endif 244 if(iramp.eq.1) then 245 call initBCprofileScale(vbc_prof,BC,yold,x) 246! fix the yold values to the reset BC 247 call itrBC (yold, ac, iBC, BC, iper, ilwork) 248 isclr=1 ! fix scalar 249 call itrBCsclr(yold,ac,iBC,BC,iper,ilwork) 250 endif 251 252867 continue 253 do 2000 istp = 1, nstp 254 if(iramp.eq.1) 255 & call BCprofileScale(vbc_prof,BC,yold) 256 257c call rerun_check(stopjob) 258cc if(stopjob.ne.0) goto 2001 259c 260c Decay of scalars 261c 262 if(nsclr.gt.0 .and. tdecay.ne.1) then 263 yold(:,6:ndof)=y(:,6:ndof)*tdecay 264 BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*tdecay 265 endif 266 267 if(nosource.eq.1) BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*0.8 268 269 270c xi=istp*one/nstp 271c datmat(1,2,1)=rmub*(1.0-xi)+xi*rmue 272 datmat(1,2,1)=xmulfact*datmat(1,2,1) 273 274c.... if we have time varying boundary conditions update the values of BC. 275c these will be for time step n+1 so use lstep+1 276c 277 if(itvn.gt.0) call BCint((lstep+1)*Delt(1), shp, shgl, 278 & shpb, shglb, x, BC, iBC) 279 280 if(iLES.gt.0) then 281c 282c.... get dynamic model coefficient 283c 284 ilesmod=iLES/10 285c 286c digit bit set filter rule, 10 bit set model 287c 288 if (ilesmod.eq.0) then ! 0 < iLES < 10 => dyn. model calculated 289 ! at nodes based on discrete filtering 290 call getdmc (yold, shgl, shp, 291 & iper, ilwork, nsons, 292 & ifath, x) 293 endif 294 if (ilesmod .eq. 1) then ! 10 < iLES < 20 => dynamic-mixed 295 ! at nodes based on discrete filtering 296 call bardmc (yold, shgl, shp, 297 & iper, ilwork, 298 & nsons, ifath, x) 299 endif 300 if (ilesmod .eq. 2) then ! 20 < iLES < 30 => dynamic at quad 301 ! pts based on lumped projection filt. 302 call projdmc (yold, shgl, shp, 303 & iper, ilwork, x) 304 endif 305c 306 endif 307 308 309c 310c.... set traction BCs for modeled walls 311c 312 if (itwmod.ne.0) then !wallfn check 313 call asbwmod(yold, acold, x, BC, iBC, 314 & iper, ilwork, ifath, velbar) 315 endif 316c 317c.... -----------------------> predictor phase <----------------------- 318c 319 call itrPredict( yold, acold, y, ac ) 320 call itrBC (y, ac, iBC, BC, iper, ilwork) 321 isclr = zero 322 if (nsclr.gt.zero) then 323 do isclr=1,nsclr 324 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 325 enddo 326 endif 327c 328c.... --------------------> multi-corrector phase <-------------------- 329c 330 iter=0 331 ilss=0 ! this is a switch thrown on first solve of LS redistance 332c 333cHACK to make it keep solving RANS until tolerance is reached 334c 335 istop=0 336 iMoreRANS=0 337c 338c find the the RANS portion of the sequence 339c 340 do istepc=1,seqsize 341 if(stepseq(istepc).eq.10) iLastRANS=istepc 342 enddo 343 344 iseqStart= 1 3459876 continue 346c 347 do istepc=iseqStart,seqsize 348 icode=stepseq(istepc) 349 if(mod(icode,10).eq.0) then ! this is a solve 350 isolve=icode/10 351 if(isolve.eq.0) then ! flow solve (encoded as 0) 352c 353 etol=epstol(1) 354 iter = iter+1 355 ifuncs(1) = ifuncs(1) + 1 356c 357c.... reset the aerodynamic forces 358c 359 Force(1) = zero 360 Force(2) = zero 361 Force(3) = zero 362 HFlux = zero 363c 364c.... form the element data and solve the matrix problem 365c 366c.... explicit solver 367c 368 if (impl(itseq) .eq. 0) then 369 if (myrank .eq. master) 370 & call error('itrdrv ','impl ',impl(itseq)) 371 endif 372 if (mod(impl(1),100)/10 .eq. 1) then ! sparse solve 373c 374c.... preconditioned sparse matrix GMRES solver 375c 376 lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 377 iprec=lhs 378 nedof = nflow*nshape 379c write(*,*) 'lhs=',lhs 380 call SolGMRs (y, ac, yold, 381 & acold, x, 382 & iBC, BC, 383 & colm, rowp, lhsk, 384 & res, 385 & BDiag, hBrg, eBrg, 386 & yBrg, Rcos, Rsin, 387 & iper, ilwork, 388 & shp, shgl, 389 & shpb, shglb, solinc, 390 & rerr) 391 else if (mod(impl(1),100)/10 .eq. 2) then ! mfg solve 392c 393c.... preconditioned matrix-free GMRES solver 394c 395 lhs=0 396 iprec = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 397 nedof = 0 398 call SolMFG (y, ac, yold, 399 & acold, x, 400 & iBC, BC, 401 & res, 402 & BDiag, HBrg, eBrg, 403 & yBrg, Rcos, Rsin, 404 & iper, ilwork, 405 & shp, shgl, 406 & shpb, shglb, solinc, 407 & rerr) 408c 409 else if (mod(impl(1),100)/10 .eq. 3) then ! ebe solve 410c.... preconditioned ebe matrix GMRES solver 411c 412 413 lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 414 iprec = lhs 415 nedof = nflow*nshape 416c write(*,*) 'lhs=',lhs 417 call SolGMRe (y, ac, yold, 418 & acold, x, 419 & iBC, BC, 420 & EGmass, res, 421 & BDiag, HBrg, eBrg, 422 & yBrg, Rcos, Rsin, 423 & iper, ilwork, 424 & shp, shgl, 425 & shpb, shglb, solinc, 426 & rerr) 427 endif 428c 429 else ! solve a scalar (encoded at isclr*10) 430 ifuncs(isclr+2) = ifuncs(isclr+2) + 1 431 etol=epstol(isclr+1) 432 isclr=isolve 433 if((iLSet.eq.2).and.(ilss.eq.0) 434 & .and.(isclr.eq.2)) then 435 ilss=1 ! throw switch (once per step) 436c 437c... copy the first scalar at t=n+1 into the second scalar of the 438c... level sets 439c 440 y(:,6) = yold(:,6) + (y(:,6)-yold(:,6))/alfi 441 y(:,7) = y(:,6) 442 yold(:,7) = y(:,7) 443 ac(:,7) = zero 444c 445 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 446c 447c....store the flow alpha, gamma parameter values and assigm them the 448c....Backward Euler parameters to solve the second levelset scalar 449c 450 alfit=alfi 451 gamit=gami 452 almit=almi 453 alfi = 1 454 gami = 1 455 almi = 1 456 endif 457c 458 lhs = 1 - min(1,mod(ifuncs(isclr+2)-1, 459 & LHSupd(isclr+2))) 460 iprec = lhs 461 istop=0 462 call SolGMRSclr(y, ac, yold, 463 & acold, EGmasst(1,isclr), 464 & x, elDw, 465 & iBC, BC, 466 & rest, 467 & HBrg, eBrg, 468 & yBrg, Rcos, Rsin, 469 & iper, ilwork, 470 & shp, shgl, 471 & shpb, shglb, solinc(1,isclr+5)) 472c 473 endif ! end of scalar type solve 474c 475c 476c.... end of the multi-corrector loop 477c 478 1000 continue !check this 479 480 else ! this is an update (mod did not equal zero) 481 iupdate=icode/10 ! what to update 482 if(iupdate.eq.0) then !update flow 483 call itrCorrect ( y, ac, yold, acold, solinc) 484 call itrBC (y, ac, iBC, BC, iper, ilwork) 485 call tnanq(y, 5, 'y_updbc') 486c Elaine-SPEBC 487 if((irscale.ge.0).and.(myrank.eq.master)) then 488 call genscale(y, x, iBC) 489c call itrBC (y, ac, iBC, BC, iper, ilwork) 490 endif 491 else ! update scalar 492 isclr=iupdate !unless 493 if(iupdate.eq.nsclr+1) isclr=0 494 call itrCorrectSclr ( y, ac, yold, acold, 495 & solinc(1,isclr+5)) 496 if (ilset.eq.2 .and. isclr.eq.2) then 497 fct2=one/almi 498 fct3=one/alfi 499 acold(:,7) = acold(:,7) 500 & + (ac(:,7)-acold(:,7))*fct2 501 yold(:,7) = yold(:,7) 502 & + (y(:,7)-yold(:,7))*fct3 503 call itrBCSclr ( yold, acold, iBC, BC, 504 & iper, ilwork) 505 ac(:,7) = acold(:,7)*(one-almi/gami) 506 y(:,7) = yold(:,7) 507 ac(:,7) = zero 508 if (ivconstraint .eq. 1) then 509 & 510c ... applying the volume constraint 511c 512 call solvecon (y, x, iBC, BC, 513 & iper, ilwork, shp, shgl) 514c 515 endif ! end of volume constraint calculations 516 endif 517 call itrBCSclr ( y, ac, iBC, BC, iper, ilwork) 518 endif 519 endif !end of switch between solve or update 520 enddo ! loop over sequence in step 521 if((istop.lt.0).and.(iMoreRANS.lt.5)) then 522 iMoreRANS=iMoreRANS+1 523 if(myrank.eq.0) write(*,*) 'istop =', istop 524 iseqStart=iLastRANS 525 goto 9876 526 endif 527c 528c Find the solution at the end of the timestep and move it to old 529c 530c.... First to reassign the parameters for the original time integrator scheme 531c 532 if((iLSet.eq.2).and.(ilss.eq.1)) then 533 alfi =alfit 534 gami =gamit 535 almi =almit 536 endif 537 call itrUpdate( yold, acold, y, ac) 538 call itrBC (yold, acold, iBC, BC, iper,ilwork) 539c Elaine-SPEBC 540 if((irscale.ge.0).and.(myrank.eq.master)) then 541 call genscale(yold, x, iBC) 542c call itrBC (y, ac, iBC, BC, iper, ilwork) 543 endif 544 do isclr=1,nsclr 545 call itrBCSclr (yold, acold, iBC, BC, iper, ilwork) 546 enddo 547c 548 istep = istep + 1 549c 550c.... compute boundary fluxes and print out 551c 552 lstep = lstep + 1 553 ntoutv=max(ntout,100) 554 if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then 555c 556c here is where we save our averaged field. In some cases we want to 557c write it less frequently 558 if( (mod(lstep, ntoutv) .eq. 0) .and. 559 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 560 & ((nsonmax.eq.1).and.(iLES.gt.0)))) 561 & call rwvelb ('out ', velbar ,ifail) 562 563 call Bflux (yold, acold, x, 564 & shp, shgl, shpb, 565 & shglb, nodflx, ilwork) 566 endif 567c 568c.... end of the NSTEP and NTSEQ loops 569c 570 571c... dump TIME SERIES 572 573 if (exts) then 574 if (mod(lstep-1,freq).eq.0) then 575 576 if (numpe > 1) then 577 do jj = 1, ntspts 578 vartssoln((jj-1)*ndof+1:jj*ndof)=varts(jj,:) 579 ivarts=zero 580 enddo 581 do k=1,ndof*ntspts 582 if(vartssoln(k).ne.zero) ivarts(k)=1 583 enddo 584 call MPI_REDUCE(vartssoln, vartssolng, ndof*ntspts, 585 & MPI_DOUBLE_PRECISION, MPI_SUM, master, 586 & MPI_COMM_WORLD, ierr) 587 588 call MPI_REDUCE(ivarts, ivartsg, ndof*ntspts, 589 & MPI_INTEGER, MPI_SUM, master, 590 & MPI_COMM_WORLD, ierr) 591 592 if (myrank.eq.zero) then 593 do jj = 1, ntspts 594 595 indxvarts = (jj-1)*ndof 596 do k=1,ndof 597 if(ivartsg(indxvarts+k).ne.0) then ! none of the vartssoln(parts) were non zero 598 varts(jj,k)=vartssolng(indxvarts+k)/ 599 & ivartsg(indxvarts+k) 600 endif 601 enddo 602 enddo 603 endif !only on master 604 endif !only if numpe > 1 605 606 if( istp.eq. nstp) then !make sure incomplete buffers get purged. 607 icheck=mod(nstp,nbuff) 608 if(icheck.ne.0) nbuff=icheck 609 endif 610 611 if (myrank.eq.zero) then 612 k=mod(lstep,nbuff) 613 if(k.eq.0) k=nbuff 614 do jj = 1, ntspts 615 vartsbuff(jj,1:5,k)=varts(jj,1:5) 616 enddo 617 if(k.eq. nbuff) then 618 do jj = 1, ntspts 619 ifile = 1000+jj 620 do ibuf=1,nbuff 621 write(ifile,555) lstep-1 -nbuff+ibuf, 622 & (vartsbuff(jj,k,ibuf), k=1,5) ! assuming ndof to be 5 623 enddo ! buff empty 624 625c call flush(ifile) 626 enddo ! jj ntspts 627 endif !only dump when buffer full 628 endif !only on master 629 630 varts(:,:) = zero ! reset the array for next step 631 632 555 format(i6,5(2x,E12.5e2)) 633 634 endif 635 endif 636 637c 638c.... -------------------> error calculation <----------------- 639 if(ierrcalc.eq.1.or.ioybar.eq.1) then 640 tfact=one/istep 641 do idof=1,ndof 642 ybar(:,idof) =tfact*yold(:,idof) + 643 & (one-tfact)*ybar(:,idof) 644 enddo 645c....compute average 646c... ybar(:,ndof+1:ndof+8) is for avg. of square as uu, vv, ww, pp, TT, uv, uw, and vw 647 do idof=1,5 ! avg. of square for 5 flow variables 648 ybar(:,ndof+idof) = tfact*yold(:,idof)**2 + 649 & (one-tfact)*ybar(:,ndof+idof) 650 enddo 651 ybar(:,ndof+6) = tfact*yold(:,1)*yold(:,2) + !uv 652 & (one-tfact)*ybar(:,ndof+6) 653 ybar(:,ndof+7) = tfact*yold(:,1)*yold(:,3) + !uw 654 & (one-tfact)*ybar(:,ndof+7) 655 ybar(:,ndof+8) = tfact*yold(:,2)*yold(:,3) + !vw 656 & (one-tfact)*ybar(:,ndof+8) 657c... compute err 658c rerr(:, 7)=rerr(:, 7)+(yold(:,1)-ybar(:,1))**2 659c rerr(:, 8)=rerr(:, 8)+(yold(:,2)-ybar(:,2))**2 660c rerr(:, 9)=rerr(:, 9)+(yold(:,3)-ybar(:,3))**2 661c rerr(:,10)=rerr(:,10)+(yold(:,4)-ybar(:,4))**2 662 endif 663c........................................................................... 664 2000 continue 665 2001 continue 666 667 ttim(1) = REAL(secs(0.0)) / 100. - ttim(1) 668 ttim(2) = secs(0.0) - ttim(2) 669 670c tcorecp2 = REAL(secs(0.0)) / 100. 671c tcorewc2 = secs(0.0) 672 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 673 if(myrank.eq.0) then 674 tcorecp2 = TMRC() 675 write(6,*) 'T(core) cpu = ',tcorecp2-tcorecp1 676 endif 677 678c call wtime 679 680 3000 continue 681c 682c.... ----------------------> Post Processing <---------------------- 683c 684c.... print out the last step 685c 686 if ((irs .ge. 1) .and. ((mod(lstep, ntout) .ne. 0) .or. 687 & (nstp .eq. 0))) then 688 if( (mod(lstep, ntoutv) .eq. 0) .and. 689 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 690 & ((nsonmax.eq.1).and.(iLES.gt.0)))) 691 & call rwvelb ('out ', velbar ,ifail) 692 693 call Bflux (yold, acold, x, 694 & shp, shgl, shpb, 695 & shglb, nodflx, ilwork) 696 endif 697c 698 if(ierrcalc.eq.1) then 699c 700c.....smooth the error indicators 701c 702c ! errsmooth is currently available only for incompressible code 703c do i=1,ierrsmooth 704c call errsmooth( rerr, x, iper, ilwork, shp, shgl, iBC ) 705c end do 706 707 call write_error(myrank, lstep, nshg, 10, rerr ) 708 endif 709 710 if(ioybar.eq.1) then 711 call write_field(myrank,'a','ybar',4, 712 & ybar,'d',nshg,ndof+8,lstep) 713 endif 714 715 if(iRANS.lt.0 .and. idistcalc.eq.1) then 716 call write_field(myrank,'a','DESd',4, 717 & elDw,'d',numel,1,lstep) 718 719 call write_field(myrank,'a','dwal',4, 720 & d2wall,'d',nshg,1,lstep) 721 endif 722 723c 724c.... close history and aerodynamic forces files 725c 726 if (myrank .eq. master) then 727 close (ihist) 728 close (iforce) 729 if(exts) then 730 deallocate(ivarts) 731 deallocate(ivartsg) 732 deallocate(vartssoln) 733 deallocate(vartssolng) 734 do jj=1,ntspts 735 close(1000+jj) 736 enddo 737 endif 738 endif 739 close (iecho) 740 if(iabc==1) deallocate(acs) 741c 742c.... end 743c 744 return 745 end 746