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 :: uold(:,:) 21 real*8, allocatable :: acold(:,:) 22 integer, allocatable :: iBCtmp(:) 23 real*8, allocatable :: BCinp(:,:) 24 25 integer, allocatable :: point2ilwork(:) 26 integer, allocatable :: nBC(:) 27 integer, allocatable :: point2iper(:) 28 integer, allocatable :: point2ifath(:) 29 integer, allocatable :: point2nsons(:) 30 31 end module 32 33 34 35 subroutine readnblk 36c 37 use readarrays 38 include "commonAcuStat.h" 39 include "mpif.h" 40c 41 real*8, allocatable :: xread(:,:),qread(:,:) 42 real*8, allocatable :: qread1(:) 43 real*8, allocatable :: uread(:,:), acread(:,:) 44 real*8, allocatable :: BCinpread(:,:) 45 real*8 globmax,globmin 46 integer, allocatable :: iperread(:), iBCtmpread(:) 47 integer, allocatable :: ilworkread(:), nBCread(:) 48 character*10 cname2 49 character*30 fmt1 50 character*255 fname1,fnamer,fnamelr 51 character*255 warning 52 53 integer :: descriptor, color, nfiles, nfields 54 integer :: numparts, nppf 55 character*255 fname2, temp2 56 character*64 temp1 57 integer :: igeom, ibndc, irstin, ierr 58 integer :: ndof, ndoferrors, ndofybar, ndofyphbar 59 integer :: itmp, itmp2 60 61 integer intfromfile(50) ! integers read from headers 62 logical exinput 63 64 integer :: numts, numphavg, its, isize, sumts, totsumts, debug 65 integer :: ierror, idwal, iphavg, ivort, ndofvort 66 integer :: descriptorG,GPID 67 integer :: timestep, ithree 68 integer :: nshg2,ndof2 69 integer, allocatable :: tablets(:,:) 70 real*8, allocatable :: qybarcumul(:,:),qerrorcumul(:,:) 71 real*8, allocatable :: qyphbarcumul(:,:,:), qvorticity(:,:) 72 73 debug = 0 74 75 numts=-1 76 77 if(myrank == 0) then 78 fnamer='AcuStat_input.dat' 79 fnamer = trim(fnamer) // char(0) 80 inquire(file=fnamer,exist=exinput) 81 if(exinput) then 82 open(unit=24,file=fnamer,status='old') 83 read(24,*) numts 84 read(24,*) ierror 85 read(24,*) numphavg 86 read(24,*) ivort 87 read(24,*) idwal 88 else 89 write(*,*) 'ERROR: Input file ', 90 & trim(fnamer),' does not exist!' 91 endif 92 endif 93 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 ! AcuStat_input.dat does not exist. Quit 97 return 98 endif 99 100 call MPI_Bcast(numts,1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) 101 allocate(tablets(numts,2)) ; tablets=-1 102 103 if(myrank == 0 ) then 104 ! Read the series of time steps 105 do its=1,numts 106 read(24,*) tablets(its,1), tablets(its,2) 107 enddo 108 close(24) 109 110 ! Print some info 111 write(*,*) 'There are ',numts, 'restart files to read for ybar' 112 write(*,*) 'The solution field from the last time step will be', 113 & ' added' 114 115 if(ierror .gt. 0) then 116 write(*,*) 'The error field will also be accumulated' 117 ierror = 1 ! security 118 else 119 write(*,*) 'The error field will NOT be accumulated' 120 endif 121 122 if(numphavg .gt. 0) then 123 write(*,*) 'The phase average fields (',numphavg, 124 & ') will also be accumulated' 125 else 126 write(*,*) 'The phase average fields will NOT be accumulated' 127 endif 128 129 if(ivort .gt. 0) then 130 write(*,*) 'The vorticity field will also be added' 131 ivort = 1 ! security 132 else 133 write(*,*) 'The vorticity field will NOT be added' 134 endif 135 136 if(idwal .gt. 0) then 137 write(*,*) 'The dwal field will also be added' 138 idwal = 1 ! security 139 else 140 write(*,*) 'The dwal field will NOT be added' 141 endif 142 143 endif 144 145 isize=2*numts 146 call MPI_Bcast(tablets,isize,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) 147 call MPI_Bcast(ierror,1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) 148 call MPI_Bcast(numphavg,1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) 149 call MPI_Bcast(idwal,1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) 150 call MPI_Bcast(ivort,1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr) 151 152! write(*,*) 'my rank is: ',myrank,numts 153 154! do irank=0,7 155! if (irank == myrank) then 156! do its=1,numts 157! write(*,*) myrank,tablets(its,1:2) 158! enddo 159! endif 160! call mpi_barrier(MPI_COMM_WORLD,ierr) 161! enddo 162 163 nfiles = nsynciofiles 164 numparts = numpe 165 color = int(myrank/(numparts/nfiles)) 166 sumts = 0 167 totsumts = 0 168 169 ithree=3 170 171! Loop over the restart files to read and accumulate the ybar field 172 do its=1,numts 173 174! 175! Read the table 176! 177 timestep = tablets(its,1) 178 sumts = tablets(its,2) 179 totsumts = totsumts + sumts 180 181! 182! Get the right restart files 183! 184 itmp=1 185 if (timestep .gt. 0) itmp = int(log10(float(timestep+1)))+1 186 write (fmt1,"('(''restart-dat.'',i',i1,',1x)')") itmp 187 write (fnamer,fmt1) timestep 188 fnamer = trim(fnamer) // cname2(color+1) 189 190 if(debug==1) then 191 write(*,*)myrank,timestep,nfiles,numparts,trim(adjustl(fnamer)) 192 endif 193 194! 195! query, init and open 196! 197 call queryphmpiio(fnamer//char(0), nfields, nppf); 198 if (myrank == 0) then 199 write(*,*) '' 200 write(*,*) 'Considering time step',timestep,sumts 201 write(*,*) 'Number of fields in restart-dat: ',nfields 202 write(*,*) 'Number of parts per file restart-dat: ',nppf 203 endif 204 205 call initphmpiio(nfields,nppf,nfiles,descriptor,'read'//char(0)) 206 call openfile( fnamer//char(0), 'read'//char(0), descriptor ) 207 208! 209! Read the ybar field 210! 211 itmp = int(log10(float(myrank+1)))+1 212 write (temp1,"('(''ybar@'',i',i1,',A1)')") itmp 213 write (fname1,temp1) (myrank+1),'?' 214 fname1 = trim(fname1) 215 216 intfromfile=0 217 call readheader(descriptor,fname1//char(0),intfromfile,ithree, 218 & 'integer'//char(0),iotype) 219 if(intfromfile(1).ne.0) then 220 nshg2=intfromfile(1) !Do not use nshg and ndof from common.h here! 221 ndofybar=intfromfile(2) 222 !lstep2=intfromfile(3) 223 else 224 write(*,*)'ERROR: Did not find the field requested 225 & in AcuStat_input',myrank 226 endif 227 228 ! Allocate some memory for the first ts only 229 if(its==1) then 230 allocate( qybarcumul(nshg2,ndofybar) ) ; qybarcumul = 0.d0 231 endif 232 233 allocate( qread(nshg2,ndofybar) ) ; qread = 0.d0 234 iqsiz = nshg2*ndofybar 235 call readdatablock(descriptor,fname1//char(0),qread,iqsiz, 236 & 'double'//char(0),iotype) 237 238 ! Weight correctly the field and accumulate for average 239 qybarcumul(:,:) = qybarcumul(:,:) + qread (:,:)*sumts 240 deallocate(qread) 241 242! 243! Read the errors 244! 245 if(ierror == 1) then 246 write (temp1,"('(''errors@'',i',i1,',a1)')") 247 & itmp 248 write (fname1,temp1) (myrank+1),'?' 249 fname1 = trim(fname1) 250 intfromfile=0 251 call readheader(descriptor,fname1//char(0),intfromfile,ithree, 252 & 'integer'//char(0),iotype) 253 if(intfromfile(1).ne.0) then 254 nshg2=intfromfile(1) 255 ndoferrors=intfromfile(2) 256 !lstep=intfromfile(3) 257 else 258 write(*,*)'ERROR: Did not find the field requested 259 & in AcuStat_input',myrank 260 endif 261 262 ! Allocate some memory for the first ts only 263 if(its==1) then 264 allocate( qerrorcumul(nshg2,ndoferrors) ) ; qerrorcumul = 0.d0 265 endif 266 267 allocate( qread(nshg2,ndoferrors) ) ; qread = 0.d0 268 iqsiz=nshg2*ndoferrors 269 call readdatablock(descriptor,fname1//char(0),qread,iqsiz, 270 & 'double'//char(0),iotype) 271 272 ! Weight correctly the field and accumulate for average 273 qerrorcumul(:,:) = qerrorcumul(:,:) + qread(:,:)*sumts 274 deallocate(qread) 275 endif ! if ierror == 1 276 277! 278! Read the phase_average fields 279! 280 if(numphavg .gt. 0) then 281 do iphavg = 1,numphavg 282 itmp = int(log10(float(myrank+1)))+1 283 itmp2 = int(log10(float(iphavg)))+1 284! write (temp1,"('(''phase_average@'',i',i1,',A1)')") itmp 285! write(*,*) temp1 ! ('phase_average@',i1,A1) 286 write (temp1, 287 & "('(''phase_average'',i',i1,',''@'',i',i1,',A1)')") 288 & itmp2, itmp 289! write(*,*) temp1 ! ('phase_average',i1,'@',i1,A1) 290 write (fname1,temp1) iphavg,(myrank+1),'?' 291 fname1 = trim(fname1) 292 if(debug == 1) then 293 if(myrank == 0 .or. myrank == numpe-1) then 294 write(*,*)'Rank: ',myrank,'- phase: ', 295 & iphavg,' - field: ',trim(fname1) 296 endif 297 endif 298 intfromfile=0 299 call readheader(descriptor,fname1//char(0),intfromfile, 300 & ithree,'integer'//char(0),iotype) 301 if(intfromfile(1).ne.0) then 302 nshg2=intfromfile(1) !Do not use nshg and ndof from common.h here! 303 ndofyphbar=intfromfile(2) 304 !lstep2=intfromfile(3) 305 else 306 write(*,*)'ERROR: Did not find the field requested 307 & in AcuStat_input',myrank 308 endif 309 310 ! Allocate some memory for the first ts only 311 if(its==1 .and. iphavg==1) then 312 allocate( qyphbarcumul(nshg2,ndofyphbar,numphavg) ) 313 qyphbarcumul = 0.d0 314 endif 315 316 allocate( qread(nshg2,ndofyphbar) ) ; qread = 0.d0 317 iqsiz = nshg2*ndofyphbar 318 call readdatablock(descriptor,fname1//char(0),qread,iqsiz, 319 & 'double'//char(0),iotype) 320 321 ! Weight correctly the field and accumulate for average 322 qyphbarcumul(:,:,iphavg) = qyphbarcumul(:,:,iphavg) 323 & + qread (:,:)*sumts 324 deallocate(qread) 325 enddo 326 endif 327! 328! If last ts, read also the solution, error and dwal field 329! for reduction with M2N or for convenience 330! 331 332 if(its == numts) then 333! 334! Read the solution of the last time step 335! 336 itmp = int(log10(float(myrank+1)))+1 337 write (temp1,"('(''solution@'',i',i1,',A1)')") itmp 338 write (fname1,temp1) (myrank+1),'?' 339 fname1 = trim(fname1) 340 341 intfromfile=0 342 call readheader(descriptor,fname1//char(0) ,intfromfile, 343 & ithree,'integer'//char(0), iotype) 344 345 if(intfromfile(1).ne.0) then 346 nshg2=intfromfile(1) 347 ndof2=intfromfile(2) 348 allocate( qold(nshg2,ndof2) ) 349 qold(:,:) = -9.87654321e32 350 lstep=intfromfile(3) 351 allocate( qread(nshg2,ndof2) ) ; qread = 0.d0 352 iqsiz=nshg2*ndof2 353 call readdatablock(descriptor,fname1//char(0),qread,iqsiz, 354 & 'double'//char(0),iotype) 355 qold(1:nshg2,1:ndof2) = qread(1:nshg2,1:ndof2) 356 deallocate(qread) 357 else 358 write(*,*)'ERROR: Did not find the field requested 359 & in AcuStat_input',myrank 360 endif 361 362! 363! Read the vorticity field from the last time step 364! 365 if(ivort == 1) then 366 write (temp1,"('(''vorticity@'',i',i1,',a1)')") 367 & itmp 368 write (fname1,temp1) (myrank+1),'?' 369 fname1 = trim(fname1) 370 intfromfile=0 371 call readheader(descriptor,fname1//char(0),intfromfile, 372 & ithree,'integer'//char(0),iotype) 373 if(intfromfile(1).ne.0) then 374 nshg2=intfromfile(1) 375 ndofvort=intfromfile(2) 376 allocate( qvorticity(nshg2,ndofvort) ) 377 qvorticity(:,:) = -9.87654321e32 378 allocate( qread(nshg2,ndofvort) ) ; qread = 0.d0 379 iqsiz=nshg2*ndofvort 380 call readdatablock(descriptor,fname1//char(0),qread, 381 & iqsiz,'double'//char(0),iotype) 382 qvorticity(1:nshg2,1:ndofvort) = qread(1:nshg2,1:ndofvort) 383 deallocate(qread) 384 else 385 write(*,*)'ERROR: Did not find the field requested 386 & in AcuStat_input',myrank 387 endif 388 endif 389 390! 391! Read the dwal 392! 393 if(idwal == 1) then 394 write (temp1,"('(''dwal@'',i',i1,',a1)')") 395 & itmp 396 write (fname1,temp1) (myrank+1),'?' 397 fname1 = trim(fname1) 398 intfromfile=0 399 call readheader(descriptor,fname1//char(0),intfromfile, 400 & ithree,'integer'//char(0),iotype) 401 if(intfromfile(1).ne.0) then 402 nshg2=intfromfile(1) 403 allocate( dwal(nshg2) ) 404 dwal(:) = -9.87654321e32 405 allocate( qread1(nshg2) ) ; qread1 = 0.d0 406 iqsiz=nshg2*1 407 call readdatablock(descriptor,fname1//char(0),qread1, 408 & iqsiz,'double'//char(0),iotype) 409 dwal(1:nshg2) = qread1(1:nshg2) 410 deallocate(qread1) 411 else 412 write(*,*)'ERROR: Did not find the field requested 413 & in AcuStat_input',myrank 414 endif 415 endif !if idwal == 1 416 417 418 419 endif ! First time step 420! 421! Close and finalize 422! 423 call closefile( descriptor, 'read'//char(0) ) 424 call finalizephmpiio( descriptor ) 425 426 enddo !numts = number of restart files to read 427 428! 429! Weight correctly the average field 430! 431 qybarcumul(:,:) = qybarcumul(:,:)/totsumts 432 if (ierror == 1) then 433 qerrorcumul(:,:) = qerrorcumul(:,:)/totsumts 434 endif 435 if(numphavg .gt. 0) then 436 do iphavg = 1,numphavg 437 qyphbarcumul(:,:,iphavg) = qyphbarcumul(:,:,iphavg)/totsumts 438 enddo 439 endif 440 441! Write qold, dwal, errors and qybarcumul (= cumulated ybar) 442 443 nsynciofieldswriterestart = 2 + ierror + numphavg + ivort + idwal ! 2 for solution and accumulated ybar 444! write(*,*) 'nsynciofieldswriterestart:',nsynciofieldswriterestart 445 446! call Write_AcuStat(myrank, tablets(1,1), tablets(numts,1), 447! & totsumts, nshg2, ndof2, ndoferrors, ndofybar, 448! & qold, dwal, qerrorcumul, qybarcumul) 449 450 call mpi_barrier(mpi_comm_world, ierr) 451 call Write_AcuStat(myrank, tablets(1,1), tablets(numts,1), 452 & totsumts, nshg2, ndof2, qold) 453 454 if(ivort == 1) then 455 call write_field(myrank,'a','vorticity',9,qvorticity, 456 & 'd', nshg2, ndofvort, tablets(numts,1)) 457 endif 458 459 if (ierror == 1) then 460 call write_error(myrank, tablets(numts,1), nshg2, ndoferrors, 461 & qerrorcumul) 462 endif 463 464 call write_field(myrank,'a', 'ybar', 4, 465 & qybarcumul, 'd', nshg2, ndofybar, tablets(numts,1)) 466 467 if(numphavg .gt. 0) then 468 do iphavg = 1,numphavg 469 call write_phavg2(myrank,'a','phase_average',13,iphavg, 470 & numphavg,qyphbarcumul(:,:,iphavg),'d',nshg2,ndofyphbar, 471 & tablets(numts,1)) 472 enddo 473 endif 474 475 if (idwal == 1) then 476 call write_field(myrank, 'a', 'dwal', 4, dwal, 477 & 'd', nshg2, 1, tablets(numts,1)) 478 endif 479 480! 481! Free memory 482! 483 deallocate(tablets) 484 deallocate(qybarcumul, qold) 485 if(ierror == 1) then 486 deallocate(qerrorcumul) 487 endif 488 if(numphavg .gt. 0) then 489 deallocate(qyphbarcumul) 490 endif 491 if(ivort == 1) then 492 deallocate(qvorticity) 493 endif 494 if(idwal == 1) then 495 deallocate(dwal) 496 endif 497 498 return 499 500 end 501