1 subroutine itrdrv (y, ac, 2 & uold, x, 3 & iBC, BC, 4 & iper, ilwork, shp, 5 & shgl, shpb, shglb, 6 & ifath, velbar, nsons ) 7c 8c---------------------------------------------------------------------- 9c 10c This iterative driver is the semi-discrete, predictor multi-corrector 11c algorithm. It contains the Hulbert Generalized Alpha method which 12c is 2nd order accurate for Rho_inf from 0 to 1. The method can be 13c made first-order accurate by setting Rho_inf=-1. It uses CGP and 14c GMRES iterative solvers. 15c 16c working arrays: 17c y (nshg,ndof) : Y variables 18c x (nshg,nsd) : node coordinates 19c iBC (nshg) : BC codes 20c BC (nshg,ndofBC) : BC constraint parameters 21c iper (nshg) : periodicity table 22c 23c 24c Zdenek Johan, Winter 1991. (Fortran 90) 25c Alberto Figueroa, Winter 2004. CMM-FSI 26c Irene Vignon, Fall 2004. Impedance BC 27c---------------------------------------------------------------------- 28c 29 use pvsQbi !gives us splag (the spmass at the end of this run 30 use specialBC !gives us itvn 31 use timedata !allows collection of time series 32 use convolImpFlow !for Imp bc 33 use convolRCRFlow !for RCR bc 34!MR CHANGE 35 use turbsa ! used to access d2wall 36!MR CHANGE END 37 use iso_c_binding 38 39c use readarrays !reads in uold and acold 40 41 include "common.h" 42 include "mpif.h" 43 include "auxmpi.h" 44c 45 46 47 real*8 y(nshg,ndof), ac(nshg,ndof), 48 & yold(nshg,ndof), acold(nshg,ndof), 49 & u(nshg,nsd), uold(nshg,nsd), 50 & x(numnp,nsd), solinc(nshg,ndof), 51 & BC(nshg,ndofBC), tf(nshg,ndof), 52 & GradV(nshg,nsdsq) 53 54c 55 real*8 res(nshg,ndof) 56c 57 real*8 shp(MAXTOP,maxsh,MAXQPT), 58 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 59 & shpb(MAXTOP,maxsh,MAXQPT), 60 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 61c 62 integer rowp(nshg,nnz), colm(nshg+1), 63 & iBC(nshg), 64 & ilwork(nlwork), 65 & iper(nshg), ifuncs(6) 66 67 real*8 vbc_prof(nshg,3) 68 69 integer stopjob 70 character*10 cname2 71 character*5 cname 72c 73c stuff for dynamic model s.w.avg and wall model 74c 75 dimension ifath(numnp), velbar(nfath,ndof), nsons(nfath) 76 77 dimension wallubar(2),walltot(2) 78c 79c.... For Farzin's Library 80c 81 integer eqnType, prjFlag, presPrjFlag, verbose 82c 83 real*8, allocatable, dimension(:,:) :: aperm, atemp, atempS 84 real*8, allocatable, dimension(:,:,:) :: apermS 85 86 real*8, allocatable, dimension(:,:) :: lhsP, lhsK, lhsS 87 real*8 almit, alfit, gamit 88c 89 character*1024 servername 90 character*20 fname1,fmt1 91 character*20 fname2,fmt2 92 character*60 fnamepold, fvarts 93 character*4 fname4c ! 4 characters 94 integer iarray(50) ! integers for headers 95 integer isgn(ndof), isgng(ndof) 96 97!MR CHANGE 98! real*8 rerr(nshg,10), ybar(nshg,13) ! including 7 sq. terms with 3 cross terms of uv, uw and vw 99! real*8 rerr(nshg,10), ybar(nshg,12) ! including 7 sq. terms with 3 cross terms of uv, uw and vw 100 real*8 rerr(nshg,10) 101 real*8, allocatable, dimension(:,:) :: ybar, strain, vorticity 102 real*8, allocatable, dimension(:,:) :: wallssVec, wallssVecbar 103!MR CHANGE 104 105 real*8 tcorecp(2), tcorecpscal(2) 106 107 integer, allocatable, dimension(:) :: ivarts 108 integer, allocatable, dimension(:) :: ivartsg 109 integer, allocatable, dimension(:) :: iv_rank 110 real*8, allocatable, dimension(:) :: vartssoln 111 real*8, allocatable, dimension(:) :: vartssolng 112 real*8, allocatable, dimension(:,:,:) :: yphbar 113 real*8 CFLworst(numel) 114 115!MR CHANGE 116 integer :: iv_rankpernode, iv_totnodes, iv_totcores 117 integer :: iv_node, iv_core, iv_thread 118!MR CHANGE 119 120 impistat = 0 121 impistat2 = 0 122 iISend = 0 123 iISendScal = 0 124 iIRecv = 0 125 iIRecvScal = 0 126 iWaitAll = 0 127 iWaitAllScal = 0 128 iAllR = 0 129 iAllRScal = 0 130 rISend = zero 131 rISendScal = zero 132 rIRecv = zero 133 rIRecvScal = zero 134 rWaitAll = zero 135 rWaitAllScal = zero 136 rAllR = zero 137 rAllRScal = zero 138 rCommu = zero 139 rCommuScal = zero 140 141 call initmemstat() 142 143c 144c hack SA variable 145c 146cHack BC(:,7)=BC(:,7)*0.001 147cHack if(lstep.eq.0) y(:,6)=y(:,6)*0.001 148 call SolverLicenseServer(servername) 149c 150c only master should be verbose 151c 152 153 if(numpe.gt.0 .and. myrank.ne.master)iverbose=0 154c 155 156 lskeep=lstep 157 158 inquire(file='xyzts.dat',exist=exts) 159 if(exts) then 160 161 open(unit=626,file='xyzts.dat',status='old') 162 read(626,*) ntspts, freq, tolpt, iterat, varcod 163 call sTD ! sets data structures 164 165 do jj=1,ntspts ! read coordinate data where solution desired 166 read(626,*) ptts(jj,1),ptts(jj,2),ptts(jj,3) 167 enddo 168 close(626) 169 170 statptts(:,:) = 0 171 parptts(:,:) = zero 172 varts(:,:) = zero 173 174 allocate (ivarts(ntspts*ndof)) 175 allocate (ivartsg(ntspts*ndof)) 176 allocate (iv_rank(ntspts)) 177 allocate (vartssoln(ntspts*ndof)) 178 allocate (vartssolng(ntspts*ndof)) 179 180 iv_rankpernode = iv_rankpercore*iv_corepernode 181 iv_totnodes = numpe/iv_rankpernode 182 iv_totcores = iv_corepernode*iv_totnodes 183 if (myrank .eq. 0) then 184 write(*,*) 'Info for probes:' 185 write(*,*) ' Ranks per core:',iv_rankpercore 186 write(*,*) ' Cores per node:',iv_corepernode 187 write(*,*) ' Ranks per node:',iv_rankpernode 188 write(*,*) ' Total number of nodes:',iv_totnodes 189 write(*,*) ' Total number of cores',iv_totcores 190 endif 191 192! if (myrank .eq. numpe-1) then 193 do jj=1,ntspts 194 195 ! Compute the adequate rank which will take care of probe jj 196 jjm1 = jj-1 197 iv_node = (iv_totnodes-1)-mod(jjm1,iv_totnodes) 198 iv_core = (iv_corepernode-1) - mod((jjm1 - 199 & mod(jjm1,iv_totnodes))/iv_totnodes,iv_corepernode) 200 iv_thread = (iv_rankpercore-1) - mod((jjm1- 201 & (mod(jjm1,iv_totcores)))/iv_totcores,iv_rankpercore) 202 iv_rank(jj) = iv_node*iv_rankpernode 203 & + iv_core*iv_rankpercore 204 & + iv_thread 205 206 if(myrank == 0) then 207 write(*,*) ' Probe', jj, 'handled by rank', 208 & iv_rank(jj), ' on node', iv_node 209 endif 210 211 ! Verification just in case 212 if(iv_rank(jj) .lt.0 .or. iv_rank(jj) .ge. numpe) then 213 write(*,*) 'WARNING: iv_rank(',jj,') is ', iv_rank(jj), 214 & ' and reset to numpe-1' 215 iv_rank(jj) = numpe-1 216 endif 217 218 ! Open the varts files 219 if(myrank == iv_rank(jj)) then 220 fvarts='varts/varts' 221 fvarts=trim(fvarts)//trim(cname2(jj)) 222 fvarts=trim(fvarts)//trim(cname2(lstep)) 223 fvarts=trim(fvarts)//'.dat' 224 fvarts=trim(fvarts) 225 open(unit=1000+jj, file=fvarts, status='unknown') 226 endif 227 enddo 228! endif 229 230 endif 231c 232c.... open history and aerodynamic forces files 233c 234 if (myrank .eq. master) then 235 open (unit=ihist, file=fhist, status='unknown') 236 open (unit=iforce, file=fforce, status='unknown') 237 open (unit=76, file="fort.76", status='unknown') 238 if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then 239 fnamepold = 'pold' 240 fnamepold = trim(fnamepold)//trim(cname2(lstep)) 241 fnamepold = trim(fnamepold)//'.dat' 242 fnamepold = trim(fnamepold) 243 open (unit=8177, file=fnamepold, status='unknown') 244 endif 245 endif 246c 247c.... initialize 248c 249 ifuncs(:) = 0 ! func. evaluation counter 250 istep = 0 251 yold = y 252 acold = ac 253 254 rerr = zero 255 256 if(ierrcalc.eq.1 .or. ioybar.eq.1) then ! we need ybar for error too 257 if (ivort == 1) then 258 allocate(ybar(nshg,17)) ! more space for vorticity if requested 259 else 260 allocate(ybar(nshg,13)) 261 endif 262 ybar = zero ! Initialize ybar to zero, which is essential 263 endif 264 265 if(ivort == 1) then 266 allocate(strain(nshg,6)) 267 allocate(vorticity(nshg,5)) 268 endif 269 270 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 271 allocate(wallssVec(nshg,3)) 272 if (ioybar .eq. 1) then 273 allocate(wallssVecbar(nshg,3)) 274 wallssVecbar = zero ! Initialization important if mean wss computed 275 endif 276 endif 277 278! both nstepsincycle and nphasesincycle needs to be set 279 if(nstepsincycle.eq.0) nphasesincycle = 0 280 if(nphasesincycle.ne.0) then 281! & allocate(yphbar(nshg,5,nphasesincycle)) 282 if (ivort == 1) then 283 allocate(yphbar(nshg,15,nphasesincycle)) ! more space for vorticity 284 else 285 allocate(yphbar(nshg,11,nphasesincycle)) 286 endif 287 yphbar = zero 288 endif 289 290!MR CHANGE END 291 292 vbc_prof(:,1:3) = BC(:,3:5) 293 if(iramp.eq.1) then 294 call BCprofileInit(vbc_prof,x) 295 endif 296 297c 298c.... ---------------> initialize Farzin's Library <--------------- 299c 300c.... assign parameter values 301c 302 do i = 1, 100 303 numeqns(i) = i 304 enddo 305 nKvecs = Kspace 306 prjFlag = iprjFlag 307 presPrjFlag = ipresPrjFlag 308 verbose = iverbose 309c 310c.... determine how many scalar equations we are going to need to solve 311c 312 nsolt=mod(impl(1),2) ! 1 if solving temperature 313 nsclrsol=nsolt+nsclr ! total number of scalars solved At 314 ! some point we probably want to create 315 ! a map, considering stepseq(), to find 316 ! what is actually solved and only 317 ! dimension lhs to the appropriate 318 ! size. (see 1.6.1 and earlier for a 319 ! "failed" attempt at this). 320 321 322 nsolflow=mod(impl(1),100)/10 ! 1 if solving flow 323 324c 325c.... Now, call Farzin's lesNew routine to initialize 326c memory space 327c 328 call genadj(colm, rowp, icnt ) ! preprocess the adjacency list 329 330 nnz_tot=icnt ! this is exactly the number of non-zero blocks on 331 ! this proc 332 333 if (nsolflow.eq.1) then 334 lesId = numeqns(1) 335 eqnType = 1 336 nDofs = 4 337 call myfLesNew( lesId, 41994, 338 & eqnType, 339 & nDofs, minIters, maxIters, 340 & nKvecs, prjFlag, nPrjs, 341 & presPrjFlag, nPresPrjs, epstol(1), 342 & prestol, verbose, statsflow, 343 & nPermDims, nTmpDims, servername ) 344 345 allocate (aperm(nshg,nPermDims)) 346 allocate (atemp(nshg,nTmpDims)) 347 allocate (lhsP(4,nnz_tot)) 348 allocate (lhsK(9,nnz_tot)) 349 350 call readLesRestart( lesId, aperm, nshg, myrank, lstep, 351 & nPermDims ) 352 353 else 354 nPermDims = 0 355 nTempDims = 0 356 endif 357 358 359 if(nsclrsol.gt.0) then 360 do isolsc=1,nsclrsol 361 lesId = numeqns(isolsc+1) 362 eqnType = 2 363 nDofs = 1 364 presPrjflag = 0 365 nPresPrjs = 0 366 prjFlag = 1 367 indx=isolsc+2-nsolt ! complicated to keep epstol(2) for 368 ! temperature followed by scalars 369 call myfLesNew( lesId, 41994, 370 & eqnType, 371 & nDofs, minIters, maxIters, 372 & nKvecs, prjFlag, nPrjs, 373 & presPrjFlag, nPresPrjs, epstol(indx), 374 & prestol, verbose, statssclr, 375 & nPermDimsS, nTmpDimsS, servername ) 376 enddo 377c 378c Assume all scalars have the same size needs 379c 380 allocate (apermS(nshg,nPermDimsS,nsclrsol)) 381 allocate (atempS(nshg,nTmpDimsS)) !they can all share this 382 allocate (lhsS(nnz_tot,nsclrsol)) 383c 384c actually they could even share with atemp but leave that for later 385c 386 else 387 nPermDimsS = 0 388 nTmpDimsS = 0 389 endif 390c 391c... prepare lumped mass if needed 392c 393 if((flmpr.ne.0).or.(flmpl.ne.0)) call genlmass(x, shp,shgl) 394c 395c.... -----------------> End of initialization <----------------- 396c 397c.....open the necessary files to gather time series 398c 399 lstep0 = lstep+1 400 nsteprcr = nstep(1)+lstep 401c 402c.... loop through the time sequences 403c 404 405 406 do 3000 itsq = 1, ntseq 407 itseq = itsq 408 409CAD tcorecp1 = second(0) 410CAD tcorewc1 = second(-1) 411c 412c.... set up the time integration parameters 413c 414 nstp = nstep(itseq) 415 nitr = niter(itseq) 416 LCtime = loctim(itseq) 417 dtol(:)= deltol(itseq,:) 418 419 call itrSetup ( y, acold ) 420c 421c...initialize the coefficients for the impedance convolution, 422c which are functions of alphaf so need to do it after itrSetup 423 if(numImpSrfs.gt.zero) then 424 call calcImpConvCoef (numImpSrfs, ntimeptpT) 425 endif 426c 427c...initialize the initial condition P(0)-RQ(0)-Pd(0) for RCR BC 428c need ndsurf so should be after initNABI 429 if(numRCRSrfs.gt.zero) then 430 call calcRCRic(y,nsrflistRCR,numRCRSrfs) 431 endif 432c 433c find the last solve of the flow in the step sequence so that we will 434c know when we are at/near end of step 435c 436c ilast=0 437 nitr=0 ! count number of flow solves in a step (# of iterations) 438 do i=1,seqsize 439 if(stepseq(i).eq.0) nitr=nitr+1 440 enddo 441 442 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 443 tcorecp(:) = zero ! used in solfar.f (solflow) 444 tcorecpscal(:) = zero ! used in solfar.f (solflow) 445 if(myrank.eq.0) then 446 tcorecp1 = TMRC() 447 endif 448 449c 450c.... loop through the time steps 451c 452 istop=0 453 rmub=datmat(1,2,1) 454 if(rmutarget.gt.0) then 455 rmue=rmutarget 456 else 457 rmue=datmat(1,2,1) ! keep constant 458 endif 459 460 if(iramp.eq.1) then 461 call BCprofileScale(vbc_prof,BC,yold) ! fix the yold values to the reset BC 462 isclr=1 ! fix scalar 463 do isclr=1,nsclr 464 call itrBCSclr(yold,ac,iBC,BC,iper,ilwork) 465 enddo 466 endif 467 468 do 2000 istp = 1, nstp 469 if(iramp.eq.1) 470 & call BCprofileScale(vbc_prof,BC,yold) 471 472 call rerun_check(stopjob) 473 if(myrank.eq.master) write(*,*) 474 & 'stopjob,lstep,istep', stopjob,lstep,istep 475 if(stopjob.eq.lstep) then 476 stopjob=-2 ! this is the code to finish 477 if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then 478 if(myrank.eq.master) write(*,*) 479 & 'line 473 says last step written so exit' 480 goto 2002 ! the step was written last step so just exit 481 else 482 if(myrank.eq.master) 483 & write(*,*) 'line 473 says last step not written' 484 istep=nstp ! have to do this so that solution will be written 485 goto 2001 486 endif 487 endif 488 489 xi=istp*1.0/nstp 490 datmat(1,2,1)=rmub*(1.0-xi)+xi*rmue 491c write(*,*) "current mol. visc = ", datmat(1,2,1) 492c.... if we have time varying boundary conditions update the values of BC. 493c these will be for time step n+1 so use lstep+1 494c 495 if(itvn.gt.0) call BCint((lstep+1)*Delt(1), shp, shgl, 496 & shpb, shglb, x, BC, iBC) 497 498c 499c ... calculate the pressure contribution that depends on the history for the Imp. BC 500c 501 if(numImpSrfs.gt.0) then 502 call pHist(poldImp,QHistImp,ImpConvCoef, 503 & ntimeptpT,numImpSrfs) 504 if (myrank.eq.master) 505 & write(8177,*) (poldImp(n), n=1,numImpSrfs) 506 endif 507c 508c ... calc the pressure contribution that depends on the history for the RCR BC 509c 510 if(numRCRSrfs.gt.0) then 511 call CalcHopRCR (Delt(itseq), lstep, numRCRSrfs) 512 call CalcRCRConvCoef(lstep,numRCRSrfs) 513 call pHist(poldRCR,QHistRCR,RCRConvCoef,nsteprcr, 514 & numRCRSrfs) 515 if (myrank.eq.master) 516 & write(8177,*) (poldRCR(n), n=1,numRCRSrfs) 517 endif 518c 519c Decay of scalars 520c 521 if(nsclr.gt.0 .and. tdecay.ne.1) then 522 yold(:,6:ndof)=y(:,6:ndof)*tdecay 523 BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*tdecay 524 endif 525 526 if(nosource.eq.1) BC(:,7:6+nsclr)= BC(:,7:6+nsclr)*0.8 527 528 529 if(iLES.gt.0) then !complicated stuff has moved to 530 !routine below 531 call lesmodels(yold, acold, shgl, shp, 532 & iper, ilwork, rowp, colm, 533 & nsons, ifath, x, 534 & iBC, BC) 535 536 537 endif 538 539c.... set traction BCs for modeled walls 540c 541 if (itwmod.ne.0) then 542 call asbwmod(yold, acold, x, BC, iBC, 543 & iper, ilwork, ifath, velbar) 544 endif 545 546!MR CHANGE 547c 548c.... Determine whether the vorticity field needs to be computed for this time step or not 549c 550 icomputevort = 0 551 if (ivort == 1) then ! Print vorticity = True in solver.inp 552 ! We then compute the vorticity only if we 553 ! 1) we write an intermediate checkpoint 554 ! 2) we reach the last time step and write the last checkpoint 555 ! 3) we accumulate statistics in ybar for every time step 556 ! BEWARE: we need here lstep+1 and istep+1 because the lstep and 557 ! istep gets incremened after the flowsolve, further below 558 if (((irs .ge. 1) .and. (mod(lstep+1, ntout) .eq. 0)) .or. 559 & istep+1.eq.nstep(itseq) .or. ioybar == 1) then 560 icomputevort = 1 561 endif 562 endif 563 564! write(*,*) 'icomputevort: ',icomputevort, ' - istep: ', 565! & istep,' - nstep(itseq):',nstep(itseq),'- lstep:', 566! & lstep, '- ntout:', ntout 567!MR CHANGE 568 569c 570c.... -----------------------> predictor phase <----------------------- 571c 572 call itrPredict(yold, y, acold, ac , uold, u, iBC) 573 call itrBC (y, ac, iBC, BC, iper,ilwork) 574 575 if(nsolt.eq.1) then 576 isclr=0 577 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 578 endif 579 do isclr=1,nsclr 580 call itrBCSclr (y, ac, iBC, BC, iper, ilwork) 581 enddo 582 iter=0 583 ilss=0 ! this is a switch thrown on first solve of LS redistance 584 do istepc=1,seqsize 585 icode=stepseq(istepc) 586 if(mod(icode,5).eq.0) then ! this is a solve 587 isolve=icode/10 588 if(icode.eq.0) then ! flow solve (encoded as 0) 589c 590 iter = iter+1 591 ifuncs(1) = ifuncs(1) + 1 592c 593 Force(1) = zero 594 Force(2) = zero 595 Force(3) = zero 596 HFlux = zero 597 lhs = 1 - min(1,mod(ifuncs(1)-1,LHSupd(1))) 598 599 call SolFlow(y, ac, u, 600 & yold, acold, uold, 601 & x, iBC, 602 & BC, res, 603 & nPermDims, nTmpDims, aperm, 604 & atemp, iper, 605 & ilwork, shp, shgl, 606 & shpb, shglb, rowp, 607 & colm, lhsK, lhsP, 608 & solinc, rerr, tcorecp, 609 & GradV) 610 611 else ! scalar type solve 612 if (icode.eq.5) then ! Solve for Temperature 613 ! (encoded as (nsclr+1)*10) 614 isclr=0 615 ifuncs(2) = ifuncs(2) + 1 616 j=1 617 else ! solve a scalar (encoded at isclr*10) 618 isclr=isolve 619 ifuncs(isclr+2) = ifuncs(isclr+2) + 1 620 j=isclr+nsolt 621 if((iLSet.eq.2).and.(ilss.eq.0) 622 & .and.(isclr.eq.2)) then 623 ilss=1 ! throw switch (once per step) 624 y(:,7)=y(:,6) ! redistance field initialized 625 ac(:,7) = zero 626 call itrBCSclr ( y, ac, iBC, BC, iper, 627 & ilwork) 628c 629c....store the flow alpha, gamma parameter values and assigm them the 630c....Backward Euler parameters to solve the second levelset scalar 631c 632 alfit=alfi 633 gamit=gami 634 almit=almi 635 Deltt=Delt(1) 636 Dtglt=Dtgl 637 alfi = 1 638 gami = 1 639 almi = 1 640c Delt(1)= Deltt ! Give a pseudo time step 641 Dtgl = one / Delt(1) 642 endif ! level set eq. 2 643 endif ! deciding between temperature and scalar 644 645 lhs = 1 - min(1,mod(ifuncs(isclr+2)-1, 646 & LHSupd(isclr+2))) 647 648 call SolSclr(y, ac, u, 649 & yold, acold, uold, 650 & x, iBC, 651 & BC, nPermDimsS,nTmpDimsS, 652 & apermS(1,1,j), atempS, iper, 653 & ilwork, shp, shgl, 654 & shpb, shglb, rowp, 655 & colm, lhsS(1,j), 656 & solinc(1,isclr+5), tcorecpscal) 657 658 659 endif ! end of scalar type solve 660 661 else ! this is an update (mod did not equal zero) 662 iupdate=icode/10 ! what to update 663 if(icode.eq.1) then !update flow 664 call itrCorrect ( y, ac, u, solinc, iBC) 665 call itrBC (y, ac, iBC, BC, iper, ilwork) 666 else ! update scalar 667 isclr=iupdate !unless 668 if(icode.eq.6) isclr=0 669 if(iRANS.lt.-100)then ! RANS 670 call itrCorrectSclrPos(y,ac,solinc(1,isclr+5)) 671 else 672 call itrCorrectSclr (y, ac, solinc(1,isclr+5)) 673 endif 674 if (ilset.eq.2 .and. isclr.eq.2) then 675 if (ivconstraint .eq. 1) then 676 call itrBCSclr ( y, ac, iBC, BC, iper, 677 & ilwork) 678c 679c ... applying the volume constraint on second level set scalar 680c 681 call solvecon (y, x, iBC, BC, 682 & iper, ilwork, shp, shgl) 683c 684 endif ! end of volume constraint calculations 685 endif ! end of redistance calculations 686c 687 call itrBCSclr ( y, ac, iBC, BC, iper, 688 & ilwork) 689 endif ! end of flow or scalar update 690 endif ! end of switch between solve or update 691 enddo ! loop over sequence in step 692c 693c 694c.... obtain the time average statistics 695c 696 if (ioform .eq. 2) then 697 698 call stsGetStats( y, yold, ac, acold, 699 & u, uold, x, 700 & shp, shgl, shpb, shglb, 701 & iBC, BC, iper, ilwork, 702 & rowp, colm, lhsK, lhsP ) 703 endif 704 705c 706c Find the solution at the end of the timestep and move it to old 707c 708c 709c ...First to reassign the parameters for the original time integrator scheme 710c 711 if((iLSet.eq.2).and.(ilss.eq.1)) then 712 alfi =alfit 713 gami =gamit 714 almi =almit 715 Delt(1)=Deltt 716 Dtgl =Dtglt 717 endif 718 call itrUpdate( yold, acold, uold, y, ac, u) 719 call itrBC (yold, acold, iBC, BC, iper,ilwork) 720 721 istep = istep + 1 722 lstep = lstep + 1 723c 724c .. Print memory consumption on BGQ 725c 726 call printmeminfo("itrdrv"//char(0)) 727 728c 729c .. Compute vorticity 730c 731 if ( icomputevort == 1) then 732 733 ! vorticity components and magnitude 734 vorticity(:,1) = GradV(:,8)-GradV(:,6) !omega_x 735 vorticity(:,2) = GradV(:,3)-GradV(:,7) !omega_y 736 vorticity(:,3) = GradV(:,4)-GradV(:,2) !omega_z 737 vorticity(:,4) = sqrt( vorticity(:,1)*vorticity(:,1) 738 & + vorticity(:,2)*vorticity(:,2) 739 & + vorticity(:,3)*vorticity(:,3) ) 740 ! Q 741 strain(:,1) = GradV(:,1) !S11 742 strain(:,2) = 0.5*(GradV(:,2)+GradV(:,4)) !S12 743 strain(:,3) = 0.5*(GradV(:,3)+GradV(:,7)) !S13 744 strain(:,4) = GradV(:,5) !S22 745 strain(:,5) = 0.5*(GradV(:,6)+GradV(:,8)) !S23 746 strain(:,6) = GradV(:,9) !S33 747 748 vorticity(:,5) = 0.25*( vorticity(:,4)*vorticity(:,4) !Q 749 & - 2.0*( strain(:,1)*strain(:,1) 750 & + 2* strain(:,2)*strain(:,2) 751 & + 2* strain(:,3)*strain(:,3) 752 & + strain(:,4)*strain(:,4) 753 & + 2* strain(:,5)*strain(:,5) 754 & + strain(:,6)*strain(:,6))) 755 756 endif 757c 758c.... update and the aerodynamic forces 759c 760 call forces ( yold, ilwork ) 761 762c 763c .. write out the instantaneous solution 764c 7652001 continue ! we could get here by 2001 label if user requested stop 766 if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or. 767 & istep.eq.nstep(itseq)) then 768 769!so that we can see progress in force file close it so that it flushes 770!and then reopen in append mode 771 772 close(iforce) 773 open (unit=iforce, file=fforce, position='append') 774 775! Call to restar() will open restart file in write mode (and not append mode) 776! that is needed as other fields are written in append mode 777 778 call restar ('out ', yold ,ac) 779 if(ideformwall == 1) then 780! call write_displ(myrank, lstep, nshg, 3, uold ) 781 if(myrank.eq.master) then 782 write(*,*) 'ideformwall is 1 and is a dead code path... exiting' 783 endif 784 stop 785 endif 786 787 if(ivort == 1) then 788 call write_field(myrank,'a','vorticity',9,vorticity, 789 & 'd',nshg,5,lstep) 790 endif 791 792 call printmeminfo("itrdrv after checkpoint"//char(0)) 793 else if(stopjob.eq.-2) then 794 if(myrank.eq.master) then 795 write(*,*) 'line 755 says no write before stopping' 796 write(*,*) 'istep,nstep,irs',istep,nstep(itseq),irs 797 endif 798 endif !just the instantaneous stuff for videos 799c 800c.... compute the consistent boundary flux 801c 802 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 803 call Bflux ( yold, acold, uold, x, 804 & shp, shgl, shpb, 805 & shglb, ilwork, iBC, 806 & BC, iper, wallssVec) 807 endif 808 809 if(stopjob.eq.-2) goto 2003 810 811 812c 813c ... update the flow history for the impedance convolution, filter it and write it out 814c 815 if(numImpSrfs.gt.zero) then 816 call UpdHistConv(y,nsrflistImp,numImpSrfs) !uses Delt(1) 817 endif 818 819c 820c ... update the flow history for the RCR convolution 821c 822 if(numRCRSrfs.gt.zero) then 823 call UpdHistConv(y,nsrflistRCR,numRCRSrfs) !uses lstep 824 endif 825 826 827c... dump TIME SERIES 828 829 if (exts) then 830 if (mod(lstep-1,freq).eq.0) then 831 832 if (numpe > 1) then 833 do jj = 1, ntspts 834 vartssoln((jj-1)*ndof+1:jj*ndof)=varts(jj,:) 835 ivarts=zero 836 enddo 837 do k=1,ndof*ntspts 838 if(vartssoln(k).ne.zero) ivarts(k)=1 839 enddo 840 841! call MPI_REDUCE(vartssoln, vartssolng, ndof*ntspts, 842! & MPI_DOUBLE_PRECISION, MPI_SUM, master, 843! & MPI_COMM_WORLD, ierr) 844 845 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 846 call MPI_ALLREDUCE(vartssoln, vartssolng, 847 & ndof*ntspts, 848 & MPI_DOUBLE_PRECISION, MPI_SUM, 849 & MPI_COMM_WORLD, ierr) 850 851! call MPI_REDUCE(ivarts, ivartsg, ndof*ntspts, 852! & MPI_INTEGER, MPI_SUM, master, 853! & MPI_COMM_WORLD, ierr) 854 855 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 856 call MPI_ALLREDUCE(ivarts, ivartsg, ndof*ntspts, 857 & MPI_INTEGER, MPI_SUM, 858 & MPI_COMM_WORLD, ierr) 859 860! if (myrank.eq.zero) then 861 do jj = 1, ntspts 862 863 if(myrank .eq. iv_rank(jj)) then 864 ! No need to update all varts components, only the one treated by the expected rank 865 ! Note: keep varts as a vector, as multiple probes could be treated by one rank 866 indxvarts = (jj-1)*ndof 867 do k=1,ndof 868 if(ivartsg(indxvarts+k).ne.0) then ! none of the vartssoln(parts) were non zero 869 varts(jj,k)=vartssolng(indxvarts+k)/ 870 & ivartsg(indxvarts+k) 871 endif 872 enddo 873 endif !only if myrank eq iv_rank(jj) 874 enddo 875! endif !only on master 876 endif !only if numpe > 1 877 878! if (myrank.eq.zero) then 879 do jj = 1, ntspts 880 if(myrank .eq. iv_rank(jj)) then 881 ifile = 1000+jj 882 write(ifile,555) lstep, (varts(jj,k),k=1,ndof) !Beware of format 555 - check ndof 883c call flush(ifile) 884 if (((irs .ge. 1) .and. 885 & (mod(lstep, ntout) .eq. 0))) then 886 close(ifile) 887 fvarts='varts/varts' 888 fvarts=trim(fvarts)//trim(cname2(jj)) 889 fvarts=trim(fvarts)//trim(cname2(lskeep)) 890 fvarts=trim(fvarts)//'.dat' 891 fvarts=trim(fvarts) 892 open(unit=ifile, file=fvarts, 893 & position='append') 894 endif !only when dumping restart 895 endif 896 enddo 897! endif !only on master 898 899 varts(:,:) = zero ! reset the array for next step 900 901 902!555 format(i6,5(2x,E12.5e2)) 903555 format(i6,6(2x,E20.12e2)) !assuming ndof = 6 here 904 905 endif 906 endif 907 908 if((irscale.ge.0).or.(itwmod.gt.0)) 909 & call getvel (yold, ilwork, iBC, 910 & nsons, ifath, velbar) 911 912 if((irscale.ge.0).and.(myrank.eq.master)) then 913 call genscale(yold, x, iper, 914 & iBC, ifath, velbar, 915 & nsons) 916 endif 917c 918c.... print out results. 919c 920 ntoutv=max(ntout,100) ! velb is not needed so often 921 if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then 922 if( (mod(lstep, ntoutv) .eq. 0) .and. 923 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 924 & ((nsonmax.eq.1).and.(iLES.gt.0)))) 925 & call rwvelb ('out ', velbar ,ifail) 926 endif 927c 928c.... end of the NSTEP and NTSEQ loops 929c 930 931 932c 933c.... -------------------> error calculation <----------------- 934c 935 if(ierrcalc.eq.1 .or. ioybar.eq.1) then 936c$$$c 937c$$$c compute average 938c$$$c 939c$$$ tfact=one/istep 940c$$$ ybar =tfact*yold + (one-tfact)*ybar 941 942c compute average 943c ybar(:,1:3) are average velocity components 944c ybar(:,4) is average pressure 945c ybar(:,5) is average speed 946c ybar(:,6:8) is average of sq. of each vel. component 947c ybar(:,9) is average of sq. of pressure 948c ybar(:,10:12) is average of cross vel. components : uv, uw and vw 949c averaging procedure justified only for identical time step sizes 950c ybar(:,13) is average of eddy viscosity 951c ybar(:,14:16) is average vorticity components 952c ybar(:,17) is average vorticity magnitude 953c istep is number of time step 954c 955 icollectybar = 0 956 if(nphasesincycle.eq.0 .or. 957 & istep.gt.ncycles_startphaseavg*nstepsincycle) then 958 icollectybar = 1 959 if((istep-1).eq.ncycles_startphaseavg*nstepsincycle) 960 & istepsinybar = 0 ! init. to zero in first cycle in avg. 961 endif 962 963 if(icollectybar.eq.1) then 964 istepsinybar = istepsinybar+1 965 tfact=one/istepsinybar 966 967 if(myrank.eq.master .and. nphasesincycle.ne.0 .and. 968 & mod((istep-1),nstepsincycle).eq.0) 969 & write(*,*)'nsamples in phase average:',istepsinybar 970 971c ybar to contain the averaged ((u,v,w),p)-fields 972c and speed average, i.e., sqrt(u^2+v^2+w^2) 973c and avg. of sq. terms including 974c u^2, v^2, w^2, p^2 and cross terms of uv, uw and vw 975 976 ybar(:,1) = tfact*yold(:,1) + (one-tfact)*ybar(:,1) 977 ybar(:,2) = tfact*yold(:,2) + (one-tfact)*ybar(:,2) 978 ybar(:,3) = tfact*yold(:,3) + (one-tfact)*ybar(:,3) 979 ybar(:,4) = tfact*yold(:,4) + (one-tfact)*ybar(:,4) 980 ybar(:,5) = tfact*sqrt(yold(:,1)**2+yold(:,2)**2+ 981 & yold(:,3)**2) + (one-tfact)*ybar(:,5) 982 ybar(:,6) = tfact*yold(:,1)**2 + 983 & (one-tfact)*ybar(:,6) 984 ybar(:,7) = tfact*yold(:,2)**2 + 985 & (one-tfact)*ybar(:,7) 986 ybar(:,8) = tfact*yold(:,3)**2 + 987 & (one-tfact)*ybar(:,8) 988 ybar(:,9) = tfact*yold(:,4)**2 + 989 & (one-tfact)*ybar(:,9) 990 ybar(:,10) = tfact*yold(:,1)*yold(:,2) + !uv 991 & (one-tfact)*ybar(:,10) 992 ybar(:,11) = tfact*yold(:,1)*yold(:,3) + !uw 993 & (one-tfact)*ybar(:,11) 994 ybar(:,12) = tfact*yold(:,2)*yold(:,3) + !vw 995 & (one-tfact)*ybar(:,12) 996!MR CHANGE 997 if(nsclr.gt.0) !nut 998 & ybar(:,13) = tfact*yold(:,6) + (one-tfact)*ybar(:,13) 999 1000 if(ivort == 1) then !vorticity 1001 ybar(:,14) = tfact*vorticity(:,1) + 1002 & (one-tfact)*ybar(:,14) 1003 ybar(:,15) = tfact*vorticity(:,2) + 1004 & (one-tfact)*ybar(:,15) 1005 ybar(:,16) = tfact*vorticity(:,3) + 1006 & (one-tfact)*ybar(:,16) 1007 ybar(:,17) = tfact*vorticity(:,4) + 1008 & (one-tfact)*ybar(:,17) 1009 endif 1010 1011 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1012 wallssVecBar(:,1) = tfact*wallssVec(:,1) 1013 & +(one-tfact)*wallssVecBar(:,1) 1014 wallssVecBar(:,2) = tfact*wallssVec(:,2) 1015 & +(one-tfact)*wallssVecBar(:,2) 1016 wallssVecBar(:,3) = tfact*wallssVec(:,3) 1017 & +(one-tfact)*wallssVecBar(:,3) 1018 endif 1019!MR CHANGE END 1020 endif 1021c 1022c compute phase average 1023c 1024 if(nphasesincycle.ne.0 .and. 1025 & istep.gt.ncycles_startphaseavg*nstepsincycle) then 1026 1027c beginning of cycle is considered as ncycles_startphaseavg*nstepsincycle+1 1028 if((istep-1).eq.ncycles_startphaseavg*nstepsincycle) 1029 & icyclesinavg = 0 ! init. to zero in first cycle in avg. 1030 1031 ! find number of steps between phases 1032 nstepsbtwphase = nstepsincycle/nphasesincycle ! integer value 1033 if(mod(istep-1,nstepsincycle).eq.0) then 1034 iphase = 1 ! init. to one in beginning of every cycle 1035 icyclesinavg = icyclesinavg + 1 1036 endif 1037 1038 icollectphase = 0 1039 istepincycle = mod(istep,nstepsincycle) 1040 if(istepincycle.eq.0) istepincycle=nstepsincycle 1041 if(istepincycle.eq.iphase*nstepsbtwphase) then 1042 icollectphase = 1 1043 iphase = iphase+1 ! use 'iphase-1' below 1044 endif 1045 1046 if(icollectphase.eq.1) then 1047 tfactphase = one/icyclesinavg 1048 1049 if(myrank.eq.master) then 1050 write(*,*) 'nsamples in phase ',iphase-1,': ', 1051 & icyclesinavg 1052 endif 1053 1054 yphbar(:,1,iphase-1) = tfactphase*yold(:,1) + 1055 & (one-tfactphase)*yphbar(:,1,iphase-1) 1056 yphbar(:,2,iphase-1) = tfactphase*yold(:,2) + 1057 & (one-tfactphase)*yphbar(:,2,iphase-1) 1058 yphbar(:,3,iphase-1) = tfactphase*yold(:,3) + 1059 & (one-tfactphase)*yphbar(:,3,iphase-1) 1060 yphbar(:,4,iphase-1) = tfactphase*yold(:,4) + 1061 & (one-tfactphase)*yphbar(:,4,iphase-1) 1062 yphbar(:,5,iphase-1) = tfactphase*sqrt(yold(:,1)**2 1063 & +yold(:,2)**2+yold(:,3)**2) + 1064 & (one-tfactphase)*yphbar(:,5,iphase-1) 1065!MR CHANGE 1066 yphbar(:,6,iphase-1) = 1067 & tfactphase*yold(:,1)*yold(:,1) 1068 & +(one-tfactphase)*yphbar(:,6,iphase-1) 1069 1070 yphbar(:,7,iphase-1) = 1071 & tfactphase*yold(:,1)*yold(:,2) 1072 & +(one-tfactphase)*yphbar(:,7,iphase-1) 1073 1074 yphbar(:,8,iphase-1) = 1075 & tfactphase*yold(:,1)*yold(:,3) 1076 & +(one-tfactphase)*yphbar(:,8,iphase-1) 1077 1078 yphbar(:,9,iphase-1) = 1079 & tfactphase*yold(:,2)*yold(:,2) 1080 & +(one-tfactphase)*yphbar(:,9,iphase-1) 1081 1082 yphbar(:,10,iphase-1) = 1083 & tfactphase*yold(:,2)*yold(:,3) 1084 & +(one-tfactphase)*yphbar(:,10,iphase-1) 1085 1086 yphbar(:,11,iphase-1) = 1087 & tfactphase*yold(:,3)*yold(:,3) 1088 & +(one-tfactphase)*yphbar(:,11,iphase-1) 1089 1090 if(ivort == 1) then 1091 yphbar(:,12,iphase-1) = 1092 & tfactphase*vorticity(:,1) 1093 & +(one-tfactphase)*yphbar(:,12,iphase-1) 1094 yphbar(:,13,iphase-1) = 1095 & tfactphase*vorticity(:,2) 1096 & +(one-tfactphase)*yphbar(:,13,iphase-1) 1097 yphbar(:,14,iphase-1) = 1098 & tfactphase*vorticity(:,3) 1099 & +(one-tfactphase)*yphbar(:,14,iphase-1) 1100 yphbar(:,15,iphase-1) = 1101 & tfactphase*vorticity(:,4) 1102 & +(one-tfactphase)*yphbar(:,15,iphase-1) 1103 endif 1104!MR CHANGE END 1105 endif 1106 endif 1107c 1108c compute rms 1109c 1110 if(icollectybar.eq.1) then 1111 rerr(:, 7)=rerr(:, 7)+(yold(:,1)-ybar(:,1))**2 1112 rerr(:, 8)=rerr(:, 8)+(yold(:,2)-ybar(:,2))**2 1113 rerr(:, 9)=rerr(:, 9)+(yold(:,3)-ybar(:,3))**2 1114 rerr(:,10)=rerr(:,10)+(yold(:,4)-ybar(:,4))**2 1115 endif 1116 endif 1117 2003 continue ! we get here if stopjob equals lstep and this jumped over 1118! the statistics computation because we have no new data to average in 1119! rather we are just trying to output the last state that was not already 1120! written 1121c 1122c.... ----------------------> Complete Restart Processing <---------------------- 1123c 1124! for now it is the same frequency but need to change this 1125! soon.... but don't forget to change the field counter in 1126! new_interface.cc 1127! 1128 if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or. 1129 & istep.eq.nstep(itseq)) then 1130 if ((irs .ge. 1) .and. ((mod(lstep, ntout) .ne. 0) .or. 1131 & (nstp .eq. 0))) then 1132 if( 1133 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 1134 & ((nsonmax.eq.1).and.iLES.gt.0))) 1135 & call rwvelb ('out ', velbar ,ifail) 1136 endif 1137 1138 lesId = numeqns(1) 1139 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1140 if(myrank.eq.0) then 1141 tcormr1 = TMRC() 1142 endif 1143 call saveLesRestart( lesId, aperm , nshg, myrank, lstep, 1144 & nPermDims ) 1145 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1146 if(myrank.eq.0) then 1147 tcormr2 = TMRC() 1148 write(6,*) 'call saveLesRestart for projection and'// 1149 & 'pressure projection vectors', tcormr2-tcormr1 1150 endif 1151 1152 if(ierrcalc.eq.1) then 1153c 1154c.....smooth the error indicators 1155c 1156 do i=1,ierrsmooth 1157 call errsmooth( rerr, x, iper, ilwork, shp, shgl, iBC ) 1158 end do 1159 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1160 if(myrank.eq.0) then 1161 tcormr1 = TMRC() 1162 endif 1163 call write_error(myrank, lstep, nshg, 10, rerr ) 1164 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1165 if(myrank.eq.0) then 1166 tcormr2 = TMRC() 1167 write(6,*) 'Time to write the error fields to the disks', 1168 & tcormr2-tcormr1 1169 endif 1170 endif ! ierrcalc 1171 1172 if(ioybar.eq.1) then 1173 if(ivort == 1) then 1174 call write_field(myrank,'a','ybar',4, 1175 & ybar,'d',nshg,17,lstep) 1176 else 1177 call write_field(myrank,'a','ybar',4, 1178 & ybar,'d',nshg,13,lstep) 1179 endif 1180 1181 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1182 call write_field(myrank,'a','wssbar',6, 1183 & wallssVecBar,'d',nshg,3,lstep) 1184 endif 1185 1186 if(nphasesincycle .gt. 0) then 1187 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1188 if(myrank.eq.0) then 1189 tcormr1 = TMRC() 1190 endif 1191 do iphase=1,nphasesincycle 1192 if(ivort == 1) then 1193 call write_phavg2(myrank,'a','phase_average',13,iphase, 1194 & nphasesincycle,yphbar(:,:,iphase),'d',nshg,15,lstep) 1195 else 1196 call write_phavg2(myrank,'a','phase_average',13,iphase, 1197 & nphasesincycle,yphbar(:,:,iphase),'d',nshg,11,lstep) 1198 endif 1199 end do 1200 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1201 if(myrank.eq.0) then 1202 tcormr2 = TMRC() 1203 write(6,*) 'write all phase avg to the disks = ', 1204 & tcormr2-tcormr1 1205 endif 1206 endif !nphasesincyle 1207 endif !ioybar 1208 1209 if ( ( ihessian .eq. 1 ) .and. ( numpe < 2 ) )then 1210 uhess = zero 1211 gradu = zero 1212 tf = zero 1213 1214 do ku=1,nshg 1215 tf(ku,1) = x(ku,1)**3 1216 end do 1217 call hessian( yold, x, shp, shgl, iBC, 1218 & shpb, shglb, iper, ilwork, uhess, gradu ) 1219 1220 call write_hessian( uhess, gradu, nshg ) 1221 endif 1222 1223 if(iRANS.lt.0) then 1224 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1225 if(myrank.eq.0) then 1226 tcormr1 = TMRC() 1227 endif 1228 call write_field(myrank,'a','dwal',4,d2wall,'d', 1229 & nshg,1,lstep) 1230 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1231 if(myrank.eq.0) then 1232 tcormr2 = TMRC() 1233 write(6,*) 'Time to write dwal to the disks = ', 1234 & tcormr2-tcormr1 1235 endif 1236 endif !iRANS 1237 1238 endif ! write out complete restart state 1239 !next 2 lines are two ways to end early 1240 if(stopjob.eq.-2) goto 2002 1241 if(istop.eq.1000) goto 2002 ! stop when delta small (see rstatic) 1242 2000 continue 1243 2002 continue 1244 1245! done with time stepping so deallocate fields already written 1246! 1247 if(ioybar.eq.1) then 1248 deallocate(ybar) 1249 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1250 deallocate(wallssVecbar) 1251 endif 1252 if(nphasesincycle .gt. 0) then 1253 deallocate(yphbar) 1254 endif !nphasesincyle 1255 endif !ioybar 1256 if(ivort == 1) then 1257 deallocate(strain,vorticity) 1258 endif 1259 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1260 deallocate(wallssVec) 1261 endif 1262 if(iRANS.lt.0) then 1263 deallocate(d2wall) 1264 endif 1265 1266 1267CAD tcorecp2 = second(0) 1268CAD tcorewc2 = second(-1) 1269 1270CAD write(6,*) 'T(core) cpu-wallclock = ',tcorecp2-tcorecp1, 1271CAD & tcorewc2-tcorewc1 1272 1273 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1274 if(myrank.eq.0) then 1275 tcorecp2 = TMRC() 1276 write(6,*) 'T(core) cpu = ',tcorecp2-tcorecp1 1277 write(6,*) '(Elm. form.',tcorecp(1), 1278 & ',Lin. alg. sol.',tcorecp(2),')' 1279 write(6,*) '(Elm. form. Scal.',tcorecpscal(1), 1280 & ',Lin. alg. sol. Scal.',tcorecpscal(2),')' 1281 write(6,*) '' 1282 1283 endif 1284 1285 call print_system_stats(tcorecp, tcorecpscal) 1286 call print_mesh_stats() 1287 call print_mpi_stats() 1288 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1289! return 1290c call MPI_Finalize() 1291c call MPI_ABORT(MPI_COMM_WORLD, ierr) 1292 1293 3000 continue 1294 1295 1296c 1297c.... close history and aerodynamic forces files 1298c 1299 if (myrank .eq. master) then 1300! close (ihist) 1301 close (iforce) 1302 close(76) 1303 if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then 1304 close (8177) 1305 endif 1306 endif 1307c 1308c.... close varts file for probes 1309c 1310 if(exts) then 1311 do jj=1,ntspts 1312 if (myrank == iv_rank(jj)) then 1313 close(1000+jj) 1314 endif 1315 enddo 1316 deallocate (ivarts) 1317 deallocate (ivartsg) 1318 deallocate (iv_rank) 1319 deallocate (vartssoln) 1320 deallocate (vartssolng) 1321 endif 1322 1323!MR CHANGE 1324 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1325 if(myrank.eq.0) then 1326 write(*,*) 'itrdrv - done with aerodynamic forces' 1327 endif 1328!MR CHANGE 1329 1330 do isrf = 0,MAXSURF 1331! if ( nsrflist(isrf).ne.0 ) then 1332 if ( nsrflist(isrf).ne.0 .and. 1333 & myrank.eq.irankfilesforce(isrf)) then 1334 iunit=60+isrf 1335 close(iunit) 1336 endif 1337 enddo 1338 1339!MR CHANGE 1340 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1341 if(myrank.eq.0) then 1342 write(*,*) 'itrdrv - done with MAXSURF' 1343 endif 1344!MR CHANGE 1345 1346 1347 5 format(1X,F15.10,3X,F15.10,3X,F15.10,3X,F15.10) 1348 444 format(6(2x,e14.7)) 1349c 1350c.... end 1351c 1352 if(nsolflow.eq.1) then 1353 deallocate (lhsK) 1354 deallocate (lhsP) 1355 deallocate (aperm) 1356 deallocate (atemp) 1357 endif 1358 if(nsclrsol.gt.0) then 1359 deallocate (lhsS) 1360 deallocate (apermS) 1361 deallocate (atempS) 1362 endif 1363 1364 if(iabc==1) deallocate(acs) 1365 1366!MR CHANGE 1367 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1368 if(myrank.eq.0) then 1369 write(*,*) 'itrdrv - done - BACK TO process.f' 1370 endif 1371!MR CHANGE 1372 1373 1374 1375 return 1376 end 1377 1378 subroutine lesmodels(y, ac, shgl, shp, 1379 & iper, ilwork, rowp, colm, 1380 & nsons, ifath, x, 1381 & iBC, BC) 1382 1383 include "common.h" 1384 1385 real*8 y(nshg,ndof), ac(nshg,ndof), 1386 & x(numnp,nsd), 1387 & BC(nshg,ndofBC) 1388 real*8 shp(MAXTOP,maxsh,MAXQPT), 1389 & shgl(MAXTOP,nsd,maxsh,MAXQPT) 1390 1391c 1392 integer rowp(nshg,nnz), colm(nshg+1), 1393 & iBC(nshg), 1394 & ilwork(nlwork), 1395 & iper(nshg) 1396 dimension ifath(numnp), nsons(nfath) 1397 1398 real*8, allocatable, dimension(:) :: fwr2,fwr3,fwr4 1399 real*8, allocatable, dimension(:) :: stabdis,cdelsq1 1400 real*8, allocatable, dimension(:,:) :: xavegt, xavegt2,xavegt3 1401 1402 if( (iLES.gt.1) ) then ! Allocate Stuff for advanced LES models 1403 allocate (fwr2(nshg)) 1404 allocate (fwr3(nshg)) 1405 allocate (fwr4(nshg)) 1406 allocate (xavegt(nfath,12)) 1407 allocate (xavegt2(nfath,12)) 1408 allocate (xavegt3(nfath,12)) 1409 allocate (stabdis(nfath)) 1410 endif 1411 1412c.... get dynamic model coefficient 1413c 1414 ilesmod=iLES/10 1415c 1416c digit bit set filter rule, 10 bit set model 1417c 1418 if (ilesmod.eq.0) then ! 0 < iLES< 10 => dyn. model calculated 1419 ! at nodes based on discrete filtering 1420 1421 1422 if(isubmod.eq.2) then 1423 call SUPGdis(y, ac, shgl, shp, 1424 & iper, ilwork, 1425 & nsons, ifath, x, 1426 & iBC, BC, stabdis, xavegt3) 1427 endif 1428 1429 if( ((isubmod.eq.0).or.(isubmod.eq.2)))then ! If no 1430 ! sub-model 1431 ! or SUPG 1432 ! model wanted 1433 1434 if(i2filt.eq.0)then ! If simple filter 1435 1436 if(modlstats .eq. 0) then ! If no model stats wanted 1437 call getdmc (y, shgl, shp, 1438 & iper, ilwork, nsons, 1439 & ifath, x) 1440 else ! else get model stats 1441 call stdfdmc (y, shgl, shp, 1442 & iper, ilwork, nsons, 1443 & ifath, x) 1444 endif ! end of stats if statement 1445 1446 else ! else if twice filtering 1447 1448 call widefdmc(y, shgl, shp, 1449 & iper, ilwork, nsons, 1450 & ifath, x) 1451 1452 1453 endif ! end of simple filter if statement 1454 1455 endif ! end of SUPG or no sub-model if statement 1456 1457 1458 if( (isubmod.eq.1) ) then ! If DFWR sub-model wanted 1459 call cdelBHsq (y, shgl, shp, 1460 & iper, ilwork, nsons, 1461 & ifath, x, cdelsq1) 1462 call FiltRat (y, shgl, shp, 1463 & iper, ilwork, nsons, 1464 & ifath, x, cdelsq1, 1465 & fwr4, fwr3) 1466 1467 1468 if (i2filt.eq.0) then ! If simple filter wanted 1469 call DFWRsfdmc(y, shgl, shp, 1470 & iper, ilwork, nsons, 1471 & ifath, x, fwr2, fwr3) 1472 else ! else if twice filtering wanted 1473 call DFWRwfdmc(y, shgl, shp, 1474 & iper, ilwork, nsons, 1475 & ifath, x, fwr4, fwr4) 1476 endif ! end of simple filter if statement 1477 1478 endif ! end of DFWR sub-model if statement 1479 1480 if( (isubmod.eq.2) )then ! If SUPG sub-model wanted 1481 call dmcSUPG (y, ac, shgl, 1482 & shp, iper, ilwork, 1483 & nsons, ifath, x, 1484 & iBC, BC, rowp, colm, 1485 & xavegt2, stabdis) 1486 endif 1487 1488 if(idis.eq.1)then ! If SUPG/Model dissipation wanted 1489 call ediss (y, ac, shgl, 1490 & shp, iper, ilwork, 1491 & nsons, ifath, x, 1492 & iBC, BC, xavegt) 1493 endif 1494 1495 endif ! end of ilesmod 1496 1497 if (ilesmod .eq. 1) then ! 10 < iLES < 20 => dynamic-mixed 1498 ! at nodes based on discrete filtering 1499 call bardmc (y, shgl, shp, 1500 & iper, ilwork, 1501 & nsons, ifath, x) 1502 endif 1503 1504 if (ilesmod .eq. 2) then ! 20 < iLES < 30 => dynamic at quad 1505 ! pts based on lumped projection filt. 1506 1507 if(isubmod.eq.0)then 1508 call projdmc (y, shgl, shp, 1509 & iper, ilwork, x) 1510 else 1511 call cpjdmcnoi (y, shgl, shp, 1512 & iper, ilwork, x, 1513 & rowp, colm, 1514 & iBC, BC) 1515 endif 1516 1517 endif 1518 1519 if( (iLES.gt.1) ) then ! Deallocate Stuff for advanced LES models 1520 deallocate (fwr2) 1521 deallocate (fwr3) 1522 deallocate (fwr4) 1523 deallocate (xavegt) 1524 deallocate (xavegt2) 1525 deallocate (xavegt3) 1526 deallocate (stabdis) 1527 endif 1528 return 1529 end 1530 1531c 1532c...initialize the coefficients for the impedance convolution 1533c 1534 subroutine CalcImpConvCoef (numISrfs, numTpoints) 1535 1536 use convolImpFlow !uses flow history and impedance for convolution 1537 1538 include "common.h" !for alfi 1539 1540 integer numISrfs, numTpoints 1541 1542 allocate (ConvCoef(numTpoints+2,3)) !same time discret. for all imp. BC 1543 do j=1,numTpoints+2 1544 ConvCoef(j,:)=0.5/numTpoints !dt/2 divided by period T=N*dt 1545 ConvCoef(j,1)=ConvCoef(j,1)*(1.0-alfi)*(1.0-alfi) 1546 ConvCoef(j,2)=ConvCoef(j,2)*(1.0+2*alfi*(1.0-alfi)) 1547 ConvCoef(j,3)=ConvCoef(j,3)*alfi*alfi 1548 enddo 1549 ConvCoef(1,2)=zero 1550 ConvCoef(1,3)=zero 1551 ConvCoef(2,3)=zero 1552 ConvCoef(numTpoints+1,1)=zero 1553 ConvCoef(numTpoints+2,2)=zero 1554 ConvCoef(numTpoints+2,1)=zero 1555c 1556c...calculate the coefficients for the impedance convolution 1557c 1558 allocate (ImpConvCoef(numTpoints+2,numISrfs)) 1559 1560c..coefficients below assume Q linear in time step, Z constant 1561c do j=3,numTpoints 1562c ImpConvCoef(j,:) = ValueListImp(j-1,:)*ConvCoef(j,3) 1563c & + ValueListImp(j,:)*ConvCoef(j,2) 1564c & + ValueListImp(j+1,:)*ConvCoef(j,1) 1565c enddo 1566c ImpConvCoef(1,:) = ValueListImp(2,:)*ConvCoef(1,1) 1567c ImpConvCoef(2,:) = ValueListImp(2,:)*ConvCoef(2,2) 1568c & + ValueListImp(3,:)*ConvCoef(2,1) 1569c ImpConvCoef(numTpoints+1,:) = 1570c & ValueListImp(numTpoints,:)*ConvCoef(numTpoints+1,3) 1571c & + ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+1,2) 1572c ImpConvCoef(numTpoints+2,:) = 1573c & ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+2,3) 1574 1575c..try easiest convolution Q and Z constant per time step 1576 do j=3,numTpoints+1 1577 ImpConvCoef(j,:) = ValueListImp(j-1,:)/numTpoints 1578 enddo 1579 ImpConvCoef(1,:) =zero 1580 ImpConvCoef(2,:) =zero 1581 ImpConvCoef(numTpoints+2,:) = 1582 & ValueListImp(numTpoints+1,:)/numTpoints 1583c compensate for yalpha passed not y in Elmgmr() 1584 ImpConvCoef(numTpoints+1,:)= ImpConvCoef(numTpoints+1,:) 1585 & - ImpConvCoef(numTpoints+2,:)*(1.0-alfi)/alfi 1586 ImpConvCoef(numTpoints+2,:)= ImpConvCoef(numTpoints+2,:)/alfi 1587 return 1588 end 1589 1590c 1591c ... update the flow rate history for the impedance convolution, filter it and write it out 1592c 1593 subroutine UpdHistConv(y,nsrfIdList,numSrfs) 1594 1595 use convolImpFlow !brings ntimeptpT, QHistImp, QHistTry, QHistTryF, numImpSrfs 1596 use convolRCRFlow !brings QHistRCR, numRCRSrfs 1597 1598 include "common.h" 1599 1600 integer nsrfIdList(0:MAXSURF), numSrfs 1601 character*20 fname1 1602 character*10 cname2 1603 character*5 cname 1604 real*8 y(nshg,3) !velocity at time n+1 1605 real*8 NewQ(0:MAXSURF) !temporary unknown for the flow rate 1606 !that needs to be added to the flow history 1607 1608 call GetFlowQ(NewQ,y,nsrfIdList,numSrfs) !new flow at time n+1 1609c 1610c... for imp BC: shift QHist, add new constribution, filter and write out 1611c 1612 if(numImpSrfs.gt.zero) then 1613 do j=1, ntimeptpT 1614 QHistImp(j,1:numSrfs)=QHistImp(j+1,1:numSrfs) 1615 enddo 1616 QHistImp(ntimeptpT+1,1:numSrfs) = NewQ(1:numSrfs) 1617 1618c 1619c....filter the flow rate history 1620c 1621 cutfreq = 10 !hardcoded cutting frequency of the filter 1622 do j=1, ntimeptpT 1623 QHistTry(j,:)=QHistImp(j+1,:) 1624 enddo 1625 call Filter(QHistTryF,QHistTry,ntimeptpT,Delt(1),cutfreq) 1626c.... if no filter applied then uncomment next three lines 1627c do j=1, ntimeptpT 1628c QHistTryF(j,:)=QHistTry(j,:) 1629c enddo 1630 1631c QHistImp(1,:)=zero ! why do we do this? for beta(1,:) = zero it does not really matters 1632 do j=1, ntimeptpT 1633 QHistImp(j+1,:)=QHistTryF(j,:) 1634 enddo 1635c 1636c.... write out the new history of flow rates to Qhistor.dat 1637c 1638 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1639 & (istep .eq. nstep(1)))) .and. 1640 & (myrank .eq. master)) then 1641 open(unit=816, file='Qhistor.dat',status='replace') 1642 write(816,*) ntimeptpT 1643 do j=1,ntimeptpT+1 1644 write(816,*) (QHistImp(j,n),n=1, numSrfs) 1645 enddo 1646 close(816) 1647c... write out a copy with step number to be able to restart 1648 fname1 = 'Qhistor' 1649 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1650 open(unit=8166,file=trim(fname1),status='unknown') 1651 write(8166,*) ntimeptpT 1652 do j=1,ntimeptpT+1 1653 write(8166,*) (QHistImp(j,n),n=1, numSrfs) 1654 enddo 1655 close(8166) 1656 endif 1657 endif 1658 1659c 1660c... for RCR bc just add the new contribution 1661c 1662 if(numRCRSrfs.gt.zero) then 1663 QHistRCR(lstep+1,1:numSrfs) = NewQ(1:numSrfs) 1664c 1665c.... write out the new history of flow rates to Qhistor.dat 1666c 1667 if ((irs .ge. 1) .and. (myrank .eq. master)) then 1668 if(istep.eq.1) then 1669 open(unit=816,file='Qhistor.dat',status='unknown') 1670 else 1671 open(unit=816,file='Qhistor.dat',position='append') 1672 endif 1673 if(istep.eq.1) then 1674 do j=1,lstep 1675 write(816,*) j, (QHistRCR(j,n),n=1,numSrfs) ! read from file of previous run 1676 enddo 1677 endif 1678 write(816,*) lstep+1, (QHistRCR(lstep+1,n),n=1, numSrfs) 1679 close(816) 1680c... write out a copy with step number to be able to restart 1681 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1682 & (istep .eq. nstep(1)))) .and. 1683 & (myrank .eq. master)) then 1684 fname1 = 'Qhistor' 1685 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1686 open(unit=8166,file=trim(fname1),status='unknown') 1687 write(8166,*) lstep+1 1688 do j=1,lstep+1 1689 write(8166,*) (QHistRCR(j,n),n=1, numSrfs) 1690 enddo 1691 close(8166) 1692 endif 1693 endif 1694 endif 1695 1696 return 1697 end 1698 1699c 1700c...calculate the time varying coefficients for the RCR convolution 1701c 1702 subroutine CalcRCRConvCoef (stepn, numSrfs) 1703 1704 use convolRCRFlow !brings in ValueListRCR, dtRCR 1705 1706 include "common.h" !brings alfi 1707 1708 integer numSrfs, stepn 1709 1710 RCRConvCoef = zero 1711 if (stepn .eq. 0) then 1712 RCRConvCoef(1,:) = ValueListRCR(1,:)*(1.0-alfi) + 1713 & ValueListRCR(3,:)*(-alfi + 1.0 + 1/dtRCR(:) 1714 & - exp(-alfi*dtRCR(:))*(1 + 1/dtRCR(:))) 1715 RCRConvCoef(2,:) = ValueListRCR(1,:)*alfi 1716 & + ValueListRCR(3,:) 1717 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1718 endif 1719 if (stepn .ge. 1) then 1720 RCRConvCoef(1,:) =-ValueListRCR(3,:)*exp(-dtRCR(:)*(stepn+alfi)) 1721 & *(1 + (1 - exp(dtRCR(:)))/dtRCR(:)) 1722 RCRConvCoef(stepn+1,:) = ValueListRCR(1,:)*(1-alfi) 1723 & - ValueListRCR(3,:)*(alfi - 1 - 1/dtRCR(:) 1724 & + exp(-alfi*dtRCR(:))/dtRCR(:)*(2 - exp(-dtRCR(:)))) 1725 RCRConvCoef(stepn+2,:) = ValueListRCR(1,:)*alfi 1726 & + ValueListRCR(3,:) 1727 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1728 endif 1729 if (stepn .ge. 2) then 1730 do j=2,stepn 1731 RCRConvCoef(j,:) = ValueListRCR(3,:)/dtRCR(:)* 1732 & exp(-dtRCR(:)*(stepn + alfi + 2 - j))* 1733 & (1 - exp(dtRCR(:)))**2 1734 enddo 1735 endif 1736 1737c compensate for yalpha passed not y in Elmgmr() 1738 RCRConvCoef(stepn+1,:)= RCRConvCoef(stepn+1,:) 1739 & - RCRConvCoef(stepn+2,:)*(1.0-alfi)/alfi 1740 RCRConvCoef(stepn+2,:)= RCRConvCoef(stepn+2,:)/alfi 1741 1742 return 1743 end 1744 1745c 1746c...calculate the time dependent H operator for the RCR convolution 1747c 1748 subroutine CalcHopRCR (timestepRCR, stepn, numSrfs) 1749 1750 use convolRCRFlow !brings in HopRCR, dtRCR 1751 1752 include "common.h" 1753 1754 integer numSrfs, stepn 1755 real*8 PdistCur(0:MAXSURF), timestepRCR 1756 1757 HopRCR=zero 1758 call RCRint(timestepRCR*(stepn + alfi),PdistCur) 1759 HopRCR(1:numSrfs) = RCRic(1:numSrfs) 1760 & *exp(-dtRCR(1:numSrfs)*(stepn + alfi)) + PdistCur(1:numSrfs) 1761 return 1762 end 1763c 1764c ... initialize the influence of the initial conditions for the RCR BC 1765c 1766 subroutine calcRCRic(y,srfIdList,numSrfs) 1767 1768 use convolRCRFlow !brings RCRic, ValueListRCR, ValuePdist 1769 1770 include "common.h" 1771 1772 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled 1773 real*8 y(nshg,4) !need velocity and pressure 1774 real*8 Qini(0:MAXSURF) !initial flow rate 1775 real*8 PdistIni(0:MAXSURF) !initial distal pressure 1776 real*8 Pini(0:MAXSURF),CoupleArea(0:MAXSURF) ! initial pressure 1777 real*8 VelOnly(nshg,3), POnly(nshg) 1778 1779 allocate (RCRic(0:MAXSURF)) 1780 1781 if(lstep.eq.0) then 1782 VelOnly(:,1:3)=y(:,1:3) 1783 call GetFlowQ(Qini,VelOnly,srfIdList,numSrfs) !get initial flow 1784 QHistRCR(1,1:numSrfs)=Qini(1:numSrfs) !initialize QHistRCR 1785 POnly(:)=y(:,4) ! pressure 1786 call integrScalar(Pini,POnly,srfIdList,numSrfs) !get initial pressure integral 1787 POnly(:)=one ! one to get area 1788 call integrScalar(CoupleArea,POnly,srfIdList,numSrfs) !get surf area 1789 Pini(1:numSrfs) = Pini(1:numSrfs)/CoupleArea(1:numSrfs) 1790 else 1791 Qini(1:numSrfs)=QHistRCR(1,1:numSrfs) 1792 Pini(1:numSrfs)=zero ! hack 1793 endif 1794 call RCRint(istep,PdistIni) !get initial distal P (use istep) 1795 RCRic(1:numSrfs) = Pini(1:numSrfs) 1796 & - ValueListRCR(1,:)*Qini(1:numSrfs)-PdistIni(1:numSrfs) 1797 return 1798 end 1799 1800c.........function that integrates a scalar over a boundary 1801 subroutine integrScalar(scalInt,scal,srfIdList,numSrfs) 1802 1803 use pvsQbi !brings ndsurf, NASC 1804 1805 include "common.h" 1806 include "mpif.h" 1807 1808 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled, i, k 1809 real*8 scal(nshg), scalInt(0:MAXSURF), scalIntProc(0:MAXSURF) 1810 1811 scalIntProc = zero 1812 do i = 1,nshg 1813 if(numSrfs.gt.zero) then 1814 do k = 1,numSrfs 1815 irankCoupled = 0 1816 if (srfIdList(k).eq.ndsurf(i)) then 1817 irankCoupled=k 1818 scalIntProc(irankCoupled) = scalIntProc(irankCoupled) 1819 & + NASC(i)*scal(i) 1820 exit 1821 endif 1822 enddo 1823 endif 1824 enddo 1825c 1826c at this point, each scalint has its "nodes" contributions to the scalar 1827c accumulated into scalIntProc. Note, because NASC is on processor this 1828c will NOT be the scalar for the surface yet 1829c 1830c.... reduce integrated scalar for each surface, push on scalInt 1831c 1832 npars=MAXSURF+1 1833 call MPI_ALLREDUCE (scalIntProc, scalInt(:), npars, 1834 & MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr) 1835 1836 return 1837 end 1838 1839 subroutine writeTimingMessage(key,iomode,timing) 1840 use iso_c_binding 1841 use phstr 1842 implicit none 1843 1844 character(len=*) :: key 1845 integer :: iomode 1846 real*8 :: timing 1847 character(len=1024) :: timing_msg 1848 character(len=*), parameter :: 1849 & streamModeString = c_char_"stream"//c_null_char, 1850 & fileModeString = c_char_"disk"//c_null_char 1851 1852 timing_msg = c_char_"Time to write "//c_null_char 1853 call phstr_appendStr(timing_msg,key) 1854 if ( iomode .eq. -1 ) then 1855 call phstr_appendStr(timing_msg, streamModeString) 1856 else 1857 call phstr_appendStr(timing_msg, fileModeString) 1858 endif 1859 call phstr_appendStr(timing_msg, c_char_' = '//c_null_char) 1860 call phstr_appendDbl(timing_msg, timing) 1861 write(6,*) trim(timing_msg) 1862 return 1863 end subroutine 1864 1865