1c readnblk.f (pronounce "Reed and Block Dot Eff") contains: 2c 3c module readarrays ("Red Arrays") -- contains the arrays that 4c are read in from binary files but not immediately blocked 5c through pointers. 6c 7c subroutine readnblk ("Reed and Block") -- allocates space for 8c and reads data to be contained in module readarrays. Reads 9c all remaining data and blocks them with pointers. 10c 11 12 13 module readarrays 14 15 real*8, allocatable :: point2x(:,:) 16 real*8, allocatable :: qold(:,:) 17 real*8, allocatable :: uold(:,:) 18 real*8, allocatable :: acold(:,:) 19 integer, allocatable :: iBCtmp(:) 20 real*8, allocatable :: BCinp(:,:) 21 22 integer, allocatable :: point2ilwork(:) 23 integer, allocatable :: nBC(:) 24 integer, allocatable :: point2iper(:) 25 integer, allocatable :: point2ifath(:) 26 integer, allocatable :: point2nsons(:) 27 28 end module 29 30 31 32 subroutine readnblk 33c 34 use readarrays 35 include "common.h" 36c 37 real*8, allocatable :: xread(:,:), qread(:,:), acread(:,:) 38 real*8, allocatable :: uread(:,:) 39 real*8, allocatable :: BCinpread(:,:) 40 integer, allocatable :: iperread(:), iBCtmpread(:) 41 integer, allocatable :: ilworkread(:), nBCread(:) 42 character*10 cname2 43 character*8 mach2 44!MR CHANGE 45! character*20 fmt1 46 character*30 fmt1 47!MR CHANGE END 48 character*255 fname1,fnamer,fnamelr 49 character*255 warning 50 integer igeomBAK, ibndc, irstin, ierr 51 integer intfromfile(50) ! integers read from headers 52 53cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 54ccccccccccccccccccccccccccccccccccccccccc New PhastaIO Definition Part ccccccccccccccccccccccccccccccccccccccccc 55 56 integer :: descriptor, descriptorG, GPID, color, nfiles, nfields 57 integer :: numparts, nppf 58 integer :: ierr_io, numprocs, itmp, itmp2 59 integer :: ignored 60 character*255 fname2, temp2 61 character*64 temp1 62 63cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 64cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 65 66 67c 68c.... determine the step number to start with 69c 70 open(unit=72,file='numstart.dat',status='old') 71 read(72,*) irstart 72 close(72) 73 lstep=irstart ! in case restart files have no fields 74 75c 76 fname1='geombc.dat' 77 fname1= trim(fname1) // cname2(myrank+1) 78c fnamelr='restart.latest' 79 80 itmp=1 81 if (irstart .gt. 0) itmp = int(log10(float(irstart)))+1 82 write (fmt1,"('(''restart.'',i',i1,',1x)')") itmp 83 write (fnamer,fmt1) irstart 84 fnamer = trim(fnamer) // cname2(myrank+1) 85c fnamelr = trim(fnamelr) // cname2(myrank+1) 86 87c 88c.... open input files 89c 90 91 92c call openfile( fname1, 'read?', igeomBAK ); 93 94 95c 96c.... try opening restart.latest.proc before trying restart.stepno.proc 97c 98c call openfile( fnamelr, 'read?', irstin ); 99c if ( irstin .eq. 0 ) 100 101!MR CHANGE 102! call openfile( fnamer, 'read?', irstin ); 103!MR CHANGE END 104 105! either one will work 106c 107c.... input the geometry parameters 108c 109 110cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 111!MR CHANGE 112 113! nfiles = 2 114! nfields = 31 115! numparts = 8 116! nppp = numparts/numpe 117! nppf = numparts/nfiles 118 119 nfiles = nsynciofiles 120! nfields = nsynciofieldsreadgeombc 121 numparts = numpe !This is the common settings. Beware if you try to compute several parts per process 122 123! nppp = numparts/numpe 124! nppf = numparts/nfiles 125!MR CHANGE END 126 127 itwo=2 128 ione=1 129 ieleven=11 130 itmp = int(log10(float(myrank+1)))+1 131 132 call phio_openfile_read('geombc-dat.' // char(0), nfiles, igeom); 133 134 call phio_readheader(igeom,'number of nodes' // char(0),numnp,ione, 135 & 'integer' // char(0), iotype) 136 137 call phio_readheader(igeom,'number of modes' // char(0),nshg,ione, 138 & 'integer' // char(0), iotype) 139 140 call phio_readheader(igeom,'number of interior elements' // char(0), 141 & numel,ione,'integer' // char(0), iotype) 142 143 call phio_readheader(igeom,'number of boundary elements' // char(0), 144 & numelb,ione,'integer' // char(0),iotype) 145 146 call phio_readheader(igeom,'maximum number of element nodes' // char(0), 147 & nen,ione,'integer' // char(0),iotype) 148 149 call phio_readheader(igeom,'number of interior tpblocks' // char(0), 150 & nelblk,ione,'integer' // char(0) ,iotype) 151 152 call phio_readheader(igeom,'number of boundary tpblocks' // char(0), 153 & nelblb,ione,'integer' // char(0), iotype) 154 155 call phio_readheader(igeom, 156 & 'number of nodes with Dirichlet BCs' // char(0), 157 & numpbc,ione,'integer' // char(0),iotype) 158 159 call phio_readheader(igeom,'number of shape functions' // char(0), 160 & ntopsh,ione,'integer' // char(0),iotype) 161 162c call closefile( igeom, "read" ) 163c call finalizephmpiio( igeom ) 164 165! if(myrank==0) then 166! print *, "Reading Finished, ********* : ", trim(fname2) 167! endif 168 169 170ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 171 172c ieleven=11 173c ione=1 174c fname1='number of nodes?' 175c call readheader(igeomBAK,fname1,numnp,ione,'integer', iotype) 176c fname1='number of modes?' 177c call readheader(igeomBAK,fname1,nshg,ione,'integer', iotype) 178cc fname1='number of global modes?' 179cc call readheader(igeomBAK,fname1,nshgt,ione,'integer', iotype) 180c fname1='number of interior elements?' 181c call readheader(igeomBAK,fname1,numel,ione,'integer', iotype) 182c fname1='number of boundary elements?' 183c call readheader(igeomBAK,fname1,numelb,ione,'integer', iotype) 184c fname1='maximum number of element nodes?' 185c call readheader(igeomBAK,fname1,nen,ione,'integer', iotype) 186c fname1='number of interior tpblocks?' 187c call readheader(igeomBAK,fname1,nelblk,ione,'integer', iotype) 188c fname1='number of boundary tpblocks?' 189c call readheader(igeomBAK,fname1,nelblb,ione,'integer', iotype) 190c fname1='number of nodes with Dirichlet BCs?' 191c call readheader(igeomBAK,fname1,numpbc,ione,'integer', iotype) 192c fname1='number of shape functions?' 193c call readheader(igeomBAK,fname1,ntopsh,ione,'integer', iotype) 194 195c 196c.... calculate the maximum number of boundary element nodes 197c 198 nenb = 0 199 do i = 1, melCat 200 if (nen .eq. nenCat(i,nsd)) nenb = max(nenCat(i,nsd-1), nenb) 201 enddo 202c 203 if (myrank == master) then 204 if (nenb .eq. 0) call error ('input ','nen ',nen) 205 endif 206c 207c.... setup some useful constants 208c 209 I3nsd = nsd / 3 ! nsd=3 integer flag 210 E3nsd = float(I3nsd) ! nsd=3 real flag 211c 212 if(matflg(1,1).lt.0) then 213 nflow = nsd + 1 214 else 215 nflow = nsd + 2 216 endif 217 ndof = nsd + 2 218 nsclr=impl(1)/100 219 ndof=ndof+nsclr ! number of sclr transport equations to solve 220 221 ndofBC = ndof + I3nsd ! dimension of BC array 222 ndiBCB = 2 ! dimension of iBCB array 223 ndBCB = ndof + 1 ! dimension of BCB array 224c 225 nsymdf = (ndof*(ndof + 1)) / 2 ! symm. d.o.f.'s 226c 227c.... ----------------------> Communication tasks <-------------------- 228c 229 230cc still read in new 231 232 if(numpe > 1) then 233 234cc MR CHANGE 235 call phio_readheader(igeom,'size of ilwork array' // char(0), 236 & nlwork,ione,'integer' // char(0) ,iotype) 237 238 call phio_readheader(igeom,'ilwork' //char(0) ,nlwork,ione, 239 & 'integer' // char(0) , iotype) 240 241 allocate( point2ilwork(nlwork) ) 242 allocate( ilworkread(nlwork) ) 243 call phio_readdatablock(igeom,'ilwork' // char(0),ilworkread, 244 & nlwork,'integer' // char(0) , iotype) 245 246c call closefile( igeom, "read" ) 247c call finalizephmpiio( igeom ) 248 249c fname1='size of ilwork array?' 250c call readheader(igeomBAK,fname1,nlwork,ione,'integer', iotype) 251 252c ione=1 253c fname1='ilwork?' 254c call readheader(igeomBAK,fname1,nlwork,ione,'integer', iotype) 255 256c allocate( point2ilwork(nlwork) ) 257c allocate( ilworkread(nlwork) ) 258c call readdatablock(igeomBAK,fname1,ilworkread, 259c & nlwork,'integer', iotype) 260cc MR CHANGE 261 262 263 point2ilwork = ilworkread 264 call ctypes (point2ilwork) 265 else 266 nlwork=1 267 allocate( point2ilwork(1)) 268 nshg0 = nshg 269 endif 270 271cccccccccccccccccccccccccccccccccccccccccc 272 273 itwo=2 274 275c print *, "fname2 is", fname2 276 277cc MR CHANGE 278c call initphmpiio(nfields,nppf,nfiles,igeom,'read') 279c call openfile( fnamer, 'read', igeom ) 280CC MR CHANGE 281 282 call phio_readheader(igeom,'co-ordinates' // char(0),intfromfile,itwo, 283 & 'double' // char(0), iotype) 284 numnp=intfromfile(1) 285c print *, "read out @@@@@@ is ", numnp 286 allocate( point2x(numnp,nsd) ) 287 allocate( xread(numnp,nsd) ) 288 ixsiz=numnp*nsd 289 call phio_readdatablock(igeom,'co-ordinates' // char(0),xread,ixsiz, 290 & 'double' // char(0), iotype) 291 point2x = xread 292 293! call closefile( igeom, "read" ) 294! call finalizephmpiio( igeom ) 295 296! print *, "Finished finalize" 297 298c deallocate (point2x) 299c deallocate (xread) 300 301cccccccccccccccccccccccccccccccccccccccccc 302 303c 304c.... read the node coordinates 305c 306 307c itwo=2 308c fname1='co-ordinates?' 309c call readheader(igeomBAK,fname1,intfromfile,itwo, 'double', iotype) 310c numnp=intfromfile(1) 311cc nsd=intfromfile(2) 312c allocate( point2x(numnp,nsd) ) 313c allocate( xread(numnp,nsd) ) 314c ixsiz=numnp*nsd 315c call readdatablock(igeomBAK,fname1,xread,ixsiz, 'double',iotype) 316c point2x = xread 317 318c 319c.... read in and block out the connectivity 320c 321 322! !MR CHANGE 323! This is not necessary but this avoids to have the geombc files opend two times. 324! A better way consists in pasisng the file handler to genblk or make it global or use igeomBAK instead of igeom 325! call closefile( igeom, "read" ) 326! call finalizephmpiio( igeom ) 327! !MR CHANGE END 328 329 call genblk (IBKSIZ) 330 331! !MR CHANGE 332! call initphmpiio( nfields, nppf, nfiles, igeom, 'read') 333! call openfile( fnamer, 'read', igeom ) 334! !MR CHANGE END 335 336c 337c.... read the boundary condition mapping array 338c 339 340cc MR CHANGE 341! call initphmpiio(nfields,nppf,nfiles,igeom, 'read') 342! call openfile( fnamer, 'read', igeom ) 343cc MR CHANGE 344 345 ione=1 346 call phio_readheader(igeom,'bc mapping array' // char(0),nshg,ione, 347 & 'integer' // char(0),iotype) 348 349c fname1='bc mapping array?' 350c call readheader(igeomBAK,fname1,nshg, 351c & ione,'integer', iotype) 352 353 allocate( nBC(nshg) ) 354 355 allocate( nBCread(nshg) ) 356 357c call readdatablock(igeomBAK,fname1,nBCread,nshg,'integer',iotype) 358 call phio_readdatablock(igeom,'bc mapping array' // char(0), 359 & nBCread,nshg,'integer' // char(0),iotype) 360 361 nBC=nBCread 362 363c 364c.... read the temporary iBC array 365c 366 ione=1 367 call phio_readheader(igeom,'bc codes array' // char(0) ,numpbc,ione, 368 & 'integer' // char(0),iotype) 369 370c ione = 1 371c fname1='bc codes array?' 372c call readheader(igeomBAK,fname1,numpbc, 373c & ione, 'integer', iotype) 374 375!MR CHANGE 376! if ( numpbc > 0 ) then 377! allocate( iBCtmp(numpbc) ) 378! allocate( iBCtmpread(numpbc) ) 379! c call readdatablock(igeomBAK,fname1,iBCtmpread,numpbc, 380! c & 'integer',iotype) 381! call readdatablock(igeom,fname2,iBCtmpread,numpbc, 382! & 'integer',iotype) 383! iBCtmp=iBCtmpread 384! else ! sometimes a partition has no BC's 385! allocate( iBCtmp(1) ) 386! iBCtmp=0 387! endif 388 389 if ( numpbc > 0 ) then 390 allocate( iBCtmp(numpbc) ) 391 allocate( iBCtmpread(numpbc) ) 392 else 393 allocate( iBCtmp(1) ) 394 allocate( iBCtmpread(1) ) 395 endif 396c call readdatablock(igeomBAK,fname1,iBCtmpread,numpbc, 397c & 'integer',iotype) 398 call phio_readdatablock(igeom,'bc codes array' // char(0), 399 & iBCtmpread,numpbc,'integer' // char(0),iotype) 400 401 if ( numpbc > 0 ) then 402 iBCtmp=iBCtmpread 403 else ! sometimes a partition has no BC's 404 deallocate( iBCtmpread) 405 iBCtmp=0 406 endif 407!MR CHANGE END 408 409c 410c.... read boundary condition data 411c 412 413 ione=1 414 415c ione=1 416c fname1='boundary condition array?' 417c call readheader(igeomBAK,fname1,intfromfile, 418c & ione, 'double', iotype) 419 call phio_readheader(igeom,'boundary condition array' // char(0), 420 & intfromfile,ione, 'double' // char(0), iotype) 421c here intfromfile(1) contains (ndof+7)*numpbc 422!MR CHANGE 423! if ( numpbc > 0 ) then 424! if(intfromfile(1).ne.(ndof+7)*numpbc) then 425! warning='WARNING more data in BCinp than needed: keeping 1st' 426! write(*,*) warning, ndof+7 427! endif 428! allocate( BCinp(numpbc,ndof+7) ) 429! nsecondrank=intfromfile(1)/numpbc 430! allocate( BCinpread(numpbc,nsecondrank) ) 431! iBCinpsiz=intfromfile(1) 432! c call readdatablock(igeomBAK,fname1,BCinpread,iBCinpsiz, 433! c & 'double',iotype) 434! call readdatablock(igeom,fname2,BCinpread,iBCinpsiz, 435! & 'double',iotype) 436! BCinp(:,1:(ndof+7))=BCinpread(:,1:(ndof+7)) 437! else ! sometimes a partition has no BC's 438! allocate( BCinp(1,ndof+7) ) 439! BCinp=0 440! endif 441 442 if ( numpbc > 0 ) then 443! if(intfromfile(1).ne.(ndof+7)*numpbc) then 444! warning='WARNING more data in BCinp than needed: keeping 1st' 445! write(*,*) warning, ndof+7 446! endif 447 allocate( BCinp(numpbc,ndof+7) ) 448 nsecondrank=intfromfile(1)/numpbc 449 allocate( BCinpread(numpbc,nsecondrank) ) 450 iBCinpsiz=intfromfile(1) 451 else 452 allocate( BCinp(1,ndof+7) ) 453 allocate( BCinpread(0,0) ) !dummy 454 iBCinpsiz=intfromfile(1) 455 endif 456c call readdatablock(igeomBAK,fname1,BCinpread,iBCinpsiz, 457c & 'double',iotype) 458 459 call phio_readdatablock(igeom,'boundary condition array' // char(0), 460 & BCinpread,iBCinpsiz,'double' // char(0) ,iotype) 461 462 if ( numpbc > 0 ) then 463 BCinp(:,1:(ndof+7))=BCinpread(:,1:(ndof+7)) 464 else ! sometimes a partition has no BC's 465 deallocate(BCinpread) 466 BCinp=0 467 endif 468!MR CHANGE END 469 470c 471c.... read periodic boundary conditions 472c 473 474 ione=1 475c fname1='periodic masters array?' 476c call readheader(igeomBAK,fname1,nshg, 477c & ione, 'integer', iotype) 478 call phio_readheader(igeom,'periodic masters array' // char(0) ,nshg, 479 & ione,'integer' // char(0), iotype) 480 allocate( point2iper(nshg) ) 481 allocate( iperread(nshg) ) 482c call readdatablock(igeomBAK,fname1,iperread,nshg, 483c & 'integer',iotype) 484 call phio_readdatablock(igeom,'periodic masters array' // char(0), 485 & iperread,nshg,'integer' // char(0),iotype) 486 point2iper=iperread 487 488 489! !MR CHANGE 490! call closefile( igeom, "read" ) 491! call finalizephmpiio( igeom ) 492! !MR CHANGE END 493 494c 495c.... generate the boundary element blocks 496c 497 call genbkb (ibksiz) 498 499 500! !MR CHANGE 501! write(*,*) 'HELLO 12 from ', myrank 502! !MR CHANGE END 503 504c 505c Read in the nsons and ifath arrays if needed 506c 507c There is a fundamental shift in the meaning of ifath based on whether 508c there exist homogenous directions in the flow. 509c 510c HOMOGENOUS DIRECTIONS EXIST: Here nfath is the number of inhomogenous 511c points in the TOTAL mesh. That is to say that each partition keeps a 512c link to ALL inhomogenous points. This link is furthermore not to the 513c sms numbering but to the original structured grid numbering. These 514c inhomogenous points are thought of as fathers, with their sons being all 515c the points in the homogenous directions that have this father's 516c inhomogeneity. The array ifath takes as an arguement the sms numbering 517c and returns as a result the father. 518c 519c In this case nsons is the number of sons that each father has and ifath 520c is an array which tells the 521c 522c NO HOMOGENOUS DIRECTIONS. In this case the mesh would grow to rapidly 523c if we followed the above strategy since every partition would index its 524c points to the ENTIRE mesh. Furthermore, there would never be a need 525c to average to a node off processor since there is no spatial averaging. 526c Therefore, to properly account for this case we must recognize it and 527c inerrupt certain actions (i.e. assembly of the average across partitions). 528c This case is easily identified by noting that maxval(nsons) =1 (i.e. no 529c father has any sons). Reiterating to be clear, in this case ifath does 530c not point to a global numbering but instead just points to itself. 531c 532 nfath=1 ! some architectures choke on a zero or undeclared 533 ! dimension variable. This sets it to a safe, small value. 534 if(((iLES .lt. 20) .and. (iLES.gt.0)) 535 & .or. (itwmod.gt.0) ) then ! don't forget same 536 ! conditional in proces.f 537 538c read (igeomBAK) nfath ! nfath already read in input.f, 539 ! needed for alloc 540 ione=1 541c call creadlist(igeomBAK,ione,nfath) 542c fname1='keyword sonfath?' 543 if(nohomog.gt.0) then 544 545 546! fname1='number of father-nodes?' 547! call readheader(igeomBAK,fname1,nfath,ione,'integer', iotype) 548 549 call phio_readheader(igeom,'number of father-nodes' // char(0), 550 & nfath,ione,'integer' // char(0), iotype) 551 552c 553c fname1='keyword nsons?' 554! fname1='number of son-nodes for each father?' 555! call readheader(igeomBAK,fname1,nfath,ione,'integer', iotype) 556 557 call phio_readheader(igeom, 558 & 'number of son-nodes for each father' // char(0), 559 & nfath,ione,'integer' // char(0), iotype) 560 561 allocate (point2nsons(nfath)) 562 563! call readdatablock(igeomBAK,fname1,point2nsons,nfath, 564! & 'integer',iotype) 565 call phio_readdatablock(igeom, 566 & 'number of son-nodes for each father' // char(0), 567 & point2nsons,nfath,'integer' // char(0), iotype) 568 569c 570! fname1='keyword ifath?' 571! call readheader(igeomBAK,fname1,nshg,ione,'integer', iotype) 572 573 call phio_readheader(igeom,'keyword ifath' // char(0),nshg,ione, 574 & 'integer' // char(0), iotype) 575 576 allocate (point2ifath(nshg)) 577 578! call readdatablock(igeomBAK,fname1,point2ifath,nshg, 579! & 'integer',iotype) 580 call phio_readdatablock(igeom,'keyword ifath' // char(0),point2ifath, 581 & nshg,'integer' // char(0) , iotype) 582 583c 584 nsonmax=maxval(point2nsons) 585c 586 else ! this is the case where there is no homogeneity 587 ! therefore ever node is a father (too itself). sonfath 588 ! (a routine in NSpre) will set this up but this gives 589 ! you an option to avoid that. 590 nfath=nshg 591 allocate (point2nsons(nfath)) 592 point2nsons=1 593 allocate (point2ifath(nshg)) 594 do i=1,nshg 595 point2ifath(i)=i 596 enddo 597 nsonmax=1 598c 599 endif 600 else 601 allocate (point2nsons(1)) 602 allocate (point2ifath(1)) 603 endif 604 605 call phio_closefile_read(igeom); 606 607! !MR CHANGE 608! write(*,*) 'HELLO 13 from ', myrank 609! !MR CHANGE END 610 611c 612c renumber the master partition for SPEBC 613c 614c if((myrank.eq.master).and.(irscale.ge.0)) then 615c call setSPEBC(numnp, nfath, nsonmax) 616c call renum(point2x,point2ifath,point2nsons) 617c endif 618c 619c.... Read restart files 620 621c$$$ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 622c 623c nfields = 11 624c numparts = 512 625c nppp = numparts/numpe 626c startpart = myrank * nppp +1 627c int endpart = startpart + nppp - 1 628c nppf = numparts/nfiles 629cc fnamer = "/users/nliu/PIG4/4-procs_case/restart-file" 630cc fnamer="./4-procs_case/restart-file" 631 632 call phio_restartname(irstart, fnamer) 633 call phio_openfile_read(fnamer, nfiles, descriptor) 634 635 ithree=3 636c call creadlist(irstin,ithree,nshg2,ndof2,lstep) 637 638 itmp = int(log10(float(myrank+1)))+1 639 640c print *, "Solution is : ", fname1 641 642 intfromfile=0 643 call phio_readheader(descriptor,'solution' // char(0) ,intfromfile, 644 & ithree,'integer' // char(0), iotype) 645c 646c.... read the values of primitive variables into q 647c 648 649c print *, "intfromfile(1) is ", intfromfile(1) 650c print *, "intfromfile(2) is ", intfromfile(2) 651c print *, "intfromfile(3) is ", intfromfile(3) 652 653 allocate( qold(nshg,ndof) ) 654 if(intfromfile(1).ne.0) then 655 nshg2=intfromfile(1) 656 ndof2=intfromfile(2) 657 lstep=intfromfile(3) 658 if(ndof2.ne.ndof) then 659 660 endif 661c 662 if (nshg2 .ne. nshg) 663 & call error ('restar ', 'nshg ', nshg) 664 allocate( qread(nshg,ndof2) ) 665 iqsiz=nshg*ndof2 666 call phio_readdatablock(descriptor,'solution' // char(0),qread,iqsiz, 667 & 'double' // char(0),iotype) 668 qold(:,1:ndof)=qread(:,1:ndof) 669 deallocate(qread) 670 else 671 if (myrank.eq.master) then 672 if (matflg(1,1).eq.0) then ! compressible 673 warning='Solution is set to zero (with p and T to one)' 674 else 675 warning='Solution is set to zero' 676 endif 677 write(*,*) warning 678 endif 679 qold=zero 680 if (matflg(1,1).eq.0) then ! compressible 681 qold(:,1)=one ! avoid zero pressure 682 qold(:,nflow)=one ! avoid zero temperature 683 endif 684 endif 685 686 687! !MR CHANGE 688! write(*,*) 'HELLO 16-8 from ', myrank 689! !MR CHANGE END 690 691c itmp=1 692c if (myrank .gt. 0) itmp = int(log10(float(myrank)))+1 693 intfromfile=0 694 call phio_readheader(descriptor,'time derivative of solution' // char(0), 695 & intfromfile,ithree,'integer' // char(0),iotype) 696 allocate( acold(nshg,ndof) ) 697 if(intfromfile(1).ne.0) then 698 nshg2=intfromfile(1) 699 ndof2=intfromfile(2) 700 lstep=intfromfile(3) 701 702c print *, "intfromfile(1) is ", intfromfile(1) 703c print *, "intfromfile(2) is ", intfromfile(2) 704c print *, "intfromfile(3) is ", intfromfile(3) 705 706 if (nshg2 .ne. nshg) 707 & call error ('restar ', 'nshg ', nshg) 708c 709 allocate( acread(nshg,ndof2) ) 710 acread=zero 711 iacsiz=nshg*ndof2 712 call phio_readdatablock(descriptor, 713 & 'time derivative of solution' // char(0),acread, 714 & iacsiz, 'double' // char(0),iotype) 715 acold(:,1:ndof)=acread(:,1:ndof) 716 deallocate(acread) 717 else 718 if (myrank.eq.master) then 719 warning='Time derivative of solution is set to zero (SAFE)' 720 write(*,*) warning 721 endif 722 acold=zero 723 endif 724 725cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 726cc 727c 728cc 729cc.... read the header and check it against the run data 730cc 731 732 733c ithree=3 734ccc call creadlist(irstin,ithree,nshg2,ndof2,lstep) 735c fname1='solution?' 736c intfromfile=0 737c call readheader(irstin,fname1,intfromfile, 738c & ithree,'integer', iotype) 739cc 740cc.... read the values of primitive variables into q 741cc 742c allocate( qold(nshg,ndof) ) 743c if(intfromfile(1).ne.0) then 744c nshg2=intfromfile(1) 745c ndof2=intfromfile(2) 746c lstep=intfromfile(3) 747c if(ndof2.ne.ndof) then 748c warning='WARNING more data in restart than needed: keeping 1st ' 749c write(*,*) warning , ndof 750c endif 751cc 752c if (nshg2 .ne. nshg) 753c & call error ('restar ', 'nshg ', nshg) 754c allocate( qread(nshg,ndof2) ) 755 756c iqsiz=nshg*ndof2 757c call readdatablock(irstin,fname1,qread,iqsiz, 758c & 'double',iotype) 759c qold(:,1:ndof)=qread(:,1:ndof) 760c deallocate(qread) 761c else 762c if (myrank.eq.master) then 763c if (matflg(1,1).eq.0) then ! compressible 764c warning='Solution is set to zero (with p and T to one)' 765c else 766c warning='Solution is set to zero' 767c endif 768c write(*,*) warning 769c endif 770c qold=zero 771c if (matflg(1,1).eq.0) then ! compressible 772c qold(:,1)=one ! avoid zero pressure 773c qold(:,nflow)=one ! avoid zero temperature 774c endif 775c endif 776cc 777c fname1='time derivative of solution?' 778c intfromfile=0 779c call readheader(irstin,fname1,intfromfile, 780c & ithree,'integer', iotype) 781c allocate( acold(nshg,ndof) ) 782c if(intfromfile(1).ne.0) then 783c nshg2=intfromfile(1) 784c ndof2=intfromfile(2) 785c lstep=intfromfile(3) 786c 787c if (nshg2 .ne. nshg) 788c & call error ('restar ', 'nshg ', nshg) 789cc 790c allocate( acread(nshg,ndof2) ) 791c acread=zero 792c 793c iacsiz=nshg*ndof2 794c call readdatablock(irstin,fname1,acread,iacsiz, 795c & 'double',iotype) 796c acold(:,1:ndof)=acread(:,1:ndof) 797c deallocate(acread) 798c else 799c if (myrank.eq.master) then 800c warning='Time derivative of solution is set to zero (SAFE)' 801c write(*,*) warning 802c endif 803c acold=zero 804c endif 805 806c call creadlist(irstin,ithree,nshg2,ndisp,lstep) 807 808 if (ideformwall.eq.1) then 809! fname1='displacement?' 810! call readheader(irstin,fname1,intfromfile, 811! & ithree,'integer', iotype) 812 813 intfromfile=0 814 call phio_readheader(descriptor,'displacement' // char(0), 815 & intfromfile,ithree,'integer' // char(0),iotype) 816 817 nshg2=intfromfile(1) 818 ndisp=intfromfile(2) 819 lstep=intfromfile(3) 820 if(ndisp.ne.nsd) then 821 warning='WARNING ndisp not equal nsd' 822 write(*,*) warning , ndisp 823 endif 824c 825 if (nshg2 .ne. nshg) 826 & call error ('restar ', 'nshg ', nshg) 827c 828c.... read the values of primitive variables into uold 829c 830 831 allocate( uold(nshg,nsd) ) 832 allocate( uread(nshg,nsd) ) 833 834 iusiz=nshg*nsd 835 836! call readdatablock(irstin,fname1,uread,iusiz, 837! & 'double',iotype) 838 call phio_readdatablock(descriptor,'displacement' // char(0), 839 & uread,iusiz, 'double' // char(0),iotype) 840 841 uold(:,1:nsd)=uread(:,1:nsd) 842 else 843 allocate( uold(nshg,nsd) ) 844 uold(:,1:nsd) = zero 845 endif 846 847c 848c.... close c-binary files 849c 850!MR CHANGE 851! call closefile( irstin, "read" ) 852 853 call phio_closefile_read(descriptor) 854 855!MR CHANGE 856! call closefile( igeomBAK, "read" ) 857c 858 deallocate(xread) 859 if ( numpbc > 0 ) then 860 deallocate(bcinpread) 861 deallocate(ibctmpread) 862 endif 863 deallocate(iperread) 864 if(numpe.gt.1) 865 & deallocate(ilworkread) 866 deallocate(nbcread) 867 868 return 869c 870 994 call error ('input ','opening ', igeomBAK) 871 995 call error ('input ','opening ', igeomBAK) 872 997 call error ('input ','end file', igeomBAK) 873 998 call error ('input ','end file', igeomBAK) 874c 875 end 876 877c 878c No longer called but kept around in case.... 879c 880 subroutine genpzero(iBC) 881 882 use pointer_data 883c 884 include "common.h" 885 integer iBC(nshg) 886c 887c.... check to see if any of the nodes have a dirichlet pressure 888c 889 pzero=1 890 if (any(btest(iBC,2))) pzero=0 891c 892 do iblk = 1, nelblb 893 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 894 do i=1, npro 895 iBCB1=miBCB(iblk)%p(i,1) 896c 897c.... check to see if any of the nodes have a Neumann pressure 898c but not periodic (note that 899c 900 if(btest(iBCB1,1)) pzero=0 901 enddo 902c 903c.... share results with other processors 904c 905 pzl=pzero 906 if (numpe .gt. 1) 907 & call MPI_ALLREDUCE (pzl, pzero, 1, 908 & MPI_DOUBLE_PRECISION,MPI_MIN, MPI_COMM_WORLD,ierr) 909 910 enddo 911c 912c.... return 913c 914 return 915c 916 end 917 918