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 iso_c_binding 35 use readarrays 36 use phio 37 include "common.h" 38c 39 real*8, allocatable :: xread(:,:), qread(:,:), acread(:,:) 40 real*8, allocatable :: uread(:,:) 41 real*8, allocatable :: BCinpread(:,:) 42 integer, allocatable :: iperread(:), iBCtmpread(:) 43 integer, 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 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, handle); 139 140 call phio_readheader(handle,c_char_'number of nodes' // char(0), 141 & c_loc(numnp),ione, dataInt, iotype) 142 143 call phio_readheader(handle,c_char_'number of modes' // char(0), 144 & c_loc(nshg),ione, dataInt, iotype) 145 146 call phio_readheader(handle, 147 & c_char_'number of interior elements' // char(0), 148 & c_loc(numel),ione, dataInt, iotype) 149 150 call phio_readheader(handle, 151 & c_char_'number of boundary elements' // char(0), 152 & c_loc(numelb),ione, dataInt, iotype) 153 154 call phio_readheader(handle, 155 & c_char_'maximum number of element nodes' // char(0), 156 & c_loc(nen),ione, dataInt, iotype) 157 158 call phio_readheader(handle, 159 & c_char_'number of interior tpblocks' // char(0), 160 & c_loc(nelblk),ione, dataInt, iotype) 161 162 call phio_readheader(handle, 163 & c_char_'number of boundary tpblocks' // char(0), 164 & c_loc(nelblb),ione, dataInt, iotype) 165 166 call phio_readheader(handle, 167 & c_char_'number of nodes with Dirichlet BCs' // char(0), 168 & c_loc(numpbc),ione, dataInt, iotype) 169 170 call phio_readheader(handle, 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(handle, 248 & c_char_'size of ilwork array' // char(0), 249 & c_loc(nlwork),ione, dataInt, iotype) 250 251 call phio_readheader(handle, 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(igeom,'ilwork' // char(0),ilworkread, 258 & nlwork,'integer' // char(0) , 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(handle, 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(igeom,'co-ordinates' // char(0),xread,ixsiz, 305 & 'double' // char(0), iotype) 306 point2x = xread 307 308! call closefile( igeom, "read" ) 309! call finalizephmpiio( igeom ) 310 311! print *, "Finished finalize" 312 313c deallocate (point2x) 314c deallocate (xread) 315 316cccccccccccccccccccccccccccccccccccccccccc 317 318c 319c.... read the node coordinates 320c 321 322c itwo=2 323c fname1='co-ordinates?' 324c call readheader(igeomBAK,fname1,intfromfile,itwo, 'double', iotype) 325c numnp=intfromfile(1) 326cc nsd=intfromfile(2) 327c allocate( point2x(numnp,nsd) ) 328c allocate( xread(numnp,nsd) ) 329c ixsiz=numnp*nsd 330c call readdatablock(igeomBAK,fname1,xread,ixsiz, 'double',iotype) 331c point2x = xread 332 333c 334c.... read in and block out the connectivity 335c 336 337! !MR CHANGE 338! This is not necessary but this avoids to have the geombc files opend two times. 339! A better way consists in pasisng the file handler to genblk or make it global or use igeomBAK instead of igeom 340! call closefile( igeom, "read" ) 341! call finalizephmpiio( igeom ) 342! !MR CHANGE END 343 344 call genblk (IBKSIZ) 345 346! !MR CHANGE 347! call initphmpiio( nfields, nppf, nfiles, igeom, 'read') 348! call openfile( fnamer, 'read', igeom ) 349! !MR CHANGE END 350 351c 352c.... read the boundary condition mapping array 353c 354 355cc MR CHANGE 356! call initphmpiio(nfields,nppf,nfiles,igeom, 'read') 357! call openfile( fnamer, 'read', igeom ) 358cc MR CHANGE 359 360 ione=1 361 call phio_readheader(handle, 362 & c_char_'bc mapping array' // char(0), 363 & c_loc(nshg),ione, dataInt, iotype) 364 365c fname1='bc mapping array?' 366c call readheader(igeomBAK,fname1,nshg, 367c & ione,'integer', iotype) 368 369 allocate( nBC(nshg) ) 370 371 allocate( nBCread(nshg) ) 372 373c call readdatablock(igeomBAK,fname1,nBCread,nshg,'integer',iotype) 374 call phio_readdatablock(igeom,'bc mapping array' // char(0), 375 & nBCread,nshg,'integer' // char(0),iotype) 376 377 nBC=nBCread 378 379c 380c.... read the temporary iBC array 381c 382 ione=1 383 call phio_readheader(handle, 384 & c_char_'bc codes array' // char(0), 385 & c_loc(numpbc),ione, dataInt, iotype) 386 387c ione = 1 388c fname1='bc codes array?' 389c call readheader(igeomBAK,fname1,numpbc, 390c & ione, 'integer', iotype) 391 392!MR CHANGE 393! if ( numpbc > 0 ) then 394! allocate( iBCtmp(numpbc) ) 395! allocate( iBCtmpread(numpbc) ) 396! c call readdatablock(igeomBAK,fname1,iBCtmpread,numpbc, 397! c & 'integer',iotype) 398! call readdatablock(igeom,fname2,iBCtmpread,numpbc, 399! & 'integer',iotype) 400! iBCtmp=iBCtmpread 401! else ! sometimes a partition has no BC's 402! allocate( iBCtmp(1) ) 403! iBCtmp=0 404! endif 405 406 if ( numpbc > 0 ) then 407 allocate( iBCtmp(numpbc) ) 408 allocate( iBCtmpread(numpbc) ) 409 else 410 allocate( iBCtmp(1) ) 411 allocate( iBCtmpread(1) ) 412 endif 413c call readdatablock(igeomBAK,fname1,iBCtmpread,numpbc, 414c & 'integer',iotype) 415 call phio_readdatablock(igeom,'bc codes array' // char(0), 416 & iBCtmpread,numpbc,'integer' // char(0),iotype) 417 418 if ( numpbc > 0 ) then 419 iBCtmp=iBCtmpread 420 else ! sometimes a partition has no BC's 421 deallocate( iBCtmpread) 422 iBCtmp=0 423 endif 424!MR CHANGE END 425 426c 427c.... read boundary condition data 428c 429 430 ione=1 431 432c ione=1 433c fname1='boundary condition array?' 434c call readheader(igeomBAK,fname1,intfromfile, 435c & ione, 'double', iotype) 436 call phio_readheader(handle, 437 & c_char_'boundary condition array' // char(0), 438 & c_loc(intfromfile),ione, dataDbl, iotype) 439c here intfromfile(1) contains (ndof+7)*numpbc 440!MR CHANGE 441! if ( numpbc > 0 ) then 442! if(intfromfile(1).ne.(ndof+7)*numpbc) then 443! warning='WARNING more data in BCinp than needed: keeping 1st' 444! write(*,*) warning, ndof+7 445! endif 446! allocate( BCinp(numpbc,ndof+7) ) 447! nsecondrank=intfromfile(1)/numpbc 448! allocate( BCinpread(numpbc,nsecondrank) ) 449! iBCinpsiz=intfromfile(1) 450! c call readdatablock(igeomBAK,fname1,BCinpread,iBCinpsiz, 451! c & 'double',iotype) 452! call readdatablock(igeom,fname2,BCinpread,iBCinpsiz, 453! & 'double',iotype) 454! BCinp(:,1:(ndof+7))=BCinpread(:,1:(ndof+7)) 455! else ! sometimes a partition has no BC's 456! allocate( BCinp(1,ndof+7) ) 457! BCinp=0 458! endif 459 460 if ( numpbc > 0 ) then 461! if(intfromfile(1).ne.(ndof+7)*numpbc) then 462! warning='WARNING more data in BCinp than needed: keeping 1st' 463! write(*,*) warning, ndof+7 464! endif 465 allocate( BCinp(numpbc,ndof+7) ) 466 nsecondrank=intfromfile(1)/numpbc 467 allocate( BCinpread(numpbc,nsecondrank) ) 468 iBCinpsiz=intfromfile(1) 469 else 470 allocate( BCinp(1,ndof+7) ) 471 allocate( BCinpread(0,0) ) !dummy 472 iBCinpsiz=intfromfile(1) 473 endif 474c call readdatablock(igeomBAK,fname1,BCinpread,iBCinpsiz, 475c & 'double',iotype) 476 477 call phio_readdatablock(igeom,'boundary condition array' // char(0), 478 & BCinpread,iBCinpsiz,'double' // char(0) ,iotype) 479 480 if ( numpbc > 0 ) then 481 BCinp(:,1:(ndof+7))=BCinpread(:,1:(ndof+7)) 482 else ! sometimes a partition has no BC's 483 deallocate(BCinpread) 484 BCinp=0 485 endif 486!MR CHANGE END 487 488c 489c.... read periodic boundary conditions 490c 491 492 ione=1 493c fname1='periodic masters array?' 494c call readheader(igeomBAK,fname1,nshg, 495c & ione, 'integer', iotype) 496 call phio_readheader(handle, 497 & c_char_'periodic masters array' // char(0), 498 & c_loc(nshg), ione, dataInt, iotype) 499 allocate( point2iper(nshg) ) 500 allocate( iperread(nshg) ) 501c call readdatablock(igeomBAK,fname1,iperread,nshg, 502c & 'integer',iotype) 503 call phio_readdatablock(igeom,'periodic masters array' // char(0), 504 & iperread,nshg,'integer' // char(0),iotype) 505 point2iper=iperread 506 507 508! !MR CHANGE 509! call closefile( igeom, "read" ) 510! call finalizephmpiio( igeom ) 511! !MR CHANGE END 512 513c 514c.... generate the boundary element blocks 515c 516 call genbkb (ibksiz) 517 518 519! !MR CHANGE 520! write(*,*) 'HELLO 12 from ', myrank 521! !MR CHANGE END 522 523c 524c Read in the nsons and ifath arrays if needed 525c 526c There is a fundamental shift in the meaning of ifath based on whether 527c there exist homogenous directions in the flow. 528c 529c HOMOGENOUS DIRECTIONS EXIST: Here nfath is the number of inhomogenous 530c points in the TOTAL mesh. That is to say that each partition keeps a 531c link to ALL inhomogenous points. This link is furthermore not to the 532c sms numbering but to the original structured grid numbering. These 533c inhomogenous points are thought of as fathers, with their sons being all 534c the points in the homogenous directions that have this father's 535c inhomogeneity. The array ifath takes as an arguement the sms numbering 536c and returns as a result the father. 537c 538c In this case nsons is the number of sons that each father has and ifath 539c is an array which tells the 540c 541c NO HOMOGENOUS DIRECTIONS. In this case the mesh would grow to rapidly 542c if we followed the above strategy since every partition would index its 543c points to the ENTIRE mesh. Furthermore, there would never be a need 544c to average to a node off processor since there is no spatial averaging. 545c Therefore, to properly account for this case we must recognize it and 546c inerrupt certain actions (i.e. assembly of the average across partitions). 547c This case is easily identified by noting that maxval(nsons) =1 (i.e. no 548c father has any sons). Reiterating to be clear, in this case ifath does 549c not point to a global numbering but instead just points to itself. 550c 551 nfath=1 ! some architectures choke on a zero or undeclared 552 ! dimension variable. This sets it to a safe, small value. 553 if(((iLES .lt. 20) .and. (iLES.gt.0)) 554 & .or. (itwmod.gt.0) ) then ! don't forget same 555 ! conditional in proces.f 556 557c read (igeomBAK) nfath ! nfath already read in input.f, 558 ! needed for alloc 559 ione=1 560c call creadlist(igeomBAK,ione,nfath) 561c fname1='keyword sonfath?' 562 if(nohomog.gt.0) then 563 564 565! fname1='number of father-nodes?' 566! call readheader(igeomBAK,fname1,nfath,ione,'integer', iotype) 567 568 call phio_readheader(handle, 569 & c_char_'number of father-nodes' // char(0), 570 & c_loc(nfath), ione, dataInt, iotype) 571 572c 573c fname1='keyword nsons?' 574! fname1='number of son-nodes for each father?' 575! call readheader(igeomBAK,fname1,nfath,ione,'integer', iotype) 576 577 call phio_readheader(handle, 578 & c_char_'number of son-nodes for each father' // char(0), 579 & c_loc(nfath), ione, dataInt, iotype) 580 581 allocate (point2nsons(nfath)) 582 583! call readdatablock(igeomBAK,fname1,point2nsons,nfath, 584! & 'integer',iotype) 585 call phio_readdatablock(igeom, 586 & 'number of son-nodes for each father' // char(0), 587 & point2nsons,nfath,'integer' // char(0), iotype) 588 589c 590! fname1='keyword ifath?' 591! call readheader(igeomBAK,fname1,nshg,ione,'integer', iotype) 592 593 call phio_readheader(handle, 594 & c_char_'keyword ifath' // char(0), 595 & c_loc(nshg), ione, dataInt, iotype); 596 597 allocate (point2ifath(nshg)) 598 599! call readdatablock(igeomBAK,fname1,point2ifath,nshg, 600! & 'integer',iotype) 601 call phio_readdatablock(igeom,'keyword ifath' // char(0),point2ifath, 602 & nshg,'integer' // char(0) , iotype) 603 604c 605 nsonmax=maxval(point2nsons) 606c 607 else ! this is the case where there is no homogeneity 608 ! therefore ever node is a father (too itself). sonfath 609 ! (a routine in NSpre) will set this up but this gives 610 ! you an option to avoid that. 611 nfath=nshg 612 allocate (point2nsons(nfath)) 613 point2nsons=1 614 allocate (point2ifath(nshg)) 615 do i=1,nshg 616 point2ifath(i)=i 617 enddo 618 nsonmax=1 619c 620 endif 621 else 622 allocate (point2nsons(1)) 623 allocate (point2ifath(1)) 624 endif 625 626 call phio_closefile_read(handle); 627 628! !MR CHANGE 629! write(*,*) 'HELLO 13 from ', myrank 630! !MR CHANGE END 631 632c 633c renumber the master partition for SPEBC 634c 635c if((myrank.eq.master).and.(irscale.ge.0)) then 636c call setSPEBC(numnp, nfath, nsonmax) 637c call renum(point2x,point2ifath,point2nsons) 638c endif 639c 640c.... Read restart files 641 642c$$$ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 643c 644c nfields = 11 645c numparts = 512 646c nppp = numparts/numpe 647c startpart = myrank * nppp +1 648c int endpart = startpart + nppp - 1 649c nppf = numparts/nfiles 650cc fnamer = "/users/nliu/PIG4/4-procs_case/restart-file" 651cc fnamer="./4-procs_case/restart-file" 652 653 call phio_restartname(irstart, fnamer) 654 call phio_openfile_read(fnamer, nfiles, handle) 655 656 ithree=3 657c call creadlist(irstin,ithree,nshg2,ndof2,lstep) 658 659 itmp = int(log10(float(myrank+1)))+1 660 661c print *, "Solution is : ", fname1 662 663 intfromfile=0 664 call phio_readheader(handle, 665 & c_char_'solution' // char(0), 666 & c_loc(intfromfile), ithree, dataInt, iotype) 667c 668c.... read the values of primitive variables into q 669c 670 671c print *, "intfromfile(1) is ", intfromfile(1) 672c print *, "intfromfile(2) is ", intfromfile(2) 673c print *, "intfromfile(3) is ", intfromfile(3) 674 675 allocate( qold(nshg,ndof) ) 676 if(intfromfile(1).ne.0) then 677 nshg2=intfromfile(1) 678 ndof2=intfromfile(2) 679 lstep=intfromfile(3) 680 if(ndof2.ne.ndof) then 681 682 endif 683c 684 if (nshg2 .ne. nshg) 685 & call error ('restar ', 'nshg ', nshg) 686 allocate( qread(nshg,ndof2) ) 687 iqsiz=nshg*ndof2 688 call phio_readdatablock(descriptor,'solution' // char(0),qread,iqsiz, 689 & 'double' // char(0),iotype) 690 qold(:,1:ndof)=qread(:,1:ndof) 691 deallocate(qread) 692 else 693 if (myrank.eq.master) then 694 if (matflg(1,1).eq.0) then ! compressible 695 warning='Solution is set to zero (with p and T to one)' 696 else 697 warning='Solution is set to zero' 698 endif 699 write(*,*) warning 700 endif 701 qold=zero 702 if (matflg(1,1).eq.0) then ! compressible 703 qold(:,1)=one ! avoid zero pressure 704 qold(:,nflow)=one ! avoid zero temperature 705 endif 706 endif 707 708 709! !MR CHANGE 710! write(*,*) 'HELLO 16-8 from ', myrank 711! !MR CHANGE END 712 713c itmp=1 714c if (myrank .gt. 0) itmp = int(log10(float(myrank)))+1 715 intfromfile=0 716 call phio_readheader(handle, 717 & c_char_'time derivative of solution' // char(0), 718 & c_loc(intfromfile), ithree, dataInt, iotype) 719 allocate( acold(nshg,ndof) ) 720 if(intfromfile(1).ne.0) then 721 nshg2=intfromfile(1) 722 ndof2=intfromfile(2) 723 lstep=intfromfile(3) 724 725c print *, "intfromfile(1) is ", intfromfile(1) 726c print *, "intfromfile(2) is ", intfromfile(2) 727c print *, "intfromfile(3) is ", intfromfile(3) 728 729 if (nshg2 .ne. nshg) 730 & call error ('restar ', 'nshg ', nshg) 731c 732 allocate( acread(nshg,ndof2) ) 733 acread=zero 734 iacsiz=nshg*ndof2 735 call phio_readdatablock(descriptor, 736 & 'time derivative of solution' // char(0),acread, 737 & iacsiz, 'double' // char(0),iotype) 738 acold(:,1:ndof)=acread(:,1:ndof) 739 deallocate(acread) 740 else 741 if (myrank.eq.master) then 742 warning='Time derivative of solution is set to zero (SAFE)' 743 write(*,*) warning 744 endif 745 acold=zero 746 endif 747 748cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 749cc 750c 751cc 752cc.... read the header and check it against the run data 753cc 754 755 756c ithree=3 757ccc call creadlist(irstin,ithree,nshg2,ndof2,lstep) 758c fname1='solution?' 759c intfromfile=0 760c call readheader(irstin,fname1,intfromfile, 761c & ithree,'integer', iotype) 762cc 763cc.... read the values of primitive variables into q 764cc 765c allocate( qold(nshg,ndof) ) 766c if(intfromfile(1).ne.0) then 767c nshg2=intfromfile(1) 768c ndof2=intfromfile(2) 769c lstep=intfromfile(3) 770c if(ndof2.ne.ndof) then 771c warning='WARNING more data in restart than needed: keeping 1st ' 772c write(*,*) warning , ndof 773c endif 774cc 775c if (nshg2 .ne. nshg) 776c & call error ('restar ', 'nshg ', nshg) 777c allocate( qread(nshg,ndof2) ) 778 779c iqsiz=nshg*ndof2 780c call readdatablock(irstin,fname1,qread,iqsiz, 781c & 'double',iotype) 782c qold(:,1:ndof)=qread(:,1:ndof) 783c deallocate(qread) 784c else 785c if (myrank.eq.master) then 786c if (matflg(1,1).eq.0) then ! compressible 787c warning='Solution is set to zero (with p and T to one)' 788c else 789c warning='Solution is set to zero' 790c endif 791c write(*,*) warning 792c endif 793c qold=zero 794c if (matflg(1,1).eq.0) then ! compressible 795c qold(:,1)=one ! avoid zero pressure 796c qold(:,nflow)=one ! avoid zero temperature 797c endif 798c endif 799cc 800c fname1='time derivative of solution?' 801c intfromfile=0 802c call readheader(irstin,fname1,intfromfile, 803c & ithree,'integer', iotype) 804c allocate( acold(nshg,ndof) ) 805c if(intfromfile(1).ne.0) then 806c nshg2=intfromfile(1) 807c ndof2=intfromfile(2) 808c lstep=intfromfile(3) 809c 810c if (nshg2 .ne. nshg) 811c & call error ('restar ', 'nshg ', nshg) 812cc 813c allocate( acread(nshg,ndof2) ) 814c acread=zero 815c 816c iacsiz=nshg*ndof2 817c call readdatablock(irstin,fname1,acread,iacsiz, 818c & 'double',iotype) 819c acold(:,1:ndof)=acread(:,1:ndof) 820c deallocate(acread) 821c else 822c if (myrank.eq.master) then 823c warning='Time derivative of solution is set to zero (SAFE)' 824c write(*,*) warning 825c endif 826c acold=zero 827c endif 828 829c call creadlist(irstin,ithree,nshg2,ndisp,lstep) 830 831 if (ideformwall.eq.1) then 832! fname1='displacement?' 833! call readheader(irstin,fname1,intfromfile, 834! & ithree,'integer', iotype) 835 836 intfromfile=0 837 call phio_readheader(handle, 838 & c_char_'displacement' // char(0), 839 & c_loc(intfromfile), ithree, dataInt, iotype) 840 841 nshg2=intfromfile(1) 842 ndisp=intfromfile(2) 843 lstep=intfromfile(3) 844 if(ndisp.ne.nsd) then 845 warning='WARNING ndisp not equal nsd' 846 write(*,*) warning , ndisp 847 endif 848c 849 if (nshg2 .ne. nshg) 850 & call error ('restar ', 'nshg ', nshg) 851c 852c.... read the values of primitive variables into uold 853c 854 855 allocate( uold(nshg,nsd) ) 856 allocate( uread(nshg,nsd) ) 857 858 iusiz=nshg*nsd 859 860! call readdatablock(irstin,fname1,uread,iusiz, 861! & 'double',iotype) 862 call phio_readdatablock(descriptor,'displacement' // char(0), 863 & uread,iusiz, 'double' // char(0),iotype) 864 865 uold(:,1:nsd)=uread(:,1:nsd) 866 else 867 allocate( uold(nshg,nsd) ) 868 uold(:,1:nsd) = zero 869 endif 870 871c 872c.... close c-binary files 873c 874!MR CHANGE 875! call closefile( irstin, "read" ) 876 877 call phio_closefile_read(handle) 878 879!MR CHANGE 880! call closefile( igeomBAK, "read" ) 881c 882 deallocate(xread) 883 if ( numpbc > 0 ) then 884 deallocate(bcinpread) 885 deallocate(ibctmpread) 886 endif 887 deallocate(iperread) 888 if(numpe.gt.1) 889 & deallocate(ilworkread) 890 deallocate(nbcread) 891 892 return 893c 894 994 call error ('input ','opening ', igeomBAK) 895 995 call error ('input ','opening ', igeomBAK) 896 997 call error ('input ','end file', igeomBAK) 897 998 call error ('input ','end file', igeomBAK) 898c 899 end 900 901c 902c No longer called but kept around in case.... 903c 904 subroutine genpzero(iBC) 905 906 use pointer_data 907c 908 include "common.h" 909 integer iBC(nshg) 910c 911c.... check to see if any of the nodes have a dirichlet pressure 912c 913 pzero=1 914 if (any(btest(iBC,2))) pzero=0 915c 916 do iblk = 1, nelblb 917 npro = lcblkb(1,iblk+1)-lcblkb(1,iblk) 918 do i=1, npro 919 iBCB1=miBCB(iblk)%p(i,1) 920c 921c.... check to see if any of the nodes have a Neumann pressure 922c but not periodic (note that 923c 924 if(btest(iBCB1,1)) pzero=0 925 enddo 926c 927c.... share results with other processors 928c 929 pzl=pzero 930 if (numpe .gt. 1) 931 & call MPI_ALLREDUCE (pzl, pzero, 1, 932 & MPI_DOUBLE_PRECISION,MPI_MIN, MPI_COMM_WORLD,ierr) 933 934 enddo 935c 936c.... return 937c 938 return 939c 940 end 941 942