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 .. write out the instantaneous solution 759c 7602001 continue ! we could get here by 2001 label if user requested stop 761 if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or. 762 & istep.eq.nstep(itseq)) then 763 764!so that we can see progress in force file close it so that it flushes 765!and then reopen in append mode 766 767 close(iforce) 768 open (unit=iforce, file=fforce, position='append') 769 770! Call to restar() will open restart file in write mode (and not append mode) 771! that is needed as other fields are written in append mode 772 773 call restar ('out ', yold ,ac) 774 if(ideformwall == 1) then 775 if(myrank.eq.master) then 776 write(*,*) 'ideformwall is 1 and is a dead code path... exiting' 777 endif 778 stop 779 endif 780 781 if(ivort == 1) then 782 call write_field(myrank,'a','vorticity',9,vorticity, 783 & 'd',nshg,5,lstep) 784 endif 785 786 call printmeminfo("itrdrv after checkpoint"//char(0)) 787 else if(stopjob.eq.-2) then 788 if(myrank.eq.master) then 789 write(*,*) 'line 755 says no write before stopping' 790 write(*,*) 'istep,nstep,irs',istep,nstep(itseq),irs 791 endif 792 endif !just the instantaneous stuff for videos 793c 794c.... compute the consistent boundary flux 795c 796 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 797 call Bflux ( yold, acold, uold, x, 798 & shp, shgl, shpb, 799 & shglb, ilwork, iBC, 800 & BC, iper, wallssVec) 801 endif 802 803 if(stopjob.eq.-2) goto 2003 804 805 806c 807c ... update the flow history for the impedance convolution, filter it and write it out 808c 809 if(numImpSrfs.gt.zero) then 810 call UpdHistConv(y,nsrflistImp,numImpSrfs) !uses Delt(1) 811 endif 812 813c 814c ... update the flow history for the RCR convolution 815c 816 if(numRCRSrfs.gt.zero) then 817 call UpdHistConv(y,nsrflistRCR,numRCRSrfs) !uses lstep 818 endif 819 820 821c... dump TIME SERIES 822 823 if (exts) then 824 if (mod(lstep-1,freq).eq.0) then 825 826 if (numpe > 1) then 827 do jj = 1, ntspts 828 vartssoln((jj-1)*ndof+1:jj*ndof)=varts(jj,:) 829 ivarts=zero 830 enddo 831 do k=1,ndof*ntspts 832 if(vartssoln(k).ne.zero) ivarts(k)=1 833 enddo 834 835! call MPI_REDUCE(vartssoln, vartssolng, ndof*ntspts, 836! & MPI_DOUBLE_PRECISION, MPI_SUM, master, 837! & MPI_COMM_WORLD, ierr) 838 839 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 840 call MPI_ALLREDUCE(vartssoln, vartssolng, 841 & ndof*ntspts, 842 & MPI_DOUBLE_PRECISION, MPI_SUM, 843 & MPI_COMM_WORLD, ierr) 844 845! call MPI_REDUCE(ivarts, ivartsg, ndof*ntspts, 846! & MPI_INTEGER, MPI_SUM, master, 847! & MPI_COMM_WORLD, ierr) 848 849 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 850 call MPI_ALLREDUCE(ivarts, ivartsg, ndof*ntspts, 851 & MPI_INTEGER, MPI_SUM, 852 & MPI_COMM_WORLD, ierr) 853 854! if (myrank.eq.zero) then 855 do jj = 1, ntspts 856 857 if(myrank .eq. iv_rank(jj)) then 858 ! No need to update all varts components, only the one treated by the expected rank 859 ! Note: keep varts as a vector, as multiple probes could be treated by one rank 860 indxvarts = (jj-1)*ndof 861 do k=1,ndof 862 if(ivartsg(indxvarts+k).ne.0) then ! none of the vartssoln(parts) were non zero 863 varts(jj,k)=vartssolng(indxvarts+k)/ 864 & ivartsg(indxvarts+k) 865 endif 866 enddo 867 endif !only if myrank eq iv_rank(jj) 868 enddo 869! endif !only on master 870 endif !only if numpe > 1 871 872! if (myrank.eq.zero) then 873 do jj = 1, ntspts 874 if(myrank .eq. iv_rank(jj)) then 875 ifile = 1000+jj 876 write(ifile,555) lstep, (varts(jj,k),k=1,ndof) !Beware of format 555 - check ndof 877c call flush(ifile) 878 if (((irs .ge. 1) .and. 879 & (mod(lstep, ntout) .eq. 0))) then 880 close(ifile) 881 fvarts='varts/varts' 882 fvarts=trim(fvarts)//trim(cname2(jj)) 883 fvarts=trim(fvarts)//trim(cname2(lskeep)) 884 fvarts=trim(fvarts)//'.dat' 885 fvarts=trim(fvarts) 886 open(unit=ifile, file=fvarts, 887 & position='append') 888 endif !only when dumping restart 889 endif 890 enddo 891! endif !only on master 892 893 varts(:,:) = zero ! reset the array for next step 894 895 896!555 format(i6,5(2x,E12.5e2)) 897555 format(i6,6(2x,E20.12e2)) !assuming ndof = 6 here 898 899 endif 900 endif 901 902c 903c.... update and the aerodynamic forces 904c 905 call forces ( yold, ilwork ) 906 907 if((irscale.ge.0).or.(itwmod.gt.0)) 908 & call getvel (yold, ilwork, iBC, 909 & nsons, ifath, velbar) 910 911 if((irscale.ge.0).and.(myrank.eq.master)) then 912 call genscale(yold, x, iper, 913 & iBC, ifath, velbar, 914 & nsons) 915 endif 916c 917c.... print out results. 918c 919 ntoutv=max(ntout,100) ! velb is not needed so often 920 if ((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) then 921 if( (mod(lstep, ntoutv) .eq. 0) .and. 922 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 923 & ((nsonmax.eq.1).and.(iLES.gt.0)))) 924 & call rwvelb ('out ', velbar ,ifail) 925 endif 926c 927c.... end of the NSTEP and NTSEQ loops 928c 929 930 931c 932c.... -------------------> error calculation <----------------- 933c 934 if(ierrcalc.eq.1 .or. ioybar.eq.1) then 935c$$$c 936c$$$c compute average 937c$$$c 938c$$$ tfact=one/istep 939c$$$ ybar =tfact*yold + (one-tfact)*ybar 940 941c compute average 942c ybar(:,1:3) are average velocity components 943c ybar(:,4) is average pressure 944c ybar(:,5) is average speed 945c ybar(:,6:8) is average of sq. of each vel. component 946c ybar(:,9) is average of sq. of pressure 947c ybar(:,10:12) is average of cross vel. components : uv, uw and vw 948c averaging procedure justified only for identical time step sizes 949c ybar(:,13) is average of eddy viscosity 950c ybar(:,14:16) is average vorticity components 951c ybar(:,17) is average vorticity magnitude 952c istep is number of time step 953c 954 icollectybar = 0 955 if(nphasesincycle.eq.0 .or. 956 & istep.gt.ncycles_startphaseavg*nstepsincycle) then 957 icollectybar = 1 958 if((istep-1).eq.ncycles_startphaseavg*nstepsincycle) 959 & istepsinybar = 0 ! init. to zero in first cycle in avg. 960 endif 961 962 if(icollectybar.eq.1) then 963 istepsinybar = istepsinybar+1 964 tfact=one/istepsinybar 965 966 if(myrank.eq.master .and. nphasesincycle.ne.0 .and. 967 & mod((istep-1),nstepsincycle).eq.0) 968 & write(*,*)'nsamples in phase average:',istepsinybar 969 970c ybar to contain the averaged ((u,v,w),p)-fields 971c and speed average, i.e., sqrt(u^2+v^2+w^2) 972c and avg. of sq. terms including 973c u^2, v^2, w^2, p^2 and cross terms of uv, uw and vw 974 975 ybar(:,1) = tfact*yold(:,1) + (one-tfact)*ybar(:,1) 976 ybar(:,2) = tfact*yold(:,2) + (one-tfact)*ybar(:,2) 977 ybar(:,3) = tfact*yold(:,3) + (one-tfact)*ybar(:,3) 978 ybar(:,4) = tfact*yold(:,4) + (one-tfact)*ybar(:,4) 979 ybar(:,5) = tfact*sqrt(yold(:,1)**2+yold(:,2)**2+ 980 & yold(:,3)**2) + (one-tfact)*ybar(:,5) 981 ybar(:,6) = tfact*yold(:,1)**2 + 982 & (one-tfact)*ybar(:,6) 983 ybar(:,7) = tfact*yold(:,2)**2 + 984 & (one-tfact)*ybar(:,7) 985 ybar(:,8) = tfact*yold(:,3)**2 + 986 & (one-tfact)*ybar(:,8) 987 ybar(:,9) = tfact*yold(:,4)**2 + 988 & (one-tfact)*ybar(:,9) 989 ybar(:,10) = tfact*yold(:,1)*yold(:,2) + !uv 990 & (one-tfact)*ybar(:,10) 991 ybar(:,11) = tfact*yold(:,1)*yold(:,3) + !uw 992 & (one-tfact)*ybar(:,11) 993 ybar(:,12) = tfact*yold(:,2)*yold(:,3) + !vw 994 & (one-tfact)*ybar(:,12) 995!MR CHANGE 996 if(nsclr.gt.0) !nut 997 & ybar(:,13) = tfact*yold(:,6) + (one-tfact)*ybar(:,13) 998 999 if(ivort == 1) then !vorticity 1000 ybar(:,14) = tfact*vorticity(:,1) + 1001 & (one-tfact)*ybar(:,14) 1002 ybar(:,15) = tfact*vorticity(:,2) + 1003 & (one-tfact)*ybar(:,15) 1004 ybar(:,16) = tfact*vorticity(:,3) + 1005 & (one-tfact)*ybar(:,16) 1006 ybar(:,17) = tfact*vorticity(:,4) + 1007 & (one-tfact)*ybar(:,17) 1008 endif 1009 1010 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1011 wallssVecBar(:,1) = tfact*wallssVec(:,1) 1012 & +(one-tfact)*wallssVecBar(:,1) 1013 wallssVecBar(:,2) = tfact*wallssVec(:,2) 1014 & +(one-tfact)*wallssVecBar(:,2) 1015 wallssVecBar(:,3) = tfact*wallssVec(:,3) 1016 & +(one-tfact)*wallssVecBar(:,3) 1017 endif 1018!MR CHANGE END 1019 endif 1020c 1021c compute phase average 1022c 1023 if(nphasesincycle.ne.0 .and. 1024 & istep.gt.ncycles_startphaseavg*nstepsincycle) then 1025 1026c beginning of cycle is considered as ncycles_startphaseavg*nstepsincycle+1 1027 if((istep-1).eq.ncycles_startphaseavg*nstepsincycle) 1028 & icyclesinavg = 0 ! init. to zero in first cycle in avg. 1029 1030 ! find number of steps between phases 1031 nstepsbtwphase = nstepsincycle/nphasesincycle ! integer value 1032 if(mod(istep-1,nstepsincycle).eq.0) then 1033 iphase = 1 ! init. to one in beginning of every cycle 1034 icyclesinavg = icyclesinavg + 1 1035 endif 1036 1037 icollectphase = 0 1038 istepincycle = mod(istep,nstepsincycle) 1039 if(istepincycle.eq.0) istepincycle=nstepsincycle 1040 if(istepincycle.eq.iphase*nstepsbtwphase) then 1041 icollectphase = 1 1042 iphase = iphase+1 ! use 'iphase-1' below 1043 endif 1044 1045 if(icollectphase.eq.1) then 1046 tfactphase = one/icyclesinavg 1047 1048 if(myrank.eq.master) then 1049 write(*,*) 'nsamples in phase ',iphase-1,': ', 1050 & icyclesinavg 1051 endif 1052 1053 yphbar(:,1,iphase-1) = tfactphase*yold(:,1) + 1054 & (one-tfactphase)*yphbar(:,1,iphase-1) 1055 yphbar(:,2,iphase-1) = tfactphase*yold(:,2) + 1056 & (one-tfactphase)*yphbar(:,2,iphase-1) 1057 yphbar(:,3,iphase-1) = tfactphase*yold(:,3) + 1058 & (one-tfactphase)*yphbar(:,3,iphase-1) 1059 yphbar(:,4,iphase-1) = tfactphase*yold(:,4) + 1060 & (one-tfactphase)*yphbar(:,4,iphase-1) 1061 yphbar(:,5,iphase-1) = tfactphase*sqrt(yold(:,1)**2 1062 & +yold(:,2)**2+yold(:,3)**2) + 1063 & (one-tfactphase)*yphbar(:,5,iphase-1) 1064!MR CHANGE 1065 yphbar(:,6,iphase-1) = 1066 & tfactphase*yold(:,1)*yold(:,1) 1067 & +(one-tfactphase)*yphbar(:,6,iphase-1) 1068 1069 yphbar(:,7,iphase-1) = 1070 & tfactphase*yold(:,1)*yold(:,2) 1071 & +(one-tfactphase)*yphbar(:,7,iphase-1) 1072 1073 yphbar(:,8,iphase-1) = 1074 & tfactphase*yold(:,1)*yold(:,3) 1075 & +(one-tfactphase)*yphbar(:,8,iphase-1) 1076 1077 yphbar(:,9,iphase-1) = 1078 & tfactphase*yold(:,2)*yold(:,2) 1079 & +(one-tfactphase)*yphbar(:,9,iphase-1) 1080 1081 yphbar(:,10,iphase-1) = 1082 & tfactphase*yold(:,2)*yold(:,3) 1083 & +(one-tfactphase)*yphbar(:,10,iphase-1) 1084 1085 yphbar(:,11,iphase-1) = 1086 & tfactphase*yold(:,3)*yold(:,3) 1087 & +(one-tfactphase)*yphbar(:,11,iphase-1) 1088 1089 if(ivort == 1) then 1090 yphbar(:,12,iphase-1) = 1091 & tfactphase*vorticity(:,1) 1092 & +(one-tfactphase)*yphbar(:,12,iphase-1) 1093 yphbar(:,13,iphase-1) = 1094 & tfactphase*vorticity(:,2) 1095 & +(one-tfactphase)*yphbar(:,13,iphase-1) 1096 yphbar(:,14,iphase-1) = 1097 & tfactphase*vorticity(:,3) 1098 & +(one-tfactphase)*yphbar(:,14,iphase-1) 1099 yphbar(:,15,iphase-1) = 1100 & tfactphase*vorticity(:,4) 1101 & +(one-tfactphase)*yphbar(:,15,iphase-1) 1102 endif 1103!MR CHANGE END 1104 endif 1105 endif 1106c 1107c compute rms 1108c 1109 if(icollectybar.eq.1) then 1110 rerr(:, 7)=rerr(:, 7)+(yold(:,1)-ybar(:,1))**2 1111 rerr(:, 8)=rerr(:, 8)+(yold(:,2)-ybar(:,2))**2 1112 rerr(:, 9)=rerr(:, 9)+(yold(:,3)-ybar(:,3))**2 1113 rerr(:,10)=rerr(:,10)+(yold(:,4)-ybar(:,4))**2 1114 endif 1115 endif 1116 2003 continue ! we get here if stopjob equals lstep and this jumped over 1117! the statistics computation because we have no new data to average in 1118! rather we are just trying to output the last state that was not already 1119! written 1120c 1121c.... ----------------------> Complete Restart Processing <---------------------- 1122c 1123! for now it is the same frequency but need to change this 1124! soon.... but don't forget to change the field counter in 1125! new_interface.cc 1126! 1127 if (((irs .ge. 1) .and. (mod(lstep, ntout) .eq. 0)) .or. 1128 & istep.eq.nstep(itseq)) then 1129 if ((irs .ge. 1) .and. ((mod(lstep, ntout) .ne. 0) .or. 1130 & (nstp .eq. 0))) then 1131 if( 1132 & ((irscale.ge.0).or.(itwmod.gt.0) .or. 1133 & ((nsonmax.eq.1).and.iLES.gt.0))) 1134 & call rwvelb ('out ', velbar ,ifail) 1135 endif 1136 1137 lesId = numeqns(1) 1138 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1139 if(myrank.eq.0) then 1140 tcormr1 = TMRC() 1141 endif 1142 if(nsolflow.eq.1) then 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 endif 1152 1153 if(ierrcalc.eq.1) then 1154c 1155c.....smooth the error indicators 1156c 1157 do i=1,ierrsmooth 1158 call errsmooth( rerr, x, iper, ilwork, shp, shgl, iBC ) 1159 end do 1160 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1161 if(myrank.eq.0) then 1162 tcormr1 = TMRC() 1163 endif 1164 call write_error(myrank, lstep, nshg, 10, rerr ) 1165 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1166 if(myrank.eq.0) then 1167 tcormr2 = TMRC() 1168 write(6,*) 'Time to write the error fields to the disks', 1169 & tcormr2-tcormr1 1170 endif 1171 endif ! ierrcalc 1172 1173 if(ioybar.eq.1) then 1174 if(ivort == 1) then 1175 call write_field(myrank,'a','ybar',4, 1176 & ybar,'d',nshg,17,lstep) 1177 else 1178 call write_field(myrank,'a','ybar',4, 1179 & ybar,'d',nshg,13,lstep) 1180 endif 1181 1182 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1183 call write_field(myrank,'a','wssbar',6, 1184 & wallssVecBar,'d',nshg,3,lstep) 1185 endif 1186 1187 if(nphasesincycle .gt. 0) then 1188 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1189 if(myrank.eq.0) then 1190 tcormr1 = TMRC() 1191 endif 1192 do iphase=1,nphasesincycle 1193 if(ivort == 1) then 1194 call write_phavg2(myrank,'a','phase_average',13,iphase, 1195 & nphasesincycle,yphbar(:,:,iphase),'d',nshg,15,lstep) 1196 else 1197 call write_phavg2(myrank,'a','phase_average',13,iphase, 1198 & nphasesincycle,yphbar(:,:,iphase),'d',nshg,11,lstep) 1199 endif 1200 end do 1201 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1202 if(myrank.eq.0) then 1203 tcormr2 = TMRC() 1204 write(6,*) 'write all phase avg to the disks = ', 1205 & tcormr2-tcormr1 1206 endif 1207 endif !nphasesincyle 1208 endif !ioybar 1209 1210 if ( ( ihessian .eq. 1 ) .and. ( numpe < 2 ) )then 1211 uhess = zero 1212 gradu = zero 1213 tf = zero 1214 1215 do ku=1,nshg 1216 tf(ku,1) = x(ku,1)**3 1217 end do 1218 call hessian( yold, x, shp, shgl, iBC, 1219 & shpb, shglb, iper, ilwork, uhess, gradu ) 1220 1221 call write_hessian( uhess, gradu, nshg ) 1222 endif 1223 1224 if(iRANS.lt.0) then 1225 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1226 if(myrank.eq.0) then 1227 tcormr1 = TMRC() 1228 endif 1229 call write_field(myrank,'a','dwal',4,d2wall,'d', 1230 & nshg,1,lstep) 1231 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1232 if(myrank.eq.0) then 1233 tcormr2 = TMRC() 1234 write(6,*) 'Time to write dwal to the disks = ', 1235 & tcormr2-tcormr1 1236 endif 1237 endif !iRANS 1238 1239 endif ! write out complete restart state 1240 !next 2 lines are two ways to end early 1241 if(stopjob.eq.-2) goto 2002 1242 if(istop.eq.1000) goto 2002 ! stop when delta small (see rstatic) 1243 2000 continue 1244 2002 continue 1245 1246! done with time stepping so deallocate fields already written 1247! 1248 if(ioybar.eq.1) then 1249 deallocate(ybar) 1250 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1251 deallocate(wallssVecbar) 1252 endif 1253 if(nphasesincycle .gt. 0) then 1254 deallocate(yphbar) 1255 endif !nphasesincyle 1256 endif !ioybar 1257 if(ivort == 1) then 1258 deallocate(strain,vorticity) 1259 endif 1260 if(abs(itwmod).ne.1 .and. iowflux.eq.1) then 1261 deallocate(wallssVec) 1262 endif 1263 if(iRANS.lt.0) then 1264 deallocate(d2wall) 1265 endif 1266 1267 1268CAD tcorecp2 = second(0) 1269CAD tcorewc2 = second(-1) 1270 1271CAD write(6,*) 'T(core) cpu-wallclock = ',tcorecp2-tcorecp1, 1272CAD & tcorewc2-tcorewc1 1273 1274 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1275 if(myrank.eq.0) then 1276 tcorecp2 = TMRC() 1277 write(6,*) 'T(core) cpu = ',tcorecp2-tcorecp1 1278 write(6,*) '(Elm. form.',tcorecp(1), 1279 & ',Lin. alg. sol.',tcorecp(2),')' 1280 write(6,*) '(Elm. form. Scal.',tcorecpscal(1), 1281 & ',Lin. alg. sol. Scal.',tcorecpscal(2),')' 1282 write(6,*) '' 1283 1284 endif 1285 1286 call print_system_stats(tcorecp, tcorecpscal) 1287 call print_mesh_stats() 1288 call print_mpi_stats() 1289 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1290! return 1291c call MPI_Finalize() 1292c call MPI_ABORT(MPI_COMM_WORLD, ierr) 1293 1294 3000 continue 1295 1296 1297c 1298c.... close history and aerodynamic forces files 1299c 1300 if (myrank .eq. master) then 1301! close (ihist) 1302 close (iforce) 1303 close(76) 1304 if(numImpSrfs.gt.0 .or. numRCRSrfs.gt.0) then 1305 close (8177) 1306 endif 1307 endif 1308c 1309c.... close varts file for probes 1310c 1311 if(exts) then 1312 do jj=1,ntspts 1313 if (myrank == iv_rank(jj)) then 1314 close(1000+jj) 1315 endif 1316 enddo 1317 deallocate (ivarts) 1318 deallocate (ivartsg) 1319 deallocate (iv_rank) 1320 deallocate (vartssoln) 1321 deallocate (vartssolng) 1322 endif 1323 1324!MR CHANGE 1325 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1326 if(myrank.eq.0) then 1327 write(*,*) 'itrdrv - done with aerodynamic forces' 1328 endif 1329!MR CHANGE 1330 1331 do isrf = 0,MAXSURF 1332! if ( nsrflist(isrf).ne.0 ) then 1333 if ( nsrflist(isrf).ne.0 .and. 1334 & myrank.eq.irankfilesforce(isrf)) then 1335 iunit=60+isrf 1336 close(iunit) 1337 endif 1338 enddo 1339 1340!MR CHANGE 1341 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1342 if(myrank.eq.0) then 1343 write(*,*) 'itrdrv - done with MAXSURF' 1344 endif 1345!MR CHANGE 1346 1347 1348 5 format(1X,F15.10,3X,F15.10,3X,F15.10,3X,F15.10) 1349 444 format(6(2x,e14.7)) 1350c 1351c.... end 1352c 1353 if(nsolflow.eq.1) then 1354 deallocate (lhsK) 1355 deallocate (lhsP) 1356 deallocate (aperm) 1357 deallocate (atemp) 1358 endif 1359 if(nsclrsol.gt.0) then 1360 deallocate (lhsS) 1361 deallocate (apermS) 1362 deallocate (atempS) 1363 endif 1364 1365 if(iabc==1) deallocate(acs) 1366 1367!MR CHANGE 1368 if (numpe > 1) call MPI_BARRIER(MPI_COMM_WORLD, ierr) 1369 if(myrank.eq.0) then 1370 write(*,*) 'itrdrv - done - BACK TO process.f' 1371 endif 1372!MR CHANGE 1373 1374 1375 1376 return 1377 end 1378 1379 subroutine lesmodels(y, ac, shgl, shp, 1380 & iper, ilwork, rowp, colm, 1381 & nsons, ifath, x, 1382 & iBC, BC) 1383 1384 include "common.h" 1385 1386 real*8 y(nshg,ndof), ac(nshg,ndof), 1387 & x(numnp,nsd), 1388 & BC(nshg,ndofBC) 1389 real*8 shp(MAXTOP,maxsh,MAXQPT), 1390 & shgl(MAXTOP,nsd,maxsh,MAXQPT) 1391 1392c 1393 integer rowp(nshg,nnz), colm(nshg+1), 1394 & iBC(nshg), 1395 & ilwork(nlwork), 1396 & iper(nshg) 1397 dimension ifath(numnp), nsons(nfath) 1398 1399 real*8, allocatable, dimension(:) :: fwr2,fwr3,fwr4 1400 real*8, allocatable, dimension(:) :: stabdis,cdelsq1 1401 real*8, allocatable, dimension(:,:) :: xavegt, xavegt2,xavegt3 1402 1403 if( (iLES.gt.1) ) then ! Allocate Stuff for advanced LES models 1404 allocate (fwr2(nshg)) 1405 allocate (fwr3(nshg)) 1406 allocate (fwr4(nshg)) 1407 allocate (xavegt(nfath,12)) 1408 allocate (xavegt2(nfath,12)) 1409 allocate (xavegt3(nfath,12)) 1410 allocate (stabdis(nfath)) 1411 endif 1412 1413c.... get dynamic model coefficient 1414c 1415 ilesmod=iLES/10 1416c 1417c digit bit set filter rule, 10 bit set model 1418c 1419 if (ilesmod.eq.0) then ! 0 < iLES< 10 => dyn. model calculated 1420 ! at nodes based on discrete filtering 1421 1422 1423 if(isubmod.eq.2) then 1424 call SUPGdis(y, ac, shgl, shp, 1425 & iper, ilwork, 1426 & nsons, ifath, x, 1427 & iBC, BC, stabdis, xavegt3) 1428 endif 1429 1430 if( ((isubmod.eq.0).or.(isubmod.eq.2)))then ! If no 1431 ! sub-model 1432 ! or SUPG 1433 ! model wanted 1434 1435 if(i2filt.eq.0)then ! If simple filter 1436 1437 if(modlstats .eq. 0) then ! If no model stats wanted 1438 call getdmc (y, shgl, shp, 1439 & iper, ilwork, nsons, 1440 & ifath, x) 1441 else ! else get model stats 1442 call stdfdmc (y, shgl, shp, 1443 & iper, ilwork, nsons, 1444 & ifath, x) 1445 endif ! end of stats if statement 1446 1447 else ! else if twice filtering 1448 1449 call widefdmc(y, shgl, shp, 1450 & iper, ilwork, nsons, 1451 & ifath, x) 1452 1453 1454 endif ! end of simple filter if statement 1455 1456 endif ! end of SUPG or no sub-model if statement 1457 1458 1459 if( (isubmod.eq.1) ) then ! If DFWR sub-model wanted 1460 call cdelBHsq (y, shgl, shp, 1461 & iper, ilwork, nsons, 1462 & ifath, x, cdelsq1) 1463 call FiltRat (y, shgl, shp, 1464 & iper, ilwork, nsons, 1465 & ifath, x, cdelsq1, 1466 & fwr4, fwr3) 1467 1468 1469 if (i2filt.eq.0) then ! If simple filter wanted 1470 call DFWRsfdmc(y, shgl, shp, 1471 & iper, ilwork, nsons, 1472 & ifath, x, fwr2, fwr3) 1473 else ! else if twice filtering wanted 1474 call DFWRwfdmc(y, shgl, shp, 1475 & iper, ilwork, nsons, 1476 & ifath, x, fwr4, fwr4) 1477 endif ! end of simple filter if statement 1478 1479 endif ! end of DFWR sub-model if statement 1480 1481 if( (isubmod.eq.2) )then ! If SUPG sub-model wanted 1482 call dmcSUPG (y, ac, shgl, 1483 & shp, iper, ilwork, 1484 & nsons, ifath, x, 1485 & iBC, BC, rowp, colm, 1486 & xavegt2, stabdis) 1487 endif 1488 1489 if(idis.eq.1)then ! If SUPG/Model dissipation wanted 1490 call ediss (y, ac, shgl, 1491 & shp, iper, ilwork, 1492 & nsons, ifath, x, 1493 & iBC, BC, xavegt) 1494 endif 1495 1496 endif ! end of ilesmod 1497 1498 if (ilesmod .eq. 1) then ! 10 < iLES < 20 => dynamic-mixed 1499 ! at nodes based on discrete filtering 1500 call bardmc (y, shgl, shp, 1501 & iper, ilwork, 1502 & nsons, ifath, x) 1503 endif 1504 1505 if (ilesmod .eq. 2) then ! 20 < iLES < 30 => dynamic at quad 1506 ! pts based on lumped projection filt. 1507 1508 if(isubmod.eq.0)then 1509 call projdmc (y, shgl, shp, 1510 & iper, ilwork, x) 1511 else 1512 call cpjdmcnoi (y, shgl, shp, 1513 & iper, ilwork, x, 1514 & rowp, colm, 1515 & iBC, BC) 1516 endif 1517 1518 endif 1519 1520 if( (iLES.gt.1) ) then ! Deallocate Stuff for advanced LES models 1521 deallocate (fwr2) 1522 deallocate (fwr3) 1523 deallocate (fwr4) 1524 deallocate (xavegt) 1525 deallocate (xavegt2) 1526 deallocate (xavegt3) 1527 deallocate (stabdis) 1528 endif 1529 return 1530 end 1531 1532c 1533c...initialize the coefficients for the impedance convolution 1534c 1535 subroutine CalcImpConvCoef (numISrfs, numTpoints) 1536 1537 use convolImpFlow !uses flow history and impedance for convolution 1538 1539 include "common.h" !for alfi 1540 1541 integer numISrfs, numTpoints 1542 1543 allocate (ConvCoef(numTpoints+2,3)) !same time discret. for all imp. BC 1544 do j=1,numTpoints+2 1545 ConvCoef(j,:)=0.5/numTpoints !dt/2 divided by period T=N*dt 1546 ConvCoef(j,1)=ConvCoef(j,1)*(1.0-alfi)*(1.0-alfi) 1547 ConvCoef(j,2)=ConvCoef(j,2)*(1.0+2*alfi*(1.0-alfi)) 1548 ConvCoef(j,3)=ConvCoef(j,3)*alfi*alfi 1549 enddo 1550 ConvCoef(1,2)=zero 1551 ConvCoef(1,3)=zero 1552 ConvCoef(2,3)=zero 1553 ConvCoef(numTpoints+1,1)=zero 1554 ConvCoef(numTpoints+2,2)=zero 1555 ConvCoef(numTpoints+2,1)=zero 1556c 1557c...calculate the coefficients for the impedance convolution 1558c 1559 allocate (ImpConvCoef(numTpoints+2,numISrfs)) 1560 1561c..coefficients below assume Q linear in time step, Z constant 1562c do j=3,numTpoints 1563c ImpConvCoef(j,:) = ValueListImp(j-1,:)*ConvCoef(j,3) 1564c & + ValueListImp(j,:)*ConvCoef(j,2) 1565c & + ValueListImp(j+1,:)*ConvCoef(j,1) 1566c enddo 1567c ImpConvCoef(1,:) = ValueListImp(2,:)*ConvCoef(1,1) 1568c ImpConvCoef(2,:) = ValueListImp(2,:)*ConvCoef(2,2) 1569c & + ValueListImp(3,:)*ConvCoef(2,1) 1570c ImpConvCoef(numTpoints+1,:) = 1571c & ValueListImp(numTpoints,:)*ConvCoef(numTpoints+1,3) 1572c & + ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+1,2) 1573c ImpConvCoef(numTpoints+2,:) = 1574c & ValueListImp(numTpoints+1,:)*ConvCoef(numTpoints+2,3) 1575 1576c..try easiest convolution Q and Z constant per time step 1577 do j=3,numTpoints+1 1578 ImpConvCoef(j,:) = ValueListImp(j-1,:)/numTpoints 1579 enddo 1580 ImpConvCoef(1,:) =zero 1581 ImpConvCoef(2,:) =zero 1582 ImpConvCoef(numTpoints+2,:) = 1583 & ValueListImp(numTpoints+1,:)/numTpoints 1584c compensate for yalpha passed not y in Elmgmr() 1585 ImpConvCoef(numTpoints+1,:)= ImpConvCoef(numTpoints+1,:) 1586 & - ImpConvCoef(numTpoints+2,:)*(1.0-alfi)/alfi 1587 ImpConvCoef(numTpoints+2,:)= ImpConvCoef(numTpoints+2,:)/alfi 1588 return 1589 end 1590 1591c 1592c ... update the flow rate history for the impedance convolution, filter it and write it out 1593c 1594 subroutine UpdHistConv(y,nsrfIdList,numSrfs) 1595 1596 use convolImpFlow !brings ntimeptpT, QHistImp, QHistTry, QHistTryF, numImpSrfs 1597 use convolRCRFlow !brings QHistRCR, numRCRSrfs 1598 1599 include "common.h" 1600 1601 integer nsrfIdList(0:MAXSURF), numSrfs 1602 character*20 fname1 1603 character*10 cname2 1604 character*5 cname 1605 real*8 y(nshg,3) !velocity at time n+1 1606 real*8 NewQ(0:MAXSURF) !temporary unknown for the flow rate 1607 !that needs to be added to the flow history 1608 1609 call GetFlowQ(NewQ,y,nsrfIdList,numSrfs) !new flow at time n+1 1610c 1611c... for imp BC: shift QHist, add new constribution, filter and write out 1612c 1613 if(numImpSrfs.gt.zero) then 1614 do j=1, ntimeptpT 1615 QHistImp(j,1:numSrfs)=QHistImp(j+1,1:numSrfs) 1616 enddo 1617 QHistImp(ntimeptpT+1,1:numSrfs) = NewQ(1:numSrfs) 1618 1619c 1620c....filter the flow rate history 1621c 1622 cutfreq = 10 !hardcoded cutting frequency of the filter 1623 do j=1, ntimeptpT 1624 QHistTry(j,:)=QHistImp(j+1,:) 1625 enddo 1626 call Filter(QHistTryF,QHistTry,ntimeptpT,Delt(1),cutfreq) 1627c.... if no filter applied then uncomment next three lines 1628c do j=1, ntimeptpT 1629c QHistTryF(j,:)=QHistTry(j,:) 1630c enddo 1631 1632c QHistImp(1,:)=zero ! why do we do this? for beta(1,:) = zero it does not really matters 1633 do j=1, ntimeptpT 1634 QHistImp(j+1,:)=QHistTryF(j,:) 1635 enddo 1636c 1637c.... write out the new history of flow rates to Qhistor.dat 1638c 1639 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1640 & (istep .eq. nstep(1)))) .and. 1641 & (myrank .eq. master)) then 1642 open(unit=816, file='Qhistor.dat',status='replace') 1643 write(816,*) ntimeptpT 1644 do j=1,ntimeptpT+1 1645 write(816,*) (QHistImp(j,n),n=1, numSrfs) 1646 enddo 1647 close(816) 1648c... write out a copy with step number to be able to restart 1649 fname1 = 'Qhistor' 1650 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1651 open(unit=8166,file=trim(fname1),status='unknown') 1652 write(8166,*) ntimeptpT 1653 do j=1,ntimeptpT+1 1654 write(8166,*) (QHistImp(j,n),n=1, numSrfs) 1655 enddo 1656 close(8166) 1657 endif 1658 endif 1659 1660c 1661c... for RCR bc just add the new contribution 1662c 1663 if(numRCRSrfs.gt.zero) then 1664 QHistRCR(lstep+1,1:numSrfs) = NewQ(1:numSrfs) 1665c 1666c.... write out the new history of flow rates to Qhistor.dat 1667c 1668 if ((irs .ge. 1) .and. (myrank .eq. master)) then 1669 if(istep.eq.1) then 1670 open(unit=816,file='Qhistor.dat',status='unknown') 1671 else 1672 open(unit=816,file='Qhistor.dat',position='append') 1673 endif 1674 if(istep.eq.1) then 1675 do j=1,lstep 1676 write(816,*) j, (QHistRCR(j,n),n=1,numSrfs) ! read from file of previous run 1677 enddo 1678 endif 1679 write(816,*) lstep+1, (QHistRCR(lstep+1,n),n=1, numSrfs) 1680 close(816) 1681c... write out a copy with step number to be able to restart 1682 if (((irs .ge. 1) .and. ((mod(lstep, ntout) .eq. 0) .or. 1683 & (istep .eq. nstep(1)))) .and. 1684 & (myrank .eq. master)) then 1685 fname1 = 'Qhistor' 1686 fname1 = trim(fname1)//trim(cname2(lstep))//'.dat' 1687 open(unit=8166,file=trim(fname1),status='unknown') 1688 write(8166,*) lstep+1 1689 do j=1,lstep+1 1690 write(8166,*) (QHistRCR(j,n),n=1, numSrfs) 1691 enddo 1692 close(8166) 1693 endif 1694 endif 1695 endif 1696 1697 return 1698 end 1699 1700c 1701c...calculate the time varying coefficients for the RCR convolution 1702c 1703 subroutine CalcRCRConvCoef (stepn, numSrfs) 1704 1705 use convolRCRFlow !brings in ValueListRCR, dtRCR 1706 1707 include "common.h" !brings alfi 1708 1709 integer numSrfs, stepn 1710 1711 RCRConvCoef = zero 1712 if (stepn .eq. 0) then 1713 RCRConvCoef(1,:) = ValueListRCR(1,:)*(1.0-alfi) + 1714 & ValueListRCR(3,:)*(-alfi + 1.0 + 1/dtRCR(:) 1715 & - exp(-alfi*dtRCR(:))*(1 + 1/dtRCR(:))) 1716 RCRConvCoef(2,:) = ValueListRCR(1,:)*alfi 1717 & + ValueListRCR(3,:) 1718 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1719 endif 1720 if (stepn .ge. 1) then 1721 RCRConvCoef(1,:) =-ValueListRCR(3,:)*exp(-dtRCR(:)*(stepn+alfi)) 1722 & *(1 + (1 - exp(dtRCR(:)))/dtRCR(:)) 1723 RCRConvCoef(stepn+1,:) = ValueListRCR(1,:)*(1-alfi) 1724 & - ValueListRCR(3,:)*(alfi - 1 - 1/dtRCR(:) 1725 & + exp(-alfi*dtRCR(:))/dtRCR(:)*(2 - exp(-dtRCR(:)))) 1726 RCRConvCoef(stepn+2,:) = ValueListRCR(1,:)*alfi 1727 & + ValueListRCR(3,:) 1728 & *(alfi - 1/dtRCR(:) + exp(-alfi*dtRCR(:))/dtRCR(:)) 1729 endif 1730 if (stepn .ge. 2) then 1731 do j=2,stepn 1732 RCRConvCoef(j,:) = ValueListRCR(3,:)/dtRCR(:)* 1733 & exp(-dtRCR(:)*(stepn + alfi + 2 - j))* 1734 & (1 - exp(dtRCR(:)))**2 1735 enddo 1736 endif 1737 1738c compensate for yalpha passed not y in Elmgmr() 1739 RCRConvCoef(stepn+1,:)= RCRConvCoef(stepn+1,:) 1740 & - RCRConvCoef(stepn+2,:)*(1.0-alfi)/alfi 1741 RCRConvCoef(stepn+2,:)= RCRConvCoef(stepn+2,:)/alfi 1742 1743 return 1744 end 1745 1746c 1747c...calculate the time dependent H operator for the RCR convolution 1748c 1749 subroutine CalcHopRCR (timestepRCR, stepn, numSrfs) 1750 1751 use convolRCRFlow !brings in HopRCR, dtRCR 1752 1753 include "common.h" 1754 1755 integer numSrfs, stepn 1756 real*8 PdistCur(0:MAXSURF), timestepRCR 1757 1758 HopRCR=zero 1759 call RCRint(timestepRCR*(stepn + alfi),PdistCur) 1760 HopRCR(1:numSrfs) = RCRic(1:numSrfs) 1761 & *exp(-dtRCR(1:numSrfs)*(stepn + alfi)) + PdistCur(1:numSrfs) 1762 return 1763 end 1764c 1765c ... initialize the influence of the initial conditions for the RCR BC 1766c 1767 subroutine calcRCRic(y,srfIdList,numSrfs) 1768 1769 use convolRCRFlow !brings RCRic, ValueListRCR, ValuePdist 1770 1771 include "common.h" 1772 1773 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled 1774 real*8 y(nshg,4) !need velocity and pressure 1775 real*8 Qini(0:MAXSURF) !initial flow rate 1776 real*8 PdistIni(0:MAXSURF) !initial distal pressure 1777 real*8 Pini(0:MAXSURF),CoupleArea(0:MAXSURF) ! initial pressure 1778 real*8 VelOnly(nshg,3), POnly(nshg) 1779 1780 allocate (RCRic(0:MAXSURF)) 1781 1782 if(lstep.eq.0) then 1783 VelOnly(:,1:3)=y(:,1:3) 1784 call GetFlowQ(Qini,VelOnly,srfIdList,numSrfs) !get initial flow 1785 QHistRCR(1,1:numSrfs)=Qini(1:numSrfs) !initialize QHistRCR 1786 POnly(:)=y(:,4) ! pressure 1787 call integrScalar(Pini,POnly,srfIdList,numSrfs) !get initial pressure integral 1788 POnly(:)=one ! one to get area 1789 call integrScalar(CoupleArea,POnly,srfIdList,numSrfs) !get surf area 1790 Pini(1:numSrfs) = Pini(1:numSrfs)/CoupleArea(1:numSrfs) 1791 else 1792 Qini(1:numSrfs)=QHistRCR(1,1:numSrfs) 1793 Pini(1:numSrfs)=zero ! hack 1794 endif 1795 call RCRint(istep,PdistIni) !get initial distal P (use istep) 1796 RCRic(1:numSrfs) = Pini(1:numSrfs) 1797 & - ValueListRCR(1,:)*Qini(1:numSrfs)-PdistIni(1:numSrfs) 1798 return 1799 end 1800 1801c.........function that integrates a scalar over a boundary 1802 subroutine integrScalar(scalInt,scal,srfIdList,numSrfs) 1803 1804 use pvsQbi !brings ndsurf, NASC 1805 1806 include "common.h" 1807 include "mpif.h" 1808 1809 integer srfIdList(0:MAXSURF), numSrfs, irankCoupled, i, k 1810 real*8 scal(nshg), scalInt(0:MAXSURF), scalIntProc(0:MAXSURF) 1811 1812 scalIntProc = zero 1813 do i = 1,nshg 1814 if(numSrfs.gt.zero) then 1815 do k = 1,numSrfs 1816 irankCoupled = 0 1817 if (srfIdList(k).eq.ndsurf(i)) then 1818 irankCoupled=k 1819 scalIntProc(irankCoupled) = scalIntProc(irankCoupled) 1820 & + NASC(i)*scal(i) 1821 exit 1822 endif 1823 enddo 1824 endif 1825 enddo 1826c 1827c at this point, each scalint has its "nodes" contributions to the scalar 1828c accumulated into scalIntProc. Note, because NASC is on processor this 1829c will NOT be the scalar for the surface yet 1830c 1831c.... reduce integrated scalar for each surface, push on scalInt 1832c 1833 npars=MAXSURF+1 1834 call MPI_ALLREDUCE (scalIntProc, scalInt(:), npars, 1835 & MPI_DOUBLE_PRECISION,MPI_SUM, MPI_COMM_WORLD,ierr) 1836 1837 return 1838 end 1839 1840 subroutine writeTimingMessage(key,iomode,timing) 1841 use iso_c_binding 1842 use phstr 1843 implicit none 1844 1845 character(len=*) :: key 1846 integer :: iomode 1847 real*8 :: timing 1848 character(len=1024) :: timing_msg 1849 character(len=*), parameter :: 1850 & streamModeString = c_char_"stream"//c_null_char, 1851 & fileModeString = c_char_"disk"//c_null_char 1852 1853 timing_msg = c_char_"Time to write "//c_null_char 1854 call phstr_appendStr(timing_msg,key) 1855 if ( iomode .eq. -1 ) then 1856 call phstr_appendStr(timing_msg, streamModeString) 1857 else 1858 call phstr_appendStr(timing_msg, fileModeString) 1859 endif 1860 call phstr_appendStr(timing_msg, c_char_' = '//c_null_char) 1861 call phstr_appendDbl(timing_msg, timing) 1862 write(6,*) trim(timing_msg) 1863 return 1864 end subroutine 1865 1866