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