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