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 :: dwal(:) 18 real*8, allocatable :: errors(:,:) 19 real*8, allocatable :: ybar(:,:) 20 real*8, allocatable :: yphbar(:,:,:) 21 real*8, allocatable :: vort(:,:) 22 real*8, allocatable :: uold(:,:) 23 real*8, allocatable :: acold(:,:) 24 integer, allocatable :: iBCtmp(:) 25 real*8, allocatable :: BCinp(:,:) 26 27 integer, allocatable :: point2ilwork(:) 28 integer, allocatable :: nBC(:) 29 integer, allocatable :: point2iper(:) 30 integer, allocatable :: point2ifath(:) 31 integer, allocatable :: point2nsons(:) 32 33 end module 34 35 36 37 subroutine readnblk 38c 39 use readarrays 40 include "commonM2NFixBnd.h" 41 include "mpif.h" 42c 43 real*8, allocatable :: xread(:,:), qread(:,:), qread1(:) 44 real*8, allocatable :: uread(:,:), acread(:,:) 45 real*8, allocatable :: BCinpread(:,:) 46 real*8 globmax,globmin 47 integer, allocatable :: iperread(:), iBCtmpread(:) 48 integer, allocatable :: ilworkread(:), nBCread(:) 49 character*10 cname2 50 character*30 fmt1 51 character*255 fname1,fnamer,fnamelr 52 character*255 warning 53 54 integer :: descriptor, color, nfiles, nfields 55 integer :: numparts, nppf 56 character*255 fname2, temp2 57 character*64 temp1 58 integer :: igeom, ibndc, irstin, ierr 59 integer :: ndof, ndoferrors, ndofybar 60 integer :: itmp, itmp2 61 integer :: irstart, irstartmap, iybar 62 integer :: ierror, numphavg, ivort, idwal, idebug, iphavg 63 64 integer intfromfile(50) ! integers read from headers 65 logical exinput 66c 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 74 if(myrank == 0) then 75 fnamer='M2N_input.dat' 76 fnamer = trim(fnamer) // char(0) 77 inquire(file=fnamer,exist=exinput) 78 if(exinput) then 79 open(unit=72,file=fnamer,status='old') 80 read(72,*) irstart 81 read(72,*) irstartmap 82 read(72,*) iybar 83 read(72,*) ierror 84 read(72,*) numphavg 85 read(72,*) ivort 86 read(72,*) idwal 87 read(72,*) idebug 88 close(72) 89 else 90 write(*,*) 'ERROR: Input file ', 91 & trim(fnamer),' does not exist!' 92 endif 93 endif 94 call mpi_barrier(mpi_comm_world, ierr) 95 call mpi_bcast(exinput,1,MPI_LOGICAL,0,mpi_comm_world,ierr) 96 if(.not. exinput) then ! M2NFixBnd_input.dat does not exist. Quit 97! call mpi_abort(mpi_comm_world,911,ierr) 98! call mpi_finalize(ierr) 99 return 100 else ! broadcast the information read by rank 0 101 call mpi_bcast(irstart,1,MPI_INTEGER,0,mpi_comm_world,ierr) 102 call mpi_bcast(iybar,1,MPI_INTEGER,0,mpi_comm_world,ierr) 103 call mpi_bcast(ierror,1,MPI_INTEGER,0,mpi_comm_world,ierr) 104 call mpi_bcast(numphavg,1,MPI_INTEGER,0,mpi_comm_world,ierr) 105 call mpi_bcast(ivort,1,MPI_INTEGER,0,mpi_comm_world,ierr) 106 call mpi_bcast(idwal,1,MPI_INTEGER,0,mpi_comm_world,ierr) 107 endif 108 109 if(myrank == 0 ) then 110 ! Print some info 111 write(*,*) 'The solution field is reduced by default' 112 113 if(iybar .gt. 0) then 114 write(*,*) 'The ybar field will also be reduced' 115 iybar = 1 ! security 116 else 117 write(*,*) 'The ybar field will NOT be reduced' 118 endif 119 120 if(ierror .gt. 0) then 121 write(*,*) 'The error field will also be reduced' 122 ierror = 1 ! security 123 else 124 write(*,*) 'The error field will NOT be reduced' 125 endif 126 127 if(numphavg .gt. 0) then 128 write(*,*) 'The phase average fields (',numphavg, 129 & ') will also be reduced' 130 else 131 write(*,*) 'The phase average fields will NOT be reduced' 132 endif 133 134 if(ivort .gt. 0) then 135 write(*,*) 'The vorticity field will also be reduced' 136 ivort = 1 ! security 137 else 138 write(*,*) 'The vorticity field will NOT be reduced' 139 endif 140 141 if(idwal .gt. 0) then 142 write(*,*) 'The dwal field will also be reduced' 143 idwal = 1 ! security 144 else 145 write(*,*) 'The dwal field will NOT be reduced' 146 endif 147 write(*,*) '' 148 endif 149 150 151 call mpi_barrier(mpi_comm_world, ierr) 152 if(myrank == 0) then 153 write(*,*) 'Reading the geombcRed-dat files' 154 endif 155 156 lstep=irstart ! in case restart files have no fields 157 158 nfiles = nsynciofiles 159 numparts = numpe !This is the common settings. Beware if you try to compute several parts per process 160c 161c.... input the geometry parameters 162c 163 164 color = int(myrank/(numparts/nfiles)) !Should call the color routine in SyncIO here 165 itmp2 = int(log10(float(color+1)))+1 166 write (temp2,"('(''geombcRed-dat.'',i',i1,')')") itmp2 167 write (fnamer,temp2) (color+1) 168 fnamer = trim(fnamer)//char(0) 169 170 ieleven=11 171 ione=1 172 173 itmp = int(log10(float(myrank+1)))+1 174 175 call queryphmpiio(fnamer, nfields, nppf); 176 if (myrank == 0) then 177 write(*,*) 'Number of fields in geombcRed-dat: ',nfields 178 write(*,*) 'Number of parts per file geombcRed-dat: ',nppf 179 endif 180 call initphmpiio( nfields, nppf, nfiles, igeom, 181 & 'read'//char(0)) 182 call openfile( fnamer, 'read'//char(0), igeom ) 183 184 write (temp1,"('(''number of nodes@'',i',i1,',A1)')") itmp 185 write (fname2,temp1) (myrank+1),'?' 186 call readheader(igeom,fname2//char(0),numnp,ione, 187 & 'integer'//char(0), iotype) 188 189 write (temp1,"('(''number of modes@'',i',i1,',A1)')") itmp 190 write (fname2,temp1) (myrank+1),'?' 191 call readheader(igeom,fname2//char(0),nshg,ione, 192 & 'integer'//char(0), iotype) 193 194 write (temp1,"('(''number of interior elements@'',i',i1,',A1)')") 195 & itmp 196 write (fname2,temp1) (myrank+1),'?' 197 call readheader(igeom,fname2//char(0),numel,ione, 198 & 'integer'//char(0), iotype) 199 200 write (temp1,"('(''number of boundary elements@'',i',i1,',A1)')") 201 & itmp 202 write (fname2,temp1) (myrank+1),'?' 203 call readheader(igeom,fname2//char(0),numelb,ione, 204 & 'integer'//char(0),iotype) 205 206 write (temp1, 207 & "('(''maximum number of element nodes@'',i',i1,',A1)')") itmp 208 write (fname2,temp1) (myrank+1),'?' 209 call readheader(igeom,fname2//char(0),nen,ione, 210 &'integer'//char(0),iotype) 211 212 write (temp1,"('(''number of interior tpblocks@'',i',i1,',A1)')") 213 & itmp 214 write (fname2,temp1) (myrank+1),'?' 215 call readheader(igeom,fname2//char(0),nelblk,ione, 216 & 'integer'//char(0) ,iotype) 217 218 write (temp1,"('(''number of boundary tpblocks@'',i',i1,',A1)')") 219 & itmp 220 write (fname2,temp1) (myrank+1),'?' 221 call readheader(igeom,fname2//char(0),nelblb,ione, 222 & 'integer'//char(0), iotype) 223 224 write (temp1, 225 & "('(''number of nodes with Dirichlet BCs@'',i',i1,',A1)')") itmp 226 write (fname2,temp1) (myrank+1),'?' 227 call readheader(igeom,fname2//char(0),numpbc,ione, 228 & 'integer'//char(0),iotype) 229 230 write (temp1,"('(''number of shape functions@'',i',i1,',A1)')") 231 & itmp 232 write (fname2,temp1) (myrank+1),'?' 233 call readheader(igeom,fname2//char(0),ntopsh,ione, 234 & 'integer'//char(0),iotype) 235 236 call closefile( igeom, "read"//char(0) ) 237 call finalizephmpiio( igeom ) 238 239c 240c.... calculate the maximum number of boundary element nodes 241c 242 nenb = 3 !was initialized to 0 but 243 do i = 1, melCat !melCat is 0 here 244 if (nen .eq. nenCat(i,nsd)) nenb = max(nenCat(i,nsd-1), nenb) 245 enddo 246c 247 if (myrank == master) then 248 if (nenb .eq. 0) call error ('input '//char(0), 249 & 'nen '//char(0),nen) 250 endif 251c 252c.... setup some useful constants 253c 254 I3nsd = nsd / 3 ! nsd=3 integer flag 255 E3nsd = float(I3nsd) ! nsd=3 real flag 256c 257 if(matflg(1,1).lt.0) then 258 nflow = nsd + 1 259 else 260 nflow = nsd + 2 261 endif 262 ndof = nsd + 2 263 nsclr=impl(1)/100 264 ndof=ndof+nsclr ! number of sclr transport equations to solve 265 266 ndofBC = ndof + I3nsd ! dimension of BC array 267 ndiBCB = 2 ! dimension of iBCB array 268 ndBCB = ndof + 1 ! dimension of BCB array 269c 270 nsymdf = (ndof*(ndof + 1)) / 2 ! symm. d.o.f.'s 271 272 273 274c 275c.... Read restart files 276c 277c.... read the header and check it against the run data 278c 279 call mpi_barrier(mpi_comm_world, ierr) 280 if(myrank == 0) then 281 write(*,*) 'Reading the RestartRedTmp-dat files' 282 endif 283 284! Beware in what follows! nshg2 read from the header of the solution, 285! dwal, errors and ybar can be different from nshg read from the geombc files 286! This is due to the fact that nshg2 represents the largest vID referenced 287! by the mapping and not the largest vID present in the mesh. 288! qold, dwal, errors, ybar must be allocated to nshg and initialized 289! to a large negative value. The unmapped vertices will be updated 290! through communication with ilwork. 291 292 itmp=1 293 if (irstart .gt. 0) itmp = int(log10(float(irstart+1)))+1 294 295 write (fmt1,"('(''restartRedTmp-dat.'',i',i1,',1x)')") itmp 296 write (fnamer,fmt1) irstart 297 fnamer = trim(fnamer) // cname2(color+1) 298 299 call queryphmpiio(fnamer//char(0), nfields, nppf); 300 if (myrank == 0) then 301 write(*,*) 'Number of fields in restartRedTmp-dat: ',nfields 302 write(*,*) 'Number of parts per file restartRedTmp-dat: ',nppf 303 endif 304 305 call initphmpiio(nfields,nppf,nfiles,descriptor, 306 & 'read'//char(0)) 307 call openfile( fnamer//char(0) , 308 & 'read'//char(0), descriptor ) 309 310c 311c Read the solution 312c 313 ithree=3 314 itmp = int(log10(float(myrank+1)))+1 315 write (temp1,"('(''solution@'',i',i1,',A1)')") itmp 316 write (fname1,temp1) (myrank+1),'?' 317 fname1 = trim(fname1) 318 319 intfromfile=0 320 call readheader(descriptor,fname1//char(0) ,intfromfile, 321 & ithree,'integer'//char(0), iotype) 322 323c 324c.... read the values of primitive variables into q 325c 326 if(intfromfile(1).ne.0) then 327 nshg2=intfromfile(1) 328 ndof2=intfromfile(2) 329 ndof=ndof2 !This must be the same anyway 330 allocate( qold(nshg,ndof) ) 331 qold(:,:) = -9.87654321e32 332 lstep=intfromfile(3) 333 allocate( qread(nshg2,ndof2) ) 334 335 if (nshg2 .ne. nshg) then 336 write(*,*) 'nshg from geombc and nshg2 from restart differ' 337 & //' on rank', myrank, ' :',nshg,nshg2, 338 & ' - Probably mixing phasta files' 339 endif 340 call mpi_barrier(mpi_comm_world, ierr) 341 342 iqsiz=nshg2*ndof2 343 344 call readdatablock(descriptor,fname1//char(0),qread,iqsiz, 345 & 'double'//char(0),iotype) 346 qold(1:nshg2,1:ndof2)=qread(1:nshg2,1:ndof2) 347 deallocate(qread) 348 else 349 if (myrank.eq.master) then 350 if (matflg(1,1).eq.0) then ! compressible 351 warning='Solution is set to zero (with p and T to one)' 352 else 353 warning='Solution is set to zero' 354 endif 355 write(*,*) warning// char(0) 356 endif 357 qold=zero 358 if (matflg(1,1).eq.0) then ! compressible 359 qold(:,1)=one ! avoid zero pressure 360 qold(:,nflow)=one ! avoid zero temperature 361 endif 362 endif 363 364c 365c Read the ybar 366c 367 ndofybar = 0 368 if(iybar==1) then 369 itmp = int(log10(float(myrank+1)))+1 370 write (temp1,"('(''ybar@'',i',i1,',a1)')") itmp 371 write (fname1,temp1) (myrank+1),'?' 372 fname1 = trim(fname1) 373 intfromfile=0 374 call readheader(descriptor,fname1//char(0),intfromfile, 375 & ithree,'integer'//char(0),iotype) 376 nshg2=intfromfile(1) 377 ndofybar=intfromfile(2) 378 !lstep=intfromfile(3) 379 if(ndofybar .ne. 0) then 380 allocate( ybar(nshg,ndofybar) ) 381 ybar(:,:) = -9.87654321e32 382 allocate( qread(nshg2,ndofybar) ) 383 384 iqsiz=nshg2*ndofybar 385 call readdatablock(descriptor,fname1//char(0) ,qread,iqsiz, 386 & 'double'//char(0),iotype) 387 ybar(1:nshg2,1:ndofybar)=qread(1:nshg2,1:ndofybar) 388 deallocate(qread) 389 else 390 write(*,*) 'WARNING: ybar is missing in the restart files' 391 iybar = 0 392 endif 393 endif 394 395c 396c Read the errors 397c 398 ndoferrors=0 399 if(ierror == 1) then 400 write (temp1,"('(''errors@'',i',i1,',a1)')") 401 & itmp 402 write (fname1,temp1) (myrank+1),'?' 403 fname1 = trim(fname1) 404 intfromfile=0 405 call readheader(descriptor,fname1//char(0),intfromfile, 406 & ithree,'integer'//char(0),iotype) 407 nshg2=intfromfile(1) 408 ndoferrors=intfromfile(2) 409 !lstep=intfromfile(3) 410 if(ndoferrors .ne. 0) then 411 allocate( errors(nshg,ndoferrors) ) 412 errors(:,:) = -9.87654321e32 413 allocate( qread(nshg2,ndoferrors) ) 414 iqsiz=nshg2*ndoferrors 415 call readdatablock(descriptor,fname1//char(0),qread,iqsiz, 416 & 'double'//char(0),iotype) 417 errors(1:nshg2,1:ndoferrors)=qread(1:nshg2,1:ndoferrors) 418 deallocate(qread) 419 else 420 if(myrank==0) then 421 write(*,*) 'WARNING: errors is missing in the restart files' 422 endif 423 ierror = 0 424 endif 425 endif 426 427 428! 429! Read the phase_average fields 430! 431 ndofyphbar=0 432 if(numphavg .gt. 0) then 433 do iphavg = 1,numphavg 434 itmp = int(log10(float(myrank+1)))+1 435 itmp2 = int(log10(float(iphavg)))+1 436 write (temp1, 437 & "('(''phase_average'',i',i1,',''@'',i',i1,',A1)')") 438 & itmp2, itmp 439 write (fname1,temp1) iphavg,(myrank+1),'?' 440 fname1 = trim(fname1) 441 intfromfile=0 442 call readheader(descriptor,fname1//char(0),intfromfile, 443 & ithree,'integer'//char(0),iotype) 444 445 nshg2=intfromfile(1) 446 ndofyphbar=intfromfile(2) 447 !lstep=intfromfile(3) 448 449 if(ndofyphbar.ne.0) then 450 ! Allocate some memory for the first ts only 451 if(iphavg==1) then 452 allocate( yphbar(nshg,ndofyphbar,numphavg) ) 453 yphbar(:,:,:) = -9.87654321e32 454 endif 455 456 allocate( qread(nshg2,ndofyphbar) ) 457 iqsiz = nshg2*ndofyphbar 458 call readdatablock(descriptor,fname1//char(0),qread,iqsiz, 459 & 'double'//char(0),iotype) 460 yphbar(1:nshg2,1:ndofyphbar,iphavg) = 461 & qread(1:nshg2,1:ndofyphbar) 462 deallocate(qread) 463 else 464 if(myrank==0) then 465 write(*,*)'WARNING: phase_average is missing '// 466 & 'in the restart files' 467 endif 468 numphavg = 0 469 if(iphavg > 1) then 470 deallocate(yphbar) 471 endif 472 exit 473 endif 474 enddo 475 endif 476 477! 478! follow the usual convention for loading the vorticity field 479! 480 ndofvort=0 481 if(ivort == 1) then 482 itmp = int(log10(float(myrank+1)))+1 483 write (temp1,"('(''vorticity@'',i',i1,',a1)')") itmp 484 write (fname1,temp1) (myrank+1),'?' 485 fname1 = trim(fname1) 486 intfromfile=0 487 call readheader(descriptor,fname1//char(0),intfromfile, 488 & ithree,'integer'//char(0),iotype) 489 490 nshg2=intfromfile(1) 491 ndofvort=intfromfile(2) 492 !lstep=intfromfile(3) 493 494 if(ndofvort .ne. 0) then 495 allocate(vort(nshg,ndofvort)) 496 vort(:,:) = -9.87654321e32 497 allocate(qread(nshg2,ndofvort)) 498 iqsiz = nshg2*ndofvort 499 call readdatablock(descriptor,fname1//char(0),qread,iqsiz, 500 & 'double'//char(0),iotype) 501 vort(1:nshg2,1:ndofvort)=qread(1:nshg2,1:ndofvort) 502 deallocate(qread) 503 else 504 if(myrank==0) then 505 write(*,*) 'WARNING: vorticity is missing '// 506 & 'in the restart files' 507 endif 508 ivort = 0 509 endif 510 endif 511 512c 513c Read the dwal 514c 515 ndofdwal=0 516 if(idwal==1) then 517 write (temp1,"('(''dwal@'',i',i1,',a1)')") itmp 518 write (fname1,temp1) (myrank+1),'?' 519 fname1 = trim(fname1) 520 intfromfile=0 521 call readheader(descriptor,fname1//char(0),intfromfile, 522 & ithree,'integer'//char(0),iotype) 523 524 nshg2=intfromfile(1) 525 ndofdwal=intfromfile(2) 526 if(ndofdwal .ne. 0) then 527 if(ndofdwal.ne.1) then 528 warning='WARNING: ndofdwal not equal 1' 529 write(*,*) warning, ndofdwal 530 endif 531 allocate( dwal(nshg) ) 532 dwal(:) = -9.87654321e32 533 allocate( qread1(nshg2) ) 534 iqsiz=nshg2*1 535 call readdatablock(descriptor,fname1//char(0),qread1,iqsiz, 536 & 'double'//char(0),iotype) 537 dwal(1:nshg2)=qread1(1:nshg2) 538 deallocate(qread1) 539 else 540 if(myrank==0) then 541 write(*,*) 'WARNING: dwal is missing in the restart files' 542 endif 543 idwal = 0 544 endif 545 endif 546 547! 548!.... close c-binary files 549! 550 call closefile( descriptor, "read"//char(0) ) 551 call finalizephmpiio( descriptor ) 552c.... ----------------------> Communication tasks <-------------------- 553c 554 555 if(numpe > 1) then 556 557 call mpi_barrier(mpi_comm_world, ierr) 558 if(myrank == 0) then 559 write(*,*) 'Reading the geombc-dat files again for ilwork' 560 endif 561 562 color = int(myrank/(numparts/nfiles)) !Should call the color routine in SyncIO here 563 itmp2 = int(log10(float(color+1)))+1 564 write (temp2,"('(''geombcRed-dat.'',i',i1,')')") itmp2 565 write (fnamer,temp2) (color+1) 566 fnamer = trim(fnamer)//char(0) 567 568 ieleven=11 569 ione=1 570 571 itmp = int(log10(float(myrank+1)))+1 572 573 call queryphmpiio(fnamer, nfields, nppf); 574 if (myrank == 0) then 575 write(*,*) 'Number of fields in geombcRed-dat: ',nfields 576 write(*,*) 'Number of parts per file geombcRed-dat: ',nppf 577 endif 578 call initphmpiio( nfields, nppf, nfiles, igeom, 579 & 'read'//char(0)) 580 call openfile( fnamer, 'read'//char(0), igeom ) 581 582 write (temp1,"('(''size of ilwork array@'',i',i1,',A1)')") itmp 583 write (fname2,temp1) (myrank+1),'?' 584 call readheader(igeom,fname2//char(0),nlwork,ione, 585 & 'integer' // char(0) ,iotype) 586 587 write (temp1,"('(''ilwork@'',i',i1,',A1)')") itmp 588 write (fname2,temp1) (myrank+1),'?' 589 call readheader(igeom,fname2//char(0) ,nlwork,ione, 590 & 'integer'//char(0) , iotype) 591 592 allocate( point2ilwork(nlwork) ) 593 allocate( ilworkread(nlwork) ) 594 call readdatablock(igeom,fname2//char(0),ilworkread, 595 & nlwork,'integer'//char(0) , iotype) 596 597 point2ilwork = ilworkread 598 deallocate(ilworkread) 599 600 call closefile( igeom, "read"//char(0) ) 601 call finalizephmpiio( igeom ) 602 603 call ctypes (point2ilwork) 604 605 else 606 nlwork=1 607 allocate( point2ilwork(1)) 608 nshg0 = nshg 609 endif 610 611c 612c.... --------------------> communications <------------------------- 613c 614 call mpi_barrier(mpi_comm_world, ierr) ! make sure every rank is synced here 615 if(myrank == 0) then 616 write(*,*) 'Updating the vertices on the part boundaries' 617 endif 618 619 if (numpe > 1) then 620 ! solution 621 call commuMax (qold, point2ilwork, ndof, 'in '//char(0)) 622 call commuMax (qold, point2ilwork, ndof, 'out'//char(0)) 623 call mpi_barrier(mpi_comm_world, ierr) ! make sure everybody is done with ilwork 624 625 ! ybar 626 if(iybar == 1) then 627 call commuMax (ybar, point2ilwork, ndofybar, 'in '//char(0)) 628 call commuMax (ybar, point2ilwork, ndofybar, 'out'//char(0)) 629 call mpi_barrier(mpi_comm_world, ierr) ! make sure everybody is done with ilwork 630 endif 631 632 ! errors 633 if(ierror == 1) then 634 call commuMax (errors, point2ilwork, ndoferrors, 635 & 'in '//char(0)) 636 call commuMax (errors, point2ilwork, ndoferrors, 637 & 'out'//char(0)) 638 call mpi_barrier(mpi_comm_world, ierr) ! make sure everybody is done with ilwork 639 endif 640 641 ! phase_average 642 if(numphavg .gt. 0) then 643 do iphavg = 1,numphavg 644 call commuMax (yphbar(:,:,iphavg), point2ilwork, 645 & ndofyphbar, 'in '//char(0)) 646 call commuMax (yphbar(:,:,iphavg), point2ilwork, 647 & ndofyphbar, 'out'//char(0)) 648 call mpi_barrier(mpi_comm_world, ierr) ! make sure everybody is done with ilwork 649 enddo 650 endif 651 652 ! vorticity 653 if(ivort == 1) then 654 call commuMax (vort, point2ilwork, ndofvort, 'in '//char(0)) 655 call commuMax (vort, point2ilwork, ndofvort, 'out'//char(0)) 656 call mpi_barrier(mpi_comm_world, ierr) ! make sure everybody is done with ilwork 657 endif 658 659 ! dwal 660 if(idwal == 1) then 661 call commuMax (dwal, point2ilwork, 1, 'in '//char(0)) 662 call commuMax (dwal, point2ilwork, 1, 'out'//char(0)) 663 call mpi_barrier(mpi_comm_world, ierr) ! make sure everybody is done with ilwork 664 endif 665 endif 666 667c 668c.... --------------------> Print Min and Max of each field <------------------------- 669c 670 call mpi_barrier(mpi_comm_world, ierr) 671 if (myrank == 0) then 672 write(*,*) 'Printing min and max of each field component' 673 write(*,*) '' 674 endif 675 676 ! qold 677 do idof=1,ndof 678 call mpi_allreduce(maxval(qold(:,idof)),globmax,1,MPI_DOUBLE, 679 & MPI_MAX, mpi_comm_world, ierr ) 680 call mpi_allreduce(minval(qold(:,idof)),globmin,1,MPI_DOUBLE, 681 & MPI_MIN, mpi_comm_world, ierr ) 682 if (myrank == 0) then 683 write(*,925) 'max/min qold(',idof,'): ',globmax,globmin 684 endif 685 enddo 686 if (myrank == 0) then 687 write(*,*) '' 688 endif 689 690 ! ybar 691 if(iybar == 1) then 692 do idof=1,ndofybar 693 call mpi_allreduce(maxval(ybar(:,idof)),globmax,1,MPI_DOUBLE, 694 & MPI_MAX, mpi_comm_world, ierr ) 695 call mpi_allreduce(minval(ybar(:,idof)),globmin,1,MPI_DOUBLE, 696 & MPI_MIN, mpi_comm_world, ierr ) 697 if (myrank == 0) then 698 write(*,925) 'max/min ybar(',idof,'): ',globmax,globmin 699 endif 700 enddo 701 if (myrank == 0) then 702 write(*,*) '' 703 endif 704 endif 705 706 ! errors 707 if(ierror == 1) then 708 do idof=1,ndoferrors 709 call mpi_allreduce(maxval(errors(:,idof)),globmax,1, 710 & MPI_DOUBLE,MPI_MAX, mpi_comm_world, ierr ) 711 call mpi_allreduce(minval(errors(:,idof)),globmin,1, 712 & MPI_DOUBLE,MPI_MIN, mpi_comm_world, ierr ) 713 if (myrank == 0) then 714 write(*,925) 'max/min errors(',idof,'): ',globmax,globmin 715 endif 716 enddo 717 if (myrank == 0) then 718 write(*,*) '' 719 endif 720 endif 721 722 ! phase_average 723 if(numphavg .gt. 0) then 724 do iphavg = 1,numphavg 725 do idof=1,ndofyphbar 726 call mpi_allreduce(maxval(yphbar(:,idof,iphavg)), globmax,1, 727 & MPI_DOUBLE, MPI_MAX, mpi_comm_world, ierr ) 728 call mpi_allreduce(minval(yphbar(:,idof,iphavg)), globmin,1, 729 & MPI_DOUBLE, MPI_MIN, mpi_comm_world, ierr ) 730 if (myrank == 0) then 731 write(*,926) 'max/min yphbar(',idof,iphavg,'):', 732 & globmax,globmin 733 endif 734 enddo 735 if (myrank == 0) then 736 write(*,*) '' 737 endif 738 enddo 739 endif 740 741 ! vorticity 742 if(ivort == 1) then 743 do idof=1,ndofvort 744 call mpi_allreduce(maxval(vort(:,idof)),globmax,1, 745 & MPI_DOUBLE,MPI_MAX, mpi_comm_world, ierr ) 746 call mpi_allreduce(minval(vort(:,idof)),globmin,1, 747 & MPI_DOUBLE,MPI_MIN, mpi_comm_world, ierr ) 748 if (myrank == 0) then 749 write(*,925) 'max/min vort(',idof,'): ',globmax,globmin 750 endif 751 enddo 752 if (myrank == 0) then 753 write(*,*) '' 754 endif 755 endif 756 757 ! dwal 758 if(idwal == 1) then 759 call mpi_allreduce(maxval(dwal(:)), globmax, 1, MPI_DOUBLE, 760 & MPI_MAX, mpi_comm_world, ierr ) 761 call mpi_allreduce(minval(dwal(:)), globmin, 1, MPI_DOUBLE, 762 & MPI_MIN, mpi_comm_world, ierr ) 763 if (myrank == 0) then 764 write(*,925) 'max/min dwal(',1,'):',globmax,globmin 765 write(*,*) '' 766 endif 767 endif 768 769925 format(A,i2,A,2(e24.17,1x)) 770926 format(A,i2,i2,A,2(e24.17,1x)) 771 772c 773c.... e-------------------> Write data to disks <------------------------- 774c 775 call mpi_barrier(mpi_comm_world, ierr) 776 if(myrank == 0) then 777 write(*,*) 'Writing the reduced restartRed-dat files' 778 endif 779 780! call Write_M2NFixBnd(myrank, lstep, nshg, 781! & ndof, ndofybar, ndoferrors, 782! & qold, ybar, errors, dwal) 783 784 nsynciofieldswriterestart = 1+iybar+ierror+numphavg+ivort+idwal 785 call Write_M2NFixBnd_SolOnly(myrank, lstep, nshg, 786 & ndof, qold) 787 deallocate(qold) 788 789 ! ybar 790 if(iybar == 1) then 791 call Write_Field(myrank,'a','ybar',4,ybar,'d', 792 & nshg,ndofybar,lstep) 793 deallocate(ybar) 794 endif 795 796 ! errors 797 if(ierror == 1) then 798 call Write_Field(myrank,'a','errors',6,errors,'d', 799 & nshg,ndoferrors,lstep) 800 deallocate(errors) 801 endif 802 803 ! phase_average 804 if(numphavg .gt. 0) then 805 do iphavg=1,numphavg 806 call write_phavg2(myrank,'a','phase_average',13, 807 & iphavg,numphavg,yphbar(:,:,iphavg),'d', 808 & nshg,ndofyphbar,lstep) 809 enddo 810 deallocate(yphbar) 811 endif 812 813 ! vorticity 814 if(ivort == 1) then 815 call Write_Field(myrank,'a','vorticity',9,vort,'d', 816 & nshg,ndofvort,lstep) 817 deallocate(vort) 818 endif 819 820 ! dwal 821 if(idwal == 1) then 822 call Write_Field(myrank,'a','dwal',4,dwal,'d', 823 & nshg,1,lstep) 824 deallocate(dwal) 825 endif 826 827 828! 829! Deallocate some remaining memory 830! 831 if(numpe.gt.1) then 832 deallocate(point2ilwork) 833 endif 834 835 return 836c 837 994 call error ('input ','opening ', igeom) 838 995 call error ('input ','opening ', igeom) 839 997 call error ('input ','end file', igeom) 840 998 call error ('input ','end file', igeom) 841c 842 end 843