! OMGMRGP V1.5 09 JAN 2019 ! ! Reconstruct omega file from multiple files (from parallel codes) ! ! D.M. Mitnik & N.R. Badnell ! !........ Reads partitioned omega files in row-wise or column-wise formats !.........and reconstructs a single omega file. (Eventually !........ this option should be implemented on the parallel codes). !......... If the input is column-wise then so is the output. !......... If the input is row-wise then the output can be in rowwise !.........(icolwise=0, default) or columnwise (icolwise=1) format. !........If infinite energy point present set ibige=1 !........If infinite energy point not present set ibige=0 ! Case MXE .gt.0 (e.g. R-matrix) default is ! ibige=1 (if OMEGANOTOP) ! ibige=0 (if OMEGA) ! Case MXE .lt.0 (e.g. AS-DW) default is ibige=1 (i.e. in all) ! implicit real*8 (a-h,o-z) logical bform,bformo character elas*3 NAMELIST/filedata/nproc,ibige,icolwise,elas !ielas allocatable :: einl(:),omegar(:,:),isat(:),lat(:),nwrt(:) allocatable :: omegaw(:),npos(:,:),enerp(:),norder(:) character*1 filec,filed,fileu call cpu_time(timei) iout = 8 !! Output unit number !.......default values ibige = 999 !! Default depends on RM or AS-DW input, see mxe below icolwise = 0 !! Output default is rowwise if Input is rowwise. cold ielas = 0 !! default neutral/elasic ion/non-elastic. -1 overrides. elas=' ' !! defaults to YES/NO for neutral/ion elastic present nproc = 99999 !! Search all possible files open(unit=5,file='domgmrgp',status='unknown') read(5,filedata,end=900) if(icolwise.ne.0)icolwise=1 itop=0 900 continue mxe=0 do 1000 iam=0,nproc-1 ich0 = ichar('0') !suffix for input files omega000, omega001, etc. ic = iam/100 id = (iam - 100*ic)/10 iu = (iam - 100*ic - 10*id) ic = ic + ich0 id = id + ich0 iu = iu + ich0 filec = char(ic) filed = char(id) fileu = char(iu) !.......input files generated from (omega###) -- output: omega inquire (file='omega'//filec//filed//fileu,exist=bform) if(bform)then open(10+iam,file='omega'//filec//filed//fileu, & & status='old') if (iam.eq.0) then open(iout,file='omega',form='formatted',status='unknown') write(*,*)'*** USING FORMATTED omega FILE' endif else inquire (file='omegau'//filec//filed//fileu,exist=bform) if(bform)then open(10+iam,file='omegau'//filec//filed//fileu, & & status='old',form='unformatted') bform=.false. if (iam.eq.0) then open(iout,file='omegau',form='unformatted', & & status='unknown') write(*,*)'*** USING UNFORMATTED omegau FILE' endif else !.......input files generated from (OMEGA###) -- output: OMEGA inquire (file='OMEGA'//filec//filed//fileu,exist=bform) if(bform)then open(10+iam,file='OMEGA'//filec//filed//fileu, & & status='old') if (iam.eq.0) then open(iout,file='OMEGA',form='formatted', & & status='unknown') write(*,*)'*** USING FORMATTED OMEGA FILE' endif else inquire (file='OMEGAU'//filec//filed//fileu,exist=bform) if(bform)then open(10+iam,file='OMEGAU'//filec//filed//fileu, & & status='old',form='unformatted') bform=.false. if (iam.eq.0) then open(iout,file='OMEGAU',form='unformatted', & & status='unknown') write(*,*)'*** USING UNFORMATTED OMEGAU FILE' endif else !.....input files generated from (OMEGANOTOP###) -- output: OMEGANOTOP itop=1 inquire (file='OMEGANOTOP'//filec//filed//fileu, & & exist=bform) if(bform)then open(10+iam,file='OMEGANOTOP'//filec//filed//fileu, & & status='old') if (iam.eq.0) then open(iout,file='OMEGANOTOP',form='formatted', & & status='unknown') write(*,*)'*** USING FORMATTED OMEGANOTOP FILE' endif else inquire (file='OMEGANOTOPU'//filec//filed//fileu, & & exist=bform) if(bform)then open(10+iam,file='OMEGANOTOPU'//filec//filed//fileu, & & status='old',form='unformatted') bform=.false. if (iam.eq.0) then open(iout,file='OMEGANOTOPU',form='unformatted', & & status='unknown') write(*,*)'*** USING UNFORMATTED OMEGANOTOPU FILE' endif else if(iam.eq.0)then stop'Input collision strengths files not found...' else go to 1001 endif endif endif endif endif endif endif !.......start reading omega### files bformo=bform !hold because search past last file sets it .false. if(bform)then read(10+iam,*) nzed,nelc read(10+iam,*) nlev,mxe0,nomwrt !assumes NOMWRT always present now else read(10+iam) nzed,nelc read(10+iam) nlev,mxe0,nomwrt endif if(mxe0.gt.0)then isigne=1 if(ibige.eq.999)ibige=-1 !default R-matrix else isigne=-1 if(ibige.eq.999)ibige=1 !default AS-DW mxe0=iabs(mxe0) endif if(ibige.gt.0)ibige = 1 if(ibige.lt.0)ibige=itop !set default ibige mxe=mxe+mxe0-ibige if(iam.eq.0)then !.......allocate the arrays allocate(isat(nlev),stat=ierr) allocate(lat(nlev),stat=ierr) allocate(einl(nlev),stat=ierr) endif !........read statistical weights and energies if(bform)then read(10+iam,*)(isat(i),lat(i),i=1,nlev) read(10+iam,*) (einl(i),i=1,nlev) else read(10+iam)(isat(i),lat(i),i=1,nlev) read(10+iam) (einl(i),i=1,nlev) endif 1000 continue iam=nproc 1001 nproc=iam nprocm=nproc-1 bform=bformo mxe=mxe+ibige ntr=abs(nomwrt) !.......set-up ione for elastic =0, or not =1, default. if(elas.eq.' ')then if(nzed.eq.nelc)elas='YES' if(nzed.ne.nelc)elas='NO' endif if(elas.eq.'YES')then ione=0 elseif(elas.eq.'NO')then ione=1 else write(6,*)' unrecognized option elas=',elas stop ' unrecognized elas option' endif !.......write omega header to iout isign=icolwise if(nomwrt.lt.0)isign=0 if(bform)then write(iout,*) nzed,nelc write(iout,*)nlev,mxe*isigne,nomwrt*(1-2*isign) write(iout,*)(isat(i),lat(i),i=1,nlev) write(iout,270) (einl(i),i=1,nlev) 270 format(1p5e16.6) else write(iout) nzed,nelc write(iout) nlev,mxe*isigne,nomwrt*(1-2*isign) write(iout)(isat(i),lat(i),i=1,nlev) write(iout) (einl(i),i=1,nlev) endif !.......deallocate isat,lat deallocate(isat) deallocate(lat) !........allocate omegas, energies etc allocate(omegar(ntr,0:nprocm),stat=ierr) allocate(omegaw(ntr),stat=ierr) allocate(npos(nlev,nlev),stat=ierr) allocate(enerp(0:nprocm),stat=ierr) allocate(norder(0:nprocm),stat=ierr) allocate(nwrt(0:nprocm),stat=ierr) !........initialize omegas , nrd do j=0,nprocm do i=1,ntr omegar(i,j)=0 enddo enddo do i=1,ntr omegaw(i)=0 enddo nrd=ntr !.......calculation of the position in columnwise format n=0 do i=1,nlev do j=i+1,nlev n=n+1 npos(i,j)=n enddo enddo iam0=0 iamm=nprocm nproct=nproc !s times=0.0 !.....loop over files. !.....first pass reads each file and writes data for lowest energy !.....subsequent passes only update read of file that was written 100 do iam=iam0,iamm enerp(iam)=1.d6 if(bform)then read(10+iam,*,end=200)ener else read(10+iam,end=200) ener endif enerp(iam)=ener 200 backspace(10+iam) if(enerp(iam).gt.1.d5)then nproct=nproct-1 !reached infinite energy pt or EOF if(nproct.eq.0)go to 300 !no more finite energy omegas go to 400 !to write omega already held endif !.......reading energies and omegas !...........Case 1) Input is columnwise if (nomwrt.lt.0) then !! if input is column, also output !............first determine the number of channels to read if (bform) then read(10+iam,*) ener else read(10+iam) ener endif backspace(10+iam) if (ener.lt.einl(nlev)) then nrd = nopen(ener,einl,nlev,ione) nrd = min(nrd,ntr) else nrd = ntr endif !............now read omega if (bform) then read(10+iam,380) ener,(omegar(n,iam),n=1,nrd) else read(10+iam,err=250) ener,(omegar(n,iam),n=1,nrd) endif 250 nwrt(iam)=nrd else !...........Case 2) Input is rowwise !............. read omega nwrt(iam)=nrd if (bform) then read(10+iam,380) ener,(omegar(n,iam),n=1,nrd) else read(10+iam) ener,(omegar(n,iam),n=1,nrd) endif endif 400 continue enddo !end loop over processor files !.........determine lowest energy !s call cpu_time(timeis) call order(enerp,norder,nprocm) !SORT !s call cpu_time(timefs) !s times=times+timefs-timeis iam=norder(0) !.........write appropriate omega if(nomwrt.lt.0.or.nomwrt.gt.0.and.icolwise.eq.0)then if (bform) then write(iout,380) enerp(iam),(omegar(n,iam),n=1,nwrt(iam)) else write(iout) enerp(iam),(omegar(n,iam),n=1,nwrt(iam)) endif else !...........determine the number of channels open if (enerp(iam).lt.einl(nlev)) then nwrt(iam)=nopen(enerp(iam),einl,nlev,ione) nwrt(iam)=min(nwrt(iam),ntr) else nwrt(iam)=ntr endif n=0 do j=1+ione,nlev do i=1,j-ione n=n+1 if (n.gt.nwrt(iam)) go to 40 omegaw(n)=omegar(npos(i,j),iam) enddo enddo 40 if(bform) then write(iout,380) enerp(iam),(omegaw(n),n=1,nwrt(iam)) else write(iout) enerp(iam),(omegaw(n),n=1,nwrt(iam)) endif endif iam0=iam iamm=iam go to 100 !......end of all finite energy omegas on all files !.......infinite omega 300 if (ibige.ne.0) then iam = 0 if(bform) then read(10+iam,380) ener,(omegar(n,iam),n=1,ntr) write(iout,380) ener,(omegar(n,iam),n=1,ntr) else read(10+iam) ener,(omegar(n,iam),n=1,ntr) write(iout) ener,(omegar(n,iam),n=1,ntr) endif endif 380 format(1pe14.8,6(1pe11.3)/(14x,6(e11.3))) deallocate(einl) deallocate(omegar) deallocate(omegaw) deallocate(npos) deallocate(enerp) deallocate(norder) deallocate(nwrt) !s write(*,*)'Sort time=',times/60.,' min' call cpu_time(timef) time=(timef-timei)/60. write(*,777)time 777 format(' CPU TIME=',f9.3,' MIN') STOP 'omgmrgp: normal end' END !---------------------------------------------------- integer function nopen(e,enat,nast,ione) implicit real*8 (a-h,o-z) dimension enat(nast) i0=1 do i=i0,nast if (e.lt.enat(i)) go to 1 enddo i=nast+1 1 i0=i-1 nopen=(i0*(i0-2*ione+1))/2 return end !---------------------------------------------------- SUBROUTINE ORDER(EN,NORDER,NDIM) IMPLICIT REAL*8 (A-H,O-Z) ! ! ! !----------------------------------------------------------------------- ! ! RETURNS NORDER(I)=POINTER TO I-TH ENERGY IN EN ARRAY, ! !----------------------------------------------------------------------- DIMENSION EN(0:NDIM),NORDER(0:NDIM) !----------------------------------------------------------------------- DO 40 K = 0,NDIM J = K J1 = J - 1 IF (J1.EQ.-1) GOTO 30 E = EN(J) + 1.0D-10 DO 20 I = 0,J1 IF (J.LT.K) GOTO 10 N = NORDER(I) IF (E.GT.EN(N)) GOTO 20 10 CONTINUE NORDER(J) = NORDER(J-1) J = J - 1 20 CONTINUE 30 NORDER(J) = K 40 CONTINUE ! END