text:
subroutine outbods(fname)
 include 'globals.h'
 integer,parameter :: ndigits=6
 character(len=ndigits) :: nstring
 character(len=24) :: filenaam
 character(len=*) :: fname
 if(fname(1:6).EQ.'XXXXXX') then
  call itos(nsnap,ndigits,nstring)
  filenaam=TRIM(outputfile)//'.'//nstring
  nsnap=nsnap+1
 else
  filenaam=fname
 endif
  call writebods(filenaam)
end subroutine outbods

subroutine writebods(filenaam)
 include 'globals.h'
 character(len=*) :: filenaam
 integer :: nbods,nkeep,p
 integer(kind=4), parameter :: nihead=32,nrhead=32,nphead=64
 integer(kind=4) :: ioversion
 integer(kind=4) :: header(8),ihead(nihead)
 real(kind=8) :: rhead(nrhead)
 integer(kind=4) :: phead(nphead)

 call prepareoutput

 ioversion=2

 header(1)=nihead
 header(2)=nrhead
 header(3)=nphead
 header(4)=0
 header(5)=0
 header(6)=0
 header(7)=0
 header(8)=0

 ihead(1:nihead)=0
 rhead(1:nrhead)=0
 phead(1:nphead)=0

 ihead(1)=nbodies
 ihead(2)=nsph
 ihead(3)=nstar
 ihead(4)=nbodies-nstar-nsph
 ihead(5)=totptag
 ihead(6)=nbh

 ihead(8)=syncflag
 if(syncflag.NE.0) then
  print*, ' ** warning ** async output??'
  stop
 endif

 rhead(1)=tnow
 rhead(2)=mtot
 rhead(3)=ektot
 rhead(4)=eptot
 rhead(5)=snheat
 rhead(6)=eradiate
 rhead(7)=esofttot
 rhead(8)=enerror
 rhead(9)=estar
 rhead(10)=efuvheat
 rhead(11)=eradcool
 rhead(12)=tpos
 rhead(13)=teth
 rhead(14)=trad
 rhead(15)=meanmwt
 rhead(16)=timescale
 rhead(17)=unitm_in_msun
 rhead(18)=lengthscale
 rhead(19)=densconst
 rhead(20)=massres
 rhead(21)=fhydrogn
 rhead(22)=tstarform
 rhead(23)=tsnfeedback
 rhead(24)=tbh
 rhead(25)=eps

! rhead(31)=redshift
! rhead(32)=redshift0

 phead(1:nphead)=output(1:nphead)
 if(nsph.eq.0) then
  phead(10:25)=0
 endif
 if(nstar.eq.0) then
  phead(33:35)=0
 endif

 open(unit=ubodsout,FILE=filenaam,status='UNKNOWN',form='UNFORMATTED')

 WRITE(ubodsout) ioversion
 WRITE(ubodsout) (header(p),p=1,8)
 WRITE(ubodsout) (ihead(p),p=1,nihead)
 WRITE(ubodsout) (rhead(p),p=1,nrhead)
 WRITE(ubodsout) (phead(p),p=1,nphead)

 IF(phead(1).EQ.1) WRITE(ubodsout) (mass(p),p=1,nbodies)
 IF(phead(2).EQ.1) WRITE(ubodsout) (pos(p,1),p=1,nbodies)
 IF(phead(2).EQ.1) WRITE(ubodsout) (pos(p,2),p=1,nbodies)
 IF(phead(2).EQ.1) WRITE(ubodsout) (pos(p,3),p=1,nbodies)
 IF(phead(3).EQ.1) WRITE(ubodsout) (vel(p,1),p=1,nbodies)
 IF(phead(3).EQ.1) WRITE(ubodsout) (vel(p,2),p=1,nbodies)
 IF(phead(3).EQ.1) WRITE(ubodsout) (vel(p,3),p=1,nbodies)
 IF(phead(4).EQ.1) WRITE(ubodsout) (epsgrav(p),p=1,nbodies)
 IF(phead(5).EQ.1) WRITE(ubodsout) (tform(p),p=1,nbodies)
 IF(phead(6).EQ.1) WRITE(ubodsout) (acc(p,1),p=1,nbodies)
 IF(phead(6).EQ.1) WRITE(ubodsout) (acc(p,2),p=1,nbodies)
 IF(phead(6).EQ.1) WRITE(ubodsout) (acc(p,3),p=1,nbodies)

IF(nsph.GT.0) THEN
 IF(phead(10).EQ.1) WRITE(ubodsout) (rho(p),p=1,nsph)
 IF(phead(11).EQ.1) THEN
  IF(uentropy) THEN
   WRITE(ubodsout) (tempvect(p),p=1,nsph)
  ELSE
   WRITE(ubodsout) (ethermal(p),p=1,nsph)
  ENDIF
 ENDIF
 IF(phead(12).EQ.1) THEN
  IF(uentropy) THEN
   WRITE(ubodsout) (entropy(p),p=1,nsph)
  ELSE
   WRITE(ubodsout) (tempvect(p),p=1,nsph)
  ENDIF
 ENDIF
 IF(phead(13).EQ.1) WRITE(ubodsout) (hsmooth(p),p=1,nsph)
 IF(phead(14).EQ.1) WRITE(ubodsout) (fuvheat(p),p=1,nsph)
 IF(phead(15).EQ.1) WRITE(ubodsout) (esnthdt(p),p=1,nsph)
 IF(phead(16).EQ.1) WRITE(ubodsout) (tcollaps(p),p=1,nsph)
 IF(phead(17).EQ.1) WRITE(ubodsout) (temperat(p),p=1,nsph)
 IF(phead(18).EQ.1) WRITE(ubodsout) (elecfrac(p),p=1,nsph)
 IF(phead(19).EQ.1) WRITE(ubodsout) (csound(p),p=1,nsph)
 IF(phead(20).EQ.1) WRITE(ubodsout) &
  (csound(p)**2*rho(p)/gamma,p=1,nsph)
 IF(phead(21).EQ.1) WRITE(ubodsout) (hsmdivv(p),p=1,nsph)
 IF(phead(22).EQ.1) WRITE(ubodsout) (mumaxdvh(p),p=1,nsph)
 IF(phead(23).EQ.1) WRITE(ubodsout) (hsmcurlv(p),p=1,nsph)
 IF(phead(24).EQ.1) WRITE(ubodsout) (vdisp(p),p=1,nsph)
 IF(phead(25).EQ.1) WRITE(ubodsout) (h2frac(p),p=1,nsph)
 IF(phead(26).EQ.1) THEN
  IF(uentropy) THEN
   WRITE(ubodsout) (tempvect(p),p=nsph+1,2*nsph)
  ELSE
   WRITE(ubodsout) (dethdt(p),p=1,nsph)
  ENDIF
 ENDIF
 IF(phead(27).EQ.1) THEN
  IF(uentropy) THEN
   WRITE(ubodsout) (dentdt(p),p=1,nsph)
  ELSE
   WRITE(ubodsout) (tempvect(p),p=nsph+1,2*nsph)
  ENDIF
 ENDIF
ENDIF

IF(nstar.GT.0) THEN
 IF(phead(34).EQ.1) WRITE(ubodsout) (starfuv(p),p=nbodies-nstar+1,nbodies)
 IF(phead(35).EQ.1) WRITE(ubodsout) (snentropy(p),p=nbodies-nstar+1,nbodies)
ENDIF

 IF(phead(40).EQ.1) WRITE(ubodsout) (phi(p),p=1,nbodies)
 IF(phead(41).EQ.1) WRITE(ubodsout) (phiext(p),p=1,nbodies)
 IF(phead(42).EQ.1) WRITE(ubodsout) (nbexist(p),p=1,nbodies)
 IF(phead(43).EQ.1) WRITE(ubodsout) (itimestp(p),p=1,nbodies)

 close(ubodsout)

end subroutine writebods

subroutine readbods(filenaam)
 include 'globals.h'
 character(len=*) :: filenaam
 integer :: nbods,nkeep,p,ioerr
 integer(kind=4),parameter ::  nihead=32,nrhead=32,nphead=64
 integer(kind=4) :: ioversion
 integer(kind=4) :: header(8),ihead(nihead)
 real(kind=8) :: rhead(nrhead)
 integer(kind=4) :: phead(nphead)

 print*,' ...reading data... ',filenaam
 open(unit=ubodsin,FILE=filenaam,status='OLD',form='UNFORMATTED',IOSTAT=ioerr)

 IF(ioerr.NE.0) THEN
  print*, ' stop -- error reading file ', filenaam
  stop
 ENDIF

 READ(ubodsin) ioversion
 READ(ubodsin) (header(p),p=1,8)

 ihead(1:nihead)=0
 rhead(1:nrhead)=0
 phead(1:nphead)=0

 if(header(1).GT.nihead.OR.header(2).GT.nrhead.OR. &
                 header(3).GT.nphead.OR.ioversion.GT.2) then
  print*,' warning: file input from newer file format?'
 endif

 READ(ubodsin) (ihead(p),p=1,MIN(nihead,header(1)))
 READ(ubodsin) (rhead(p),p=1,MIN(nrhead,header(2)))
 READ(ubodsin) (phead(p),p=1,MIN(nphead,header(3)))

 nbodies=ihead(1)
 nsph=ihead(2)
 nstar=ihead(3)
 IF(nbodies-nstar-nsph.NE.ihead(4)) THEN
  print*,' inconsistent body count'
  STOP
 ENDIF
 IF(nbodies.GT.nbodsmax.OR.nsph.gt.nsphmax) THEN
  PRINT*,nbodies,nbodsmax,nsph,nsphmax
  PRINT*,' particle overflow '
  PRINT*,'(increase compiled limits)'
  STOP
 ENDIF
 IF(nbodies-nstar-nsph.GT.0.AND.fixthalo) THEN
  PRINT*,'** warning ** rigid halo used and DM particles found'
 ENDIF
 totptag=ihead(5)
 nbh=ihead(6)

 syncflag=ihead(8)
 if(syncflag.NE.0) print*,' ** warning ** async input'

 tnow=rhead(1)
 mtot=rhead(2)
 ektot=rhead(3)
 eptot=rhead(4)
 snheat=rhead(5)
 eradiate=rhead(6)
 esofttot=rhead(7)
 enerror=rhead(8)
 estar=rhead(9)
 efuvheat=rhead(10)
 eradcool=rhead(11)
 tpos=rhead(12)
 teth=rhead(13)
 trad=rhead(14)
 if(meanmwt.LE.0) meanmwt=rhead(15)
 massres=rhead(20)
 ! tstarform=rhead(22)
 tstarform=tnow
 tsnfeedback=tnow
 tbh=tnow
 if(rhead(25).NE.0) eps=rhead(25)

 input(1:nphead)=phead(1:nphead)

 pordercount=pordercount+1
 ppropcount=ppropcount+1

 IF(phead(1).EQ.1) READ(ubodsin) (mass(p),p=1,nbodies)
 IF(phead(2).EQ.1) READ(ubodsin) (pos(p,1),p=1,nbodies)
 IF(phead(2).EQ.1) READ(ubodsin) (pos(p,2),p=1,nbodies)
 IF(phead(2).EQ.1) READ(ubodsin) (pos(p,3),p=1,nbodies)
 IF(phead(3).EQ.1) READ(ubodsin) (vel(p,1),p=1,nbodies)
 IF(phead(3).EQ.1) READ(ubodsin) (vel(p,2),p=1,nbodies)
 IF(phead(3).EQ.1) READ(ubodsin) (vel(p,3),p=1,nbodies)
 IF(phead(4).EQ.1) READ(ubodsin) (epsgrav(p),p=1,nbodies)
 IF(phead(5).EQ.1) READ(ubodsin) (tform(p),p=1,nbodies)
 IF(phead(6).EQ.1) READ(ubodsin) (acc(p,1),p=1,nbodies)
 IF(phead(6).EQ.1) READ(ubodsin) (acc(p,2),p=1,nbodies)
 IF(phead(6).EQ.1) READ(ubodsin) (acc(p,3),p=1,nbodies)

IF(nsph.GT.0) THEN
 IF(phead(10).EQ.1) READ(ubodsin) (rho(p),p=1,nsph)
 IF(phead(11).EQ.1) THEN
  IF(uentropy) THEN
   READ(ubodsin) (tempvect(p),p=1,nsph)
  ELSE
   READ(ubodsin) (ethermal(p),p=1,nsph)
  ENDIF
 ENDIF
 IF(phead(12).EQ.1) THEN
  IF(uentropy) THEN
   READ(ubodsin) (entropy(p),p=1,nsph)
  ELSE
   READ(ubodsin) (tempvect(p),p=1,nsph)
  ENDIF
 ENDIF
 IF(phead(13).EQ.1) READ(ubodsin) (hsmooth(p),p=1,nsph)
 IF(phead(14).EQ.1) READ(ubodsin) (fuvheat(p),p=1,nsph)
 IF(phead(15).EQ.1) READ(ubodsin) (esnthdt(p),p=1,nsph)
 IF(phead(16).EQ.1) READ(ubodsin) (tcollaps(p),p=1,nsph)
 IF(phead(17).EQ.1) READ(ubodsin) (temperat(p),p=1,nsph)
 IF(phead(18).EQ.1) READ(ubodsin) (elecfrac(p),p=1,nsph)
 IF(phead(19).EQ.1) READ(ubodsin) (csound(p),p=1,nsph)
 IF(phead(20).EQ.1) READ(ubodsin) (csound(p),p=1,nsph)
 IF(phead(21).EQ.1) READ(ubodsin) (hsmdivv(p),p=1,nsph)
 IF(phead(22).EQ.1) READ(ubodsin) (mumaxdvh(p),p=1,nsph)
 IF(phead(23).EQ.1) READ(ubodsin) (hsmcurlv(p),p=1,nsph)
 IF(phead(24).EQ.1) READ(ubodsin) (vdisp(p),p=1,nsph)
 IF(phead(25).EQ.1) READ(ubodsin) (h2frac(p),p=1,nsph)
 IF(phead(26).EQ.1) READ(ubodsin) ! ignore
 IF(phead(27).EQ.1) READ(ubodsin) ! ignore
 ENDIF

IF(nstar.GT.0) THEN
 IF(phead(34).EQ.1) READ(ubodsin) (starfuv(p),p=nbodies-nstar+1,nbodies)
 IF(phead(35).EQ.1) READ(ubodsin) (snentropy(p),p=nbodies-nstar+1,nbodies)
ENDIF

 IF(phead(40).EQ.1) READ(ubodsin) (phi(p),p=1,nbodies)
 IF(phead(41).EQ.1) READ(ubodsin) (phiext(p),p=1,nbodies)
 IF(phead(42).EQ.1) READ(ubodsin) (nbexist(p),p=1,nbodies)
 IF(phead(43).EQ.1) READ(ubodsin) (itimestp(p),p=1,nbodies)

 close(ubodsin)

end subroutine readbods

subroutine prepareoutput
 include 'globals.h'
 integer :: nkeep,nbods
 if(uentropy) then
  tempvect(1:nsph)=entropy(1:nsph)/gamma1*rho(1:nsph)**gamma1
  tempvect(nsph+1:2*nsph)=dentdt(1:nsph)/gamma1*rho(1:nsph)**gamma1
 else
  if(.NOT.isotherm) then
   tempvect(1:nsph)=ethermal(1:nsph)*gamma1/rho(1:nsph)**gamma1
   tempvect(nsph+1:2*nsph)=dethdt(1:nsph)*gamma1/rho(1:nsph)**gamma1
  else
   ethermal(1:nsph)=csound(1:nsph)**2
   tempvect(1:nsph)=csound(1:nsph)**2
   tempvect(nsph+1:2*nsph)=0.
   temperat(1:nsph)=csound(1:nsph)**2*meanmwt*mhboltz
  endif
 endif
end subroutine

subroutine itos(i,n,s)
 integer,intent(in) :: i,n
 character(len=n), intent(inout) :: s
 character(len=11) :: nstring
 data nstring/'0123456789X'/
 integer :: j,k,l

 if(i.LT.0.OR.i.GE.10**n) then
  do k=1,n
  s(k:k)=nstring(11:11)
  enddo
  return
 endif
 j=1
 do k=1,n
 l=1+mod(i,10*j)/j
 s(n-k+1:n-k+1)=nstring(l:l)
 j=j*10
 enddo

end subroutine

! insert autmatically generated code below

subroutine writedump(n)
 include 'globals.h'
 integer i,n
 real dummy
 logical ldummy
 open(unit=uboddump,file=dumpfile,status='UNKNOWN',form='UNFORMATTED')

write(uboddump) n

      write(uboddump)  output,firstsnap,stepout,steplog,verbosity,       &
   datadir,inputfile,outputfile,halofile
      write(uboddump)  pboxsize,usesph,radiate,starform,cosmo
      write(uboddump)  unitl_in_kpc,unitm_in_msun
      write(uboddump)  dtime,tstepcrit,tstpcr2,freev,freea,freevexp,   &
   freeaexp,nsteps,max_tbin,minppbin,sqrttstp,acc_tstp, &
   freetstp
      write(uboddump)  bh_tol,eps,gdgtol,nn_tol,targetnn,usequad,         &
   directsum,selfgrav,fixthalo,adaptive_eps,gdgop
      write(uboddump)  epsgas,gamma,alpha,beta,epssph,courant,removgas, &
   consthsm,nsmtol,nsmooth,smoothinput,consph,sphinit,uentropy, &
   isotherm,eps_is_h,mingaseps,balsara,hupdatemethod,sph_visc
      write(uboddump)  graineff,crionrate,heat_par1,heat_par2,cool_par,    &
   optdepth
      write(uboddump)  comove
      write(uboddump)  tcollfac,masscrit,sfeff,tbubble,sne_eff,tsnbeg,   &
   rhomax,sfmode,feedback
      write(uboddump)  rndtable,rnseed
      write(uboddump)  nbodies,nsnap
      write(uboddump)  rmin,rsize,incells,incellsg
      write(uboddump)  root
      write(uboddump)  nttot,ntmin,ntmax,ntavg,nttotfuv,ntminfuv,        &
   ntmaxfuv,ntavgfuv
      write(uboddump)  nstot,nsmin,nsmax,nsavg
      write(uboddump)  tnow,tpos
      write(uboddump)  tiny
      write(uboddump)  mtot,etot,ektot,eptot,mstar,mgas,snheat,esofttot, &
   enerror,amvec,cmpos,cmvel
      write(uboddump)  input
      write(uboddump)  npactive,nsphact
      write(uboddump)  etol,dummy,tsteppos,active_bin
      write(uboddump)  syncflag,entropyflag
      write(uboddump)  teth,ethtot
      write(uboddump)  symmetry
      write(uboddump)  dummy,dummy,gamma1,nstar,nsph
      write(uboddump)  massres
      write(uboddump)  deldr2i
      write(uboddump)  wsmooth,dwsmooth
      write(uboddump)  nntot,nnmin,nnmax,nnavg
      write(uboddump)  hboxsize
      write(uboddump)  poshalo,masshalo
      write(uboddump)  eradiate,trad,meanmwt,     &
   mhboltz,fhydrogn,mumhkbol,mumhkgam, &
   mumhkg1,efuvheat,eradcool,mcold,mluke,mwarm,mhot,estar, &
   heatconst,densconst,timescale,lengthscale,velscale,flxscale, &
   massscale
      write(uboddump)  tstarform,tsnfeedback,snenergy,totptag
      write(uboddump)  h2time
      write(uboddump)  smoothuv
      write(uboddump)  tbh,nbh
      write(uboddump)  tpm,rcut,rcut2,pmsoft,pmdr,pmpot,       &
   pmacc
      write(uboddump)  searchpos,searchh,searchdelta,searchacc4,     &
   searchn,searchreuse,reuseflag,ncalls,nsearches
 write(uboddump) ( itimestp(i),i=1,nbodies)
 write(uboddump) ( nbexist(i),i=1,nbodies)
 write(uboddump) ( order_bodlist(i),i=1,nbodies)
 write(uboddump) ( otimestp(i),i=1,nbodies)
 write(uboddump) ( pactive(i),i=1,nbodies)
 write(uboddump) ( acc(i,1:ndim+1),i=1,nbodies)
 write(uboddump) ( csound(i),i=1,nsph)
 write(uboddump) ( derad(i),i=1,nsph)
 write(uboddump) ( dethdt(i),i=1,nsph)
 write(uboddump) ( dethold(i),i=1,nsph)
 write(uboddump) ( drhodh(i),i=1,nsph)
 write(uboddump) ( elecfrac(i),i=1,nsph)
 write(uboddump) ( epsgrav(i),i=1,nbodies)
 write(uboddump) ( esnthdt(i),i=1,nsph)
 write(uboddump) ( ethermal(i),i=1,nsph)
 write(uboddump) ( ethold(i),i=1,nsph)
 write(uboddump) ( fuvheat(i),i=1,nbodies)
 write(uboddump) ( h2frac(i),i=1,nsph)
 write(uboddump) ( hsmcurlv(i),i=1,nsph)
 write(uboddump) ( hsmdivv(i),i=1,nsph)
 write(uboddump) ( hsmooth(i),i=1,nbodies)
 write(uboddump) ( mass(i),i=1,nbodies)
 write(uboddump) ( mumaxdvh(i),i=1,nsph)
 write(uboddump) ( oldderad(i),i=1,nsph)
 write(uboddump) ( phi(i),i=1,nbodies)
 write(uboddump) ( phiext(i),i=1,nbodies)
 write(uboddump) ( pos(i,1:ndim),i=1,nbodies)
 write(uboddump) ( rho(i),i=1,nsph)
 write(uboddump) ( snentropy(i),i=1,nbodies)
 write(uboddump) ( starfuv(i),i=1,nbodies)
 write(uboddump) ( tcollaps(i),i=1,nsph)
 write(uboddump) ( temperat(i),i=1,nsph)
 write(uboddump) ( tfeedb(i),i=1,nbodies)
 write(uboddump) ( tform(i),i=1,nbodies)
 write(uboddump) ( tvel(i),i=1,nbodies)
 write(uboddump) ( vdisp(i),i=1,nsph)
 write(uboddump) ( vel(i,1:ndim),i=1,nbodies)
 write(uboddump) ( veltpos(i,1:ndim),i=1,nsph)
 close(uboddump)
end subroutine

subroutine readdump(n)
 include 'globals.h'
 integer ioerror,n,i,idummy,idummy2
 real dummy
 logical ldummy
 open(unit=uboddump,file=dumpfile,status='OLD',form='UNFORMATTED',iostat=ioerror)
 if(ioerror.NE.0) then
 print*,' dump file not present'
 stop
 endif

read(uboddump) n

      read(uboddump)  output,firstsnap,stepout,steplog,verbosity,       &
   datadir,inputfile,outputfile,halofile
      read(uboddump)  pboxsize,usesph,radiate,starform,cosmo
      read(uboddump)  unitl_in_kpc,unitm_in_msun
      read(uboddump)  dtime,tstepcrit,tstpcr2,freev,freea,freevexp,   &
   freeaexp,nsteps,max_tbin,minppbin,sqrttstp,acc_tstp, &
   freetstp
      read(uboddump)  bh_tol,eps,gdgtol,nn_tol,targetnn,usequad,         &
   directsum,selfgrav,fixthalo,adaptive_eps,gdgop
      read(uboddump)  epsgas,gamma,alpha,beta,epssph,courant,removgas, &
   consthsm,nsmtol,nsmooth,smoothinput,consph,sphinit,uentropy, &
   isotherm,eps_is_h,mingaseps,balsara,hupdatemethod,sph_visc
      read(uboddump)  graineff,crionrate,heat_par1,heat_par2,cool_par,    &
   optdepth
      read(uboddump)  comove
      read(uboddump)  tcollfac,masscrit,sfeff,tbubble,sne_eff,tsnbeg,   &
   rhomax,sfmode,feedback
      read(uboddump)  rndtable,rnseed
      read(uboddump)  nbodies,nsnap
      read(uboddump)  rmin,rsize,incells,incellsg
      read(uboddump)  root
      read(uboddump)  nttot,ntmin,ntmax,ntavg,nttotfuv,ntminfuv,        &
   ntmaxfuv,ntavgfuv
      read(uboddump)  nstot,nsmin,nsmax,nsavg
      read(uboddump)  tnow,tpos
      read(uboddump)  tiny
      read(uboddump)  mtot,etot,ektot,eptot,mstar,mgas,snheat,esofttot, &
   enerror,amvec,cmpos,cmvel
      read(uboddump)  input
      read(uboddump)  npactive,nsphact
      read(uboddump)  etol,dummy,tsteppos,active_bin
      read(uboddump)  syncflag,entropyflag
      read(uboddump)  teth,ethtot
      read(uboddump)  symmetry
      read(uboddump)  dummy,dummy,gamma1,nstar,nsph
      read(uboddump)  massres
      read(uboddump)  deldr2i
      read(uboddump)  wsmooth,dwsmooth
      read(uboddump)  nntot,nnmin,nnmax,nnavg
      read(uboddump)  hboxsize
      read(uboddump)  poshalo,masshalo
      read(uboddump)  eradiate,trad,meanmwt,     &
   mhboltz,fhydrogn,mumhkbol,mumhkgam, &
   mumhkg1,efuvheat,eradcool,mcold,mluke,mwarm,mhot,estar, &
   heatconst,densconst,timescale,lengthscale,velscale,flxscale, &
   massscale
      read(uboddump)  tstarform,tsnfeedback,snenergy,totptag
      read(uboddump)  h2time
      read(uboddump)  smoothuv
      read(uboddump)  tbh,nbh
      read(uboddump)  tpm,rcut,rcut2,pmsoft,pmdr,pmpot,       &
   pmacc
      read(uboddump)  searchpos,searchh,searchdelta,searchacc4,     &
   searchn,searchreuse,reuseflag,ncalls,nsearches
 read(uboddump) ( itimestp(i),i=1,nbodies)
 read(uboddump) ( nbexist(i),i=1,nbodies)
 read(uboddump) ( order_bodlist(i),i=1,nbodies)
 read(uboddump) ( otimestp(i),i=1,nbodies)
 read(uboddump) ( pactive(i),i=1,nbodies)
 read(uboddump) ( acc(i,1:ndim+1),i=1,nbodies)
 read(uboddump) ( csound(i),i=1,nsph)
 read(uboddump) ( derad(i),i=1,nsph)
 read(uboddump) ( dethdt(i),i=1,nsph)
 read(uboddump) ( dethold(i),i=1,nsph)
 read(uboddump) ( drhodh(i),i=1,nsph)
 read(uboddump) ( elecfrac(i),i=1,nsph)
 read(uboddump) ( epsgrav(i),i=1,nbodies)
 read(uboddump) ( esnthdt(i),i=1,nsph)
 read(uboddump) ( ethermal(i),i=1,nsph)
 read(uboddump) ( ethold(i),i=1,nsph)
 read(uboddump) ( fuvheat(i),i=1,nbodies)
 read(uboddump) ( h2frac(i),i=1,nsph)
 read(uboddump) ( hsmcurlv(i),i=1,nsph)
 read(uboddump) ( hsmdivv(i),i=1,nsph)
 read(uboddump) ( hsmooth(i),i=1,nbodies)
 read(uboddump) ( mass(i),i=1,nbodies)
 read(uboddump) ( mumaxdvh(i),i=1,nsph)
 read(uboddump) ( oldderad(i),i=1,nsph)
 read(uboddump) ( phi(i),i=1,nbodies)
 read(uboddump) ( phiext(i),i=1,nbodies)
 read(uboddump) ( pos(i,1:ndim),i=1,nbodies)
 read(uboddump) ( rho(i),i=1,nsph)
 read(uboddump) ( snentropy(i),i=1,nbodies)
 read(uboddump) ( starfuv(i),i=1,nbodies)
 read(uboddump) ( tcollaps(i),i=1,nsph)
 read(uboddump) ( temperat(i),i=1,nsph)
 read(uboddump) ( tfeedb(i),i=1,nbodies)
 read(uboddump) ( tform(i),i=1,nbodies)
 read(uboddump) ( tvel(i),i=1,nbodies)
 read(uboddump) ( vdisp(i),i=1,nsph)
 read(uboddump) ( vel(i,1:ndim),i=1,nbodies)
 read(uboddump) ( veltpos(i,1:ndim),i=1,nsph)
 close(uboddump)
end subroutine

max_stars_repo_name: joshuawall/amuse

max_issues_repo_name: joshuawall/amuse

max_forks_repo_name: joshuawall/amuse


text:
!***********************************************************************
!  Integrated Water Flow Model (IWFM)
!  Copyright (C) 2005-2021
!  State of California, Department of Water Resources
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU General Public License
!  as published by the Free Software Foundation; either version 2
!  of the License, or (at your option) any later version.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  (http://www.gnu.org/copyleft/gpl.html)
!
!  You should have received a copy of the GNU General Public License
!  along with this program; if not, write to the Free Software
!  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
!
!  For tecnical support, e-mail: IWFMtechsupport@water.ca.gov
!***********************************************************************
MODULE Class_Version
  USE GeneralUtilities  ,  ONLY : FirstLocation           , &
                                  CleanSpecialCharacters
  USE MessageLogger     ,  ONLY : SetLastMessage          , &
                                  MessageArray            , &
                                  iFatal
  USE IOInterface       ,  ONLY : GenericFileType
  IMPLICIT NONE

! ******************************************************************
! ******************************************************************
! ******************************************************************
! ***
! *** VARIABLE DEFINITIONS
! ***
! ******************************************************************
! ******************************************************************
! ******************************************************************

  ! -------------------------------------------------------------
  ! --- PUBLIC ENTITIES
  ! -------------------------------------------------------------
    PRIVATE
    PUBLIC :: VersionType  , &
              ReadVersion

  ! -------------------------------------------------------------
  ! --- VERSION DATA TYPE
  ! -------------------------------------------------------------
  TYPE VersionType
    PRIVATE
    INTEGER                  :: iLenVersion = 0
    CHARACTER(:),ALLOCATABLE :: cVersion
  CONTAINS
    GENERIC,PUBLIC          :: New  => Version_New_FromComponents  , &
                                       Version_New_FromFullString
    PROCEDURE,NOPASS,PUBLIC :: Version_New_FromComponents
    PROCEDURE,NOPASS,PUBLIC :: Version_New_FromFullString
    PROCEDURE,PASS,PUBLIC   :: Kill
    PROCEDURE,PASS,PUBLIC   :: GetVersion
    PROCEDURE,PASS,PUBLIC   :: IsDefined
  END TYPE VersionType

  ! -------------------------------------------------------------
  ! --- MISC. METHODS
  ! -------------------------------------------------------------
  INTEGER,PARAMETER                   :: ModNameLen = 15
  CHARACTER(LEN=ModNameLen),PARAMETER :: ModName    = 'Class_Version::'

CONTAINS

! ******************************************************************
! ******************************************************************
! ******************************************************************
! ***
! *** METHODS
! ***
! ******************************************************************
! ******************************************************************
! ******************************************************************

  ! -------------------------------------------------------------
  ! --- CONSTRUCTOR 1
  ! -------------------------------------------------------------
  FUNCTION Version_New_FromComponents(iLenVersion,cVersion,cRevision) RESULT(Version)
    INTEGER,INTENT(IN)                    :: iLenVersion
    CHARACTER(LEN=iLenVersion),INTENT(IN) :: cVersion
    CHARACTER(LEN=*),INTENT(IN)           :: cRevision
    TYPE(VersionType)                     :: Version

    !Local variables
    INTEGER :: iLenRev

    !Initialize the version
    Version%iLenVersion = iLenVersion
    ALLOCATE(CHARACTER(iLenVersion) :: Version%cVersion)
    Version%cVersion = cVersion

    !Include the revision number in the version
    iLenRev                                             = LEN_TRIM(ADJUSTL(cRevision))
    Version%cVersion(iLenVersion-iLenRev+1:iLenVersion) = TRIM(ADJUSTL(cRevision))

  END FUNCTION Version_New_FromComponents

  ! -------------------------------------------------------------
  ! --- CONSTRUCTOR 2
  ! -------------------------------------------------------------
  FUNCTION Version_New_FromFullString(cVrs) RESULT(Version)
    CHARACTER(LEN=*),INTENT(IN) :: cVrs
    TYPE(VersionType)           :: Version

    !Local variables
    INTEGER :: iLenVersion

    iLenVersion         = LEN_TRIM(cVrs)
    Version%iLenVersion = iLenVersion
    ALLOCATE (CHARACTER(iLenVersion) :: Version%cVersion)
    Version%cVersion = cvrs

  END FUNCTION Version_New_FromFullString

  ! -------------------------------------------------------------
  ! --- DESTRUCTOR
  ! -------------------------------------------------------------
  SUBROUTINE Kill(Version)
    CLASS(VersionType) :: Version

    !Local variables
    INTEGER           :: ErrorCode

    DEALLOCATE (Version%cVersion , STAT=ErrorCode)
    Version%iLenVersion = 0

  END SUBROUTINE Kill

  ! -------------------------------------------------------------
  ! --- GETTER
  ! -------------------------------------------------------------
  FUNCTION GetVersion(Version) RESULT(cVrs)
    CLASS(VersionType)       :: Version
    CHARACTER(:),ALLOCATABLE :: cVrs

    ALLOCATE( CHARACTER(LEN=Version%iLenVersion) :: cVrs )
    cVrs = Version%cVersion

  END FUNCTION GetVersion

  ! -------------------------------------------------------------
  ! --- PREDICATE
  ! -------------------------------------------------------------
  FUNCTION IsDefined(Version) RESULT(lDefined)
    CLASS(VersionType),INTENT(IN) :: Version
    LOGICAL                       :: lDefined

    IF (Version%iLenVersion .EQ. 0) THEN
        lDefined = .FALSE.
    ELSE
        lDefined = .TRUE.
    END IF

  END FUNCTION IsDefined

  ! -------------------------------------------------------------
  ! --- READ VERSION FROM A FILE
  ! -------------------------------------------------------------
  SUBROUTINE ReadVersion(InFile,cComponent,cVersion,iStat)
    TYPE(GenericFileType)       :: InFile
    CHARACTER(LEN=*),INTENT(IN) :: cComponent
    CHARACTER(:),ALLOCATABLE    :: cVersion
    INTEGER,INTENT(OUT)         :: iStat

    !Local variables
    CHARACTER(LEN=ModNameLen+11) :: ThisProcedure = ModName // 'ReadVersion'
    CHARACTER                    :: ALine*2000
    INTEGER                      :: iStart

    !Initialize
    iStat = 0

    CALL InFile%ReadData(ALine,iStat)  ;  IF (iStat .EQ. -1) RETURN
    CALL CleanSpecialCharacters(ALine)  ;  ALine = ADJUSTL(ALine)
    iStart   = FirstLocation('#',ALine,Back=.FALSE.)
    IF (iStart .LE. 0) THEN
        MessageArray(1) = 'Error in identifying the version number of the '//TRIM(cComponent)//' component!'
        MessageArray(2) = 'Make sure that the version number is listed at the first line of the'
        MessageArray(3) = TRIM(cComponent)//' input file (see the input file template for format).'
        CALL SetLastMessage(MessageArray(1:3),iFatal,ThisProcedure)
        iStat = -1
        RETURN
    END IF
    cVersion = ALine(iStart+1:LEN_TRIM(ALine))

  END SUBROUTINE ReadVersion

END MODULE

max_stars_repo_name: giorgk/IWFM-SAFE

max_issues_repo_name: giorgk/IWFM-SAFE

max_forks_repo_name: giorgk/IWFM-SAFE


text:
module rate_type_module

  use actual_network, only : nrates, num_rate_groups

  use amrex_fort_module, only : rt => amrex_real
  implicit none

  type :: rate_t
     ! the temperature at which the rates were evaluated
     real(rt)         :: T_eval
     real(rt)         :: rates(num_rate_groups, nrates)
  end type rate_t

end module rate_type_module

max_stars_repo_name: doreenfan/Microphysics

max_issues_repo_name: doreenfan/Microphysics

max_forks_repo_name: doreenfan/Microphysics


text:
module maxvol_lib
use time_lib
use ort_lib, only : ort0
use mat_lib, only : matinv
use sort_lib
use say_lib

double precision,private,parameter :: tol_def=1.d-2
logical,private,parameter :: loud=.false.
interface maxvol
 module procedure maxvol_d,maxvol_z
end interface
interface maxvola
 module procedure maxvola_d,maxvola_z
end interface
interface volume
 module procedure volume_d,volume_z
end interface
interface domval
 module procedure domval_d,domval_z
end interface

private :: maxvol_d,volume_d,maxvol0_d
private :: maxvol_z,volume_z,maxvol0_z

contains
 subroutine maxvol0_d(a,p,z,tol)
  implicit none
  double precision,intent(in) :: a(:,:)
  integer,intent(inout) :: p(size(a,2))
  double precision,intent(inout) :: z(size(a,1),size(a,2))
  double precision,intent(in),optional  :: tol
  character(len=*),parameter :: subnam='maxvol'
  integer,parameter          :: maxiter=2**10
  integer          :: m,n,i,j,ip,ij(2),info,iter,lwork,ii,jj
  integer,allocatable         :: q(:),per(:)
  double precision :: zinv,vol,amax,err,nrm,eps,t0
  double precision,allocatable :: b(:,:),work(:),u(:),v(:)
  double precision,allocatable :: bb(:,:),tt(:,:) !debug
  character(len=4) :: verb
  double precision,external :: dnrm2
  integer,external :: idamax

  m=size(a,1); n=size(a,2); t0=timef()

  if(present(tol))then; eps=tol;else;eps=tol_def;endif
  lwork=256*n
  allocate(work(lwork),b(n,n),q(n),per(m),u(m),v(n))

  if(isperm(p,m))then
   verb="from"
   vol=volume(a(p(1:n),:))
  else
   verb="init"
   call dcopy(m*n,a,1,z,1)   ! z=a
   call dgetrf(m,n,z,m,p,info)
   if(info.ne.0)then; write(*,*)subnam,': dgetrf(mxn) info: ',info; stop; end if
   forall(i=1:n)v(i)=dlog(dabs(z(i,i))); vol=sum(v)
   q=(/ (i,i=1,n) /)
   call dgetri(n,z,m,q,work,lwork,info)
   if(info.ne.0)then; write(*,*)subnam,': dgetri(mxn) info: ',info; stop; end if
   b=z(1:n,:)
   per(1:m)=(/ (i,i=1,m) /)
   do i=1,n
    ip=per(i); per(i)=per(p(i)); per(p(i))=ip ! per(i) <-> per(p(i))
   end do
   p=per(1:n)
   call dgemm('n','n',m,n,n,1.d0,a,m,b,n,0.d0,z,m)

   !write(*,'(2a,i8,1x,i3,3a,f9.3,a,f7.1,a)') subnam, '[',shape(a),'] ',verb,' e',vol,' in',timef()-t0,' sec'
  end if

  do iter=1,maxiter
!--
!   allocate(tt(m,n)); call maxvola_d(a,p,tt); nrm=dnrm2(m*n,tt,1); tt=tt-z; err=dnrm2(m*n,tt,1); deallocate(tt)
!   write(*,'(a,i4,3(a,e9.3))')'iter: ',iter,' err: ',err,' nrm: ',nrm,' relerr: ',err/nrm
!--

   do jj=1,n
    q(jj)=idamax(m,z(1,jj),1)
    v(jj)=abs(z(q(jj),jj))
   end do
   j=maxloc(v,1); i=q(j)

   amax=dabs(z(i,j))
   if(amax.lt.1.d0+eps .or. p(j).eq.i)goto 100
   ip=p(j); p(j)=i
   !write(*,'(a,i4,a,e11.5,a,i7,a,i7)') 'iter: ',iter,' amax: ',amax,' swap: ',ip,' <-> ',i
   u=z(:,j); v=z(ip,:)-z(i,:)
   zinv=1.d0/z(i,j)

!!  call dger(m,n,zinv,u,1,v,1,z,m)

   do jj=1,n
    call daxpy(m,zinv*v(jj),u,1,z(1,jj),1)
   end do
   z(i,:)=0.d0; z(i,j)=1.d0   ! FORCE
  end do
100 continue
  deallocate(work,b,q,per,u,v)

  if(loud)write(*,'(2a,i8,1x,i3,3a,f9.3,a,f9.3,a,i5,a,f7.1,a)') subnam, '[',shape(a),&
       '] ',verb,' e',vol,' upto e',volume(a(p,:)),' in ',iter,' iters in',timef()-t0,' sec'

!  call say(z)
!  call sort(p)
!   write(*,'(a,64(i7,1x))') 'after: ',p
!  do i=1,n
!   z(p(i),:)=0.d0; z(p(i),i)=1.d0
!  end do
!  call say(z)
!--
!allocate(tt(m,n)); call maxvola_d(a,p,tt); nrm=dnrm2(m*n,tt,1); tt=tt-z; err=dnrm2(m*n,tt,1); deallocate(tt)
!  write(*,'(2a,i9,1x,i4,3a,f10.5,a,f10.5,a,i5,a,f9.1,a,e9.3)') subnam, '[',shape(a),&
!       ']:',verb,' e',vol,' upto e',volume(a(p,:)),' in ',iter,' iters in',timef()-t0,' sec; relerr: ',err/nrm
!--
  return
 end subroutine

 subroutine maxvol0_z(a,p,z,tol)
  implicit none
  complex(8),intent(in) :: a(:,:)
  integer,intent(inout) :: p(size(a,2))
  complex(8),intent(inout) :: z(size(a,1),size(a,2))
  double precision,intent(in),optional  :: tol
  character(len=*),parameter :: subnam='maxvol'
  integer,parameter          :: maxiter=2**10
  integer          :: m,n,i,j,ip,ij(2),info,iter,lwork,ii,jj
  integer,allocatable         :: q(:),per(:)
  double precision :: vol,amax,err,nrm,eps,t0
  complex(8)   :: zinv
  complex(8),allocatable :: b(:,:),work(:),u(:)
  double precision,allocatable :: v(:)
  complex(8),allocatable :: bb(:,:),tt(:,:) !debug
  character(len=4) :: verb
  double precision,external :: dznrm2
  integer,external :: izamax
  real(8) abs
  m=size(a,1); n=size(a,2); t0=timef()

  if(present(tol))then; eps=tol;else;eps=tol_def;endif
  lwork=256*n
  allocate(work(lwork),b(n,n),q(n),per(m),u(m),v(n))

  if(isperm(p,m))then
   verb="from"
   vol=volume(a(p(1:n),:))
  else
   verb="init"
   call zcopy(m*n,a,1,z,1)   ! z=a
   call zgetrf(m,n,z,m,p,info)
   if(info.ne.0)then; write(*,*)subnam,': zgetrf(mxn) info: ',info; stop; end if
   forall(i=1:n)v(i)=dlog(abs(z(i,i))); vol=sum(v)
   q=(/ (i,i=1,n) /)
   call zgetri(n,z,m,q,work,lwork,info)
   if(info.ne.0)then; write(*,*)subnam,': zgetri(mxn) info: ',info; stop; end if
   b=z(1:n,:)
   per(1:m)=(/ (i,i=1,m) /)
   do i=1,n
    ip=per(i); per(i)=per(p(i)); per(p(i))=ip ! per(i) <-> per(p(i))
   end do
   p=per(1:n)
   call zgemm('n','n',m,n,n,(1.d0,0.d0),a,m,b,n,(0.d0,0.d0),z,m)

   !write(*,'(2a,i8,1x,i3,3a,f9.3,a,f7.1,a)') subnam, '[',shape(a),'] ',verb,' e',vol,' in',timef()-t0,' sec'
  end if

  do iter=1,maxiter
   do jj=1,n
    q(jj)=izamax(m,z(1,jj),1)
    v(jj)=abs(z(q(jj),jj))
   end do
   j=maxloc(v,1); i=q(j)

   amax=abs(z(i,j))
   if(amax.lt.1.d0+eps .or. p(j).eq.i)goto 100
   ip=p(j); p(j)=i
   !write(*,'(a,i4,a,e11.5,a,i7,a,i7)') 'iter: ',iter,' amax: ',amax,' swap: ',ip,' <-> ',i
   u=z(:,j); v=z(ip,:)-z(i,:)
   zinv=(1.d0,0.d0)/z(i,j)

   do jj=1,n
    call zaxpy(m,zinv*v(jj),u,1,z(1,jj),1)
   end do
   z(i,:)=(0.d0,0.d0); z(i,j)=(1.d0,0.d0)   ! FORCE
  end do
100 continue
  deallocate(work,b,q,per,u,v)

  if(loud)write(*,'(2a,i8,1x,i3,3a,f9.3,a,f9.3,a,i5,a,f7.1,a)') subnam, '[',shape(a),&
       '] ',verb,' e',vol,' upto e',volume(a(p,:)),' in ',iter,' iters in',timef()-t0,' sec'
 end subroutine

 subroutine maxvol_d(a,p,z,tol)
  implicit none
  double precision, intent(in)            :: a(:,:)
  integer,intent(inout)                   :: p(size(a,2))
  double precision,pointer,optional       :: z(:,:)
  double precision,intent(in),optional    :: tol
  character(len=*),parameter :: subnam='maxvol'
  integer :: m,n,k
  double precision,allocatable :: d(:,:),y(:,:),f(:,:),g(:,:),dd(:,:)
  double precision,pointer     :: znew(:,:)
  double precision :: err,nrm
  double precision,external :: dnrm2

  if(size(a,1).lt.size(a,2))then
   write(*,'(2a,a,2i8,a)') subnam,': wide matrices not permited: ','a[',shape(a),']'
   stop
  end if

  if(.not.present(z))then
   allocate(y(size(a,1),size(a,2)))
   if(isperm(p,size(a,1)))call maxvola_d(a,p,y)
   call maxvol0_d(a,p,y,tol)
   deallocate(y)
   return
  end if
  if(.not.associated(z))then
   allocate(z(size(a,1),size(a,2)))
   if(isperm(p,size(a,1)))call maxvola_d(a,p,z)
   call maxvol0_d(a,p,z,tol)
   return
  end if

  if(size(z,1).ne.size(a,1))then
   write(*,'(2a,2(a,2i8),a)') subnam,': sizes mismatch: ','a[',shape(a),'], z[',shape(z),']'
   stop
  end if
  if(size(z,2).gt.size(a,2))then
   write(*,'(2a,2(a,2i8),a)') subnam,': z is wider than a:','a[',shape(a),'], z[',shape(z),']'
   stop
  end if

  m=size(a,1); k=size(a,2)-size(z,2); n=size(z,2)
  if(k.eq.0)then
   if(isperm(p,size(a,1)))call maxvola_d(a,p,z)
   call maxvol0_d(a,p,z,tol)
   return
  end if

  IF(isperm(p(1:n),m))THEN
   allocate(znew(m,n+k),d(m,k),y(m,k),f(n,k),g(k,n))
   d=a(:,n+1:n+k); f=a(p(1:n),n+1:n+k)
   call dgemm('n','n',m,k,n,-1.d0,z,m,f,n,1.d0,d,m)

   d(p(1:n),:)=0.d0     ! FORCE
   call ort0(d)
   p(n+1:n+k)=0; call maxvol0_d(d,p(n+1:n+k),y,tol)
   y(p(1:n),:)=0.d0     ! FORCE

   g=z(p(n+1:n+k),:)
   call dcopy(m*n,z,1,znew,1)
   call dcopy(m*k,y,1,znew(:,n+1:n+k),1)
   call dgemm('n','n',m,n,k,-1.d0,y,m,g,k,1.d0,znew,m)
   znew(p(n+1:n+k),1:n)=0.d0     ! FORCE
   deallocate(d,y,f,g,z)
!--
!   allocate(d(n+k,n+k),dd(n+k,n+k),y(m,n+k))
!   d=a(p(1:n+k),:); write(*,'(2(a,f10.5))')'before: vol: ',volume(d),' val: ',domval(a,p)
!   call matinv(d,dd); call dgemm('n','n',m,n+k,n+k,1.d0,a,m,dd,n+k,0.d0,y,m)
!   y=y-ab; err=dnrm2(m*(n+k),y,1); nrm=dnrm2(m*(n+k),ab,1)
!   write(*,'(3(a,e9.3))')'before: err: ',err,' nrm: ',nrm,' relerr: ',err/nrm
!   deallocate(d,dd,y)
!--
   call maxvol0_d(a,p,znew,tol)
   z=>znew
!--
!   allocate(d(n+k,n+k),dd(n+k,n+k),y(m,n+k))
!   d=a(p(1:n+k),:); write(*,'(2(a,f10.5))')'after: vol: ',volume(d),' val: ',domval(a,p)
!   call matinv(d,dd); call dgemm('n','n',m,n+k,n+k,1.d0,a,m,dd,n+k,0.d0,y,m)
!   y=y-ab; err=dnrm2(m*(n+k),y,1); nrm=dnrm2(m*(n+k),ab,1)
!   write(*,'(3(a,e9.3))')'after: err: ',err,' nrm: ',nrm,' relerr: ',err/nrm
!   deallocate(d,dd,y)
!--
  ELSE
   write(*,*)subnam,': reject invalid p on input: ',p(1:n)
   deallocate(z)
   allocate(z(m,n+k))
   p=0; call maxvol0_d(a,p,z,tol)
  ENDIF

  return
 end subroutine

 subroutine maxvol_z(a,p,z,tol)
  implicit none
  complex(8), intent(in)            :: a(:,:)
  integer,intent(inout)                   :: p(size(a,2))
  complex(8),pointer,optional       :: z(:,:)
  double precision,intent(in),optional    :: tol
  character(len=*),parameter :: subnam='maxvol'
  integer :: m,n,k
  complex(8),allocatable :: d(:,:),y(:,:),f(:,:),g(:,:),dd(:,:)
  complex(8),pointer     :: znew(:,:)
  double precision :: err,nrm
  double precision,external :: dznrm2

  if(size(a,1).lt.size(a,2))then
   write(*,'(2a,a,2i8,a)') subnam,': wide matrices not permited: ','a[',shape(a),']'
   stop
  end if

  if(.not.present(z))then
   allocate(y(size(a,1),size(a,2)))
   if(isperm(p,size(a,1)))call maxvola_z(a,p,y)
   call maxvol0_z(a,p,y,tol)
   deallocate(y)
   return
  end if
  if(.not.associated(z))then
   allocate(z(size(a,1),size(a,2)))
   if(isperm(p,size(a,1)))call maxvola_z(a,p,z)
   call maxvol0_z(a,p,z,tol)
   return
  end if

  if(size(z,1).ne.size(a,1))then
   write(*,'(2a,2(a,2i8),a)') subnam,': sizes mismatch: ','a[',shape(a),'], z[',shape(z),']'
   stop
  end if
  if(size(z,2).gt.size(a,2))then
   write(*,'(2a,2(a,2i8),a)') subnam,': z is wider than a:','a[',shape(a),'], z[',shape(z),']'
   stop
  end if

  m=size(a,1); k=size(a,2)-size(z,2); n=size(z,2)
  if(k.eq.0)then
   if(isperm(p,size(a,1)))call maxvola_z(a,p,z)
   call maxvol0_z(a,p,z,tol)
   return
  end if

  IF(isperm(p(1:n),m))THEN
   allocate(znew(m,n+k),d(m,k),y(m,k),f(n,k),g(k,n))
   d=a(:,n+1:n+k); f=a(p(1:n),n+1:n+k)
   call zgemm('n','n',m,k,n,-(1.d0,0.d0),z,m,f,n,(1.d0,0.d0),d,m)

   d(p(1:n),:)=(0.d0,0.d0)     ! FORCE
   call ort0(d)
   p(n+1:n+k)=0; call maxvol0_z(d,p(n+1:n+k),y,tol)
   y(p(1:n),:)=(0.d0,0.d0)     ! FORCE

   g=z(p(n+1:n+k),:)
   call zcopy(m*n,z,1,znew,1)
   call zcopy(m*k,y,1,znew(:,n+1:n+k),1)
   call zgemm('n','n',m,n,k,-(1.d0,0.d0),y,m,g,k,(1.d0,0.d0),znew,m)
   znew(p(n+1:n+k),1:n)=(0.d0,0.d0)     ! FORCE
   deallocate(d,y,f,g,z)

   call maxvol0_z(a,p,znew,tol)
   z=>znew

  ELSE
   write(*,*)subnam,': reject invalid p on input: ',p(1:n)
   deallocate(z)
   allocate(z(m,n+k))
   p=0; call maxvol0_z(a,p,z,tol)
  ENDIF
 end subroutine

 double precision function volume_d(a) result(vol)
  implicit none
  double precision,intent(in) :: a(:,:)
  integer :: m,n,info,i
  double precision,allocatable :: b(:,:),d(:)
  integer, allocatable         :: p(:)

  m=size(a,1); n=size(a,2)
  if(m.ne.n)then
   vol=-999.d0
   return
  end if
  allocate(b(m,n),d(n),p(n))
  b=a
  call dgetrf(n,n,b,n,p,info)
  if(info.ne.0)then
   vol=-888.d0
   return
  end if
  forall(i=1:n) d(i)=dlog(dabs(b(i,i)))
  vol=sum(d)
  deallocate(b,d,p)
  return
 end function

 double precision function volume_z(a) result(vol)
  implicit none
  complex(8),intent(in) :: a(:,:)
  integer :: m,n,info,i
  complex(8),allocatable :: b(:,:),d(:)
  integer, allocatable         :: p(:)

  m=size(a,1); n=size(a,2)
  if(m.ne.n)then
   vol=-999.d0
   return
  end if
  allocate(b(m,n),d(n),p(n))
  b=a
  call zgetrf(n,n,b,n,p,info)
  if(info.ne.0)then
   vol=-888.d0
   return
  end if
  forall(i=1:n) d(i)=dlog(abs(b(i,i)))
  vol=sum(d)
  deallocate(b,d,p)
  return
 end function

 double precision function domval_d(a,p) result(val)
  implicit none
  double precision,intent(in) :: a(:,:)
  integer,intent(in)          :: p(:)
  double precision,allocatable:: b(:,:)
  integer                     :: m,n
  m=size(a,1); n=size(a,2)
  if(size(p).ne.n)then; val=-999.d0; return; end if
  allocate(b(m,n))
  call maxvola(a,p,b)
  val=maxval(abs(b))
  deallocate(b)
  return
 end function
 double precision function domval_z(a,p) result(val)
  implicit none
  complex(8),intent(in) :: a(:,:)
  integer,intent(in)          :: p(:)
  complex(8),allocatable:: b(:,:)
  integer                     :: m,n
  m=size(a,1); n=size(a,2)
  if(size(p).ne.n)then; val=-999.d0; return; end if
  allocate(b(m,n))
  call maxvola(a,p,b)
  val=maxval(abs(b))
  deallocate(b)
  return
 end function

 subroutine maxvola_d(a,p,b)
  implicit none
  double precision,intent(in) :: a(:,:)
  integer,intent(in)          :: p(size(a,2))
  double precision,intent(out):: b(size(a,1),size(a,2))
  double precision,allocatable:: abox(:,:),ainv(:,:)
  integer                     :: m,n
  m=size(a,1); n=size(a,2)
  allocate(abox(n,n),ainv(n,n))
  abox=a(p,:)
  call matinv(abox,ainv,alg='trf')
  call dgemm('n','n',m,n,n,1.d0,a,m,ainv,n,0.d0,b,m)
  deallocate(abox,ainv)
 end subroutine
 subroutine maxvola_z(a,p,b)
  implicit none
  complex(8),intent(in) :: a(:,:)
  integer,intent(in)          :: p(size(a,2))
  complex(8),intent(out):: b(size(a,1),size(a,2))
  complex(8),allocatable:: abox(:,:),ainv(:,:)
  integer                     :: m,n
  m=size(a,1); n=size(a,2)
  allocate(abox(n,n),ainv(n,n))
  abox=a(p,:)
  call matinv(abox,ainv,alg='trf')
  call zgemm('n','n',m,n,n,(1.d0,0.d0),a,m,ainv,n,(0.d0,0.d0),b,m)
  deallocate(abox,ainv)
 end subroutine

end module

max_stars_repo_name: dimpase/tt-fort

max_issues_repo_name: dimpase/tt-fort

max_forks_repo_name: dimpase/tt-fort


text:
      SUBROUTINE CGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W,
     $                   VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK,
     $                   BWORK, INFO )
*
*  -- LAPACK driver routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1999
*
*     .. Scalar Arguments ..
      CHARACTER          JOBVS, SENSE, SORT
      INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
      REAL               RCONDE, RCONDV
*     ..
*     .. Array Arguments ..
      LOGICAL            BWORK( * )
      REAL               RWORK( * )
      COMPLEX            A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
*     ..
*     .. Function Arguments ..
      LOGICAL            SELECT
      EXTERNAL           SELECT
*     ..
*
*  Purpose
*  =======
*
*  CGEESX computes for an N-by-N complex nonsymmetric matrix A, the
*  eigenvalues, the Schur form T, and, optionally, the matrix of Schur
*  vectors Z.  This gives the Schur factorization A = Z*T*(Z**H).
*
*  Optionally, it also orders the eigenvalues on the diagonal of the
*  Schur form so that selected eigenvalues are at the top left;
*  computes a reciprocal condition number for the average of the
*  selected eigenvalues (RCONDE); and computes a reciprocal condition
*  number for the right invariant subspace corresponding to the
*  selected eigenvalues (RCONDV).  The leading columns of Z form an
*  orthonormal basis for this invariant subspace.
*
*  For further explanation of the reciprocal condition numbers RCONDE
*  and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
*  these quantities are called s and sep respectively).
*
*  A complex matrix is in Schur form if it is upper triangular.
*
*  Arguments
*  =========
*
*  JOBVS   (input) CHARACTER*1
*          = 'N': Schur vectors are not computed;
*          = 'V': Schur vectors are computed.
*
*  SORT    (input) CHARACTER*1
*          Specifies whether or not to order the eigenvalues on the
*          diagonal of the Schur form.
*          = 'N': Eigenvalues are not ordered;
*          = 'S': Eigenvalues are ordered (see SELECT).
*
*  SELECT  (input) LOGICAL FUNCTION of one COMPLEX argument
*          SELECT must be declared EXTERNAL in the calling subroutine.
*          If SORT = 'S', SELECT is used to select eigenvalues to order
*          to the top left of the Schur form.
*          If SORT = 'N', SELECT is not referenced.
*          An eigenvalue W(j) is selected if SELECT(W(j)) is true.
*
*  SENSE   (input) CHARACTER*1
*          Determines which reciprocal condition numbers are computed.
*          = 'N': None are computed;
*          = 'E': Computed for average of selected eigenvalues only;
*          = 'V': Computed for selected right invariant subspace only;
*          = 'B': Computed for both.
*          If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
*
*  N       (input) INTEGER
*          The order of the matrix A. N >= 0.
*
*  A       (input/output) COMPLEX array, dimension (LDA, N)
*          On entry, the N-by-N matrix A.
*          On exit, A is overwritten by its Schur form T.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*
*  SDIM    (output) INTEGER
*          If SORT = 'N', SDIM = 0.
*          If SORT = 'S', SDIM = number of eigenvalues for which
*                         SELECT is true.
*
*  W       (output) COMPLEX array, dimension (N)
*          W contains the computed eigenvalues, in the same order
*          that they appear on the diagonal of the output Schur form T.
*
*  VS      (output) COMPLEX array, dimension (LDVS,N)
*          If JOBVS = 'V', VS contains the unitary matrix Z of Schur
*          vectors.
*          If JOBVS = 'N', VS is not referenced.
*
*  LDVS    (input) INTEGER
*          The leading dimension of the array VS.  LDVS >= 1, and if
*          JOBVS = 'V', LDVS >= N.
*
*  RCONDE  (output) REAL
*          If SENSE = 'E' or 'B', RCONDE contains the reciprocal
*          condition number for the average of the selected eigenvalues.
*          Not referenced if SENSE = 'N' or 'V'.
*
*  RCONDV  (output) REAL
*          If SENSE = 'V' or 'B', RCONDV contains the reciprocal
*          condition number for the selected right invariant subspace.
*          Not referenced if SENSE = 'N' or 'E'.
*
*  WORK    (workspace/output) COMPLEX array, dimension (LWORK)
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
*  LWORK   (input) INTEGER
*          The dimension of the array WORK.  LWORK >= max(1,2*N).
*          Also, if SENSE = 'E' or 'V' or 'B', LWORK >= 2*SDIM*(N-SDIM),
*          where SDIM is the number of selected eigenvalues computed by
*          this routine.  Note that 2*SDIM*(N-SDIM) <= N*N/2.
*          For good performance, LWORK must generally be larger.
*
*  RWORK   (workspace) REAL array, dimension (N)
*
*  BWORK   (workspace) LOGICAL array, dimension (N)
*          Not referenced if SORT = 'N'.
*
*  INFO    (output) INTEGER
*          = 0: successful exit
*          < 0: if INFO = -i, the i-th argument had an illegal value.
*          > 0: if INFO = i, and i is
*             <= N: the QR algorithm failed to compute all the
*                   eigenvalues; elements 1:ILO-1 and i+1:N of W
*                   contain those eigenvalues which have converged; if
*                   JOBVS = 'V', VS contains the transformation which
*                   reduces A to its partially converged Schur form.
*             = N+1: the eigenvalues could not be reordered because some
*                   eigenvalues were too close to separate (the problem
*                   is very ill-conditioned);
*             = N+2: after reordering, roundoff changed values of some
*                   complex eigenvalues so that leading eigenvalues in
*                   the Schur form no longer satisfy SELECT=.TRUE.  This
*                   could also be caused by underflow due to scaling.
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            SCALEA, WANTSB, WANTSE, WANTSN, WANTST,
     $                   WANTSV, WANTVS
      INTEGER            HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
     $                   ITAU, IWRK, K, MAXB, MAXWRK, MINWRK
      REAL               ANRM, BIGNUM, CSCALE, EPS, SMLNUM
*     ..
*     .. Local Arrays ..
      REAL               DUM( 1 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           CCOPY, CGEBAK, CGEBAL, CGEHRD, CHSEQR, CLACPY,
     $                   CLASCL, CTRSEN, CUNGHR, SLABAD, SLASCL, XERBLA
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      REAL               CLANGE, SLAMCH
      EXTERNAL           LSAME, ILAENV, CLANGE, SLAMCH
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN, SQRT
*     ..
*     .. Executable Statements ..
*
*     Test the input arguments
*
      INFO = 0
      WANTVS = LSAME( JOBVS, 'V' )
      WANTST = LSAME( SORT, 'S' )
      WANTSN = LSAME( SENSE, 'N' )
      WANTSE = LSAME( SENSE, 'E' )
      WANTSV = LSAME( SENSE, 'V' )
      WANTSB = LSAME( SENSE, 'B' )
      IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
         INFO = -1
      ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
         INFO = -2
      ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
     $         ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
         INFO = -4
      ELSE IF( N.LT.0 ) THEN
         INFO = -5
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
         INFO = -11
      END IF
*
*     Compute workspace
*      (Note: Comments in the code beginning "Workspace:" describe the
*       minimal amount of real workspace needed at that point in the
*       code, as well as the preferred amount for good performance.
*       CWorkspace refers to complex workspace, and RWorkspace to real
*       workspace. NB refers to the optimal block size for the
*       immediately following subroutine, as returned by ILAENV.
*       HSWORK refers to the workspace preferred by CHSEQR, as
*       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
*       the worst case.
*       If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
*       depends on SDIM, which is computed by the routine CTRSEN later
*       in the code.)
*
      MINWRK = 1
      IF( INFO.EQ.0 .AND. ( LWORK.GE.1 ) ) THEN
         MAXWRK = N + N*ILAENV( 1, 'CGEHRD', ' ', N, 1, N, 0 )
         MINWRK = MAX( 1, 2*N )
         IF( .NOT.WANTVS ) THEN
            MAXB = MAX( ILAENV( 8, 'CHSEQR', 'SN', N, 1, N, -1 ), 2 )
            K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'CHSEQR', 'SN', N, 1,
     $          N, -1 ) ) )
            HSWORK = MAX( K*( K+2 ), 2*N )
            MAXWRK = MAX( MAXWRK, HSWORK, 1 )
         ELSE
            MAXWRK = MAX( MAXWRK, N+( N-1 )*
     $               ILAENV( 1, 'CUNGHR', ' ', N, 1, N, -1 ) )
            MAXB = MAX( ILAENV( 8, 'CHSEQR', 'SV', N, 1, N, -1 ), 2 )
            K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'CHSEQR', 'SV', N, 1,
     $          N, -1 ) ) )
            HSWORK = MAX( K*( K+2 ), 2*N )
            MAXWRK = MAX( MAXWRK, HSWORK, 1 )
         END IF
         WORK( 1 ) = MAXWRK
      END IF
      IF( LWORK.LT.MINWRK ) THEN
         INFO = -15
      END IF
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'CGEESX', -INFO )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( N.EQ.0 ) THEN
         SDIM = 0
         RETURN
      END IF
*
*     Get machine constants
*
      EPS = SLAMCH( 'P' )
      SMLNUM = SLAMCH( 'S' )
      BIGNUM = ONE / SMLNUM
      CALL SLABAD( SMLNUM, BIGNUM )
      SMLNUM = SQRT( SMLNUM ) / EPS
      BIGNUM = ONE / SMLNUM
*
*     Scale A if max element outside range [SMLNUM,BIGNUM]
*
      ANRM = CLANGE( 'M', N, N, A, LDA, DUM )
      SCALEA = .FALSE.
      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
         SCALEA = .TRUE.
         CSCALE = SMLNUM
      ELSE IF( ANRM.GT.BIGNUM ) THEN
         SCALEA = .TRUE.
         CSCALE = BIGNUM
      END IF
      IF( SCALEA )
     $   CALL CLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
*
*
*     Permute the matrix to make it more nearly triangular
*     (CWorkspace: none)
*     (RWorkspace: need N)
*
      IBAL = 1
      CALL CGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
*
*     Reduce to upper Hessenberg form
*     (CWorkspace: need 2*N, prefer N+N*NB)
*     (RWorkspace: none)
*
      ITAU = 1
      IWRK = N + ITAU
      CALL CGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
     $             LWORK-IWRK+1, IERR )
*
      IF( WANTVS ) THEN
*
*        Copy Householder vectors to VS
*
         CALL CLACPY( 'L', N, N, A, LDA, VS, LDVS )
*
*        Generate unitary matrix in VS
*        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
*        (RWorkspace: none)
*
         CALL CUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
     $                LWORK-IWRK+1, IERR )
      END IF
*
      SDIM = 0
*
*     Perform QR iteration, accumulating Schur vectors in VS if desired
*     (CWorkspace: need 1, prefer HSWORK (see comments) )
*     (RWorkspace: none)
*
      IWRK = ITAU
      CALL CHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS,
     $             WORK( IWRK ), LWORK-IWRK+1, IEVAL )
      IF( IEVAL.GT.0 )
     $   INFO = IEVAL
*
*     Sort eigenvalues if desired
*
      IF( WANTST .AND. INFO.EQ.0 ) THEN
         IF( SCALEA )
     $      CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR )
         DO 10 I = 1, N
            BWORK( I ) = SELECT( W( I ) )
   10    CONTINUE
*
*        Reorder eigenvalues, transform Schur vectors, and compute
*        reciprocal condition numbers
*        (CWorkspace: if SENSE is not 'N', need 2*SDIM*(N-SDIM)
*                     otherwise, need none )
*        (RWorkspace: none)
*
         CALL CTRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM,
     $                RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1,
     $                ICOND )
         IF( .NOT.WANTSN )
     $      MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
         IF( ICOND.EQ.-14 ) THEN
*
*           Not enough complex workspace
*
            INFO = -15
         END IF
      END IF
*
      IF( WANTVS ) THEN
*
*        Undo balancing
*        (CWorkspace: none)
*        (RWorkspace: need N)
*
         CALL CGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS,
     $                IERR )
      END IF
*
      IF( SCALEA ) THEN
*
*        Undo scaling for the Schur form of A
*
         CALL CLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
         CALL CCOPY( N, A, LDA+1, W, 1 )
         IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN
            DUM( 1 ) = RCONDV
            CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
            RCONDV = DUM( 1 )
         END IF
      END IF
*
      WORK( 1 ) = MAXWRK
      RETURN
*
*     End of CGEESX
*
      END

ext: f

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


      SUBROUTINE DCSTAT(Cerr,Citem,Istat)
!=======================================================================
!     LATEST CORRECTION BY
!
!     PURPOSE
!           Decode channels for statistical analysis
!
!     PROGRAMMED BY:  bjorn.melhus@akersolutions.com
!     CREATED......:  10.05.2021
!=======================================================================
      use clist
      IMPLICIT NONE
      INCLUDE 'redwri'
      CHARACTER*(*) Cerr , Citem
      CHARACTER*14 channel
      LOGICAL MATCH
      INTEGER iend , IFERR , imemb , ilen , intnum ,Istat
!
      intnum=0
 20   CALL MOVEIP(iend)
      IF ( iend.EQ.0 ) THEN
         CALL DCSTRG(channel,ilen,3)
!
         IF ( IFERR().LT.0 ) THEN
            WRITE (IOUT,Citem) channel
            WRITE (IOUT,Cerr) 'ERROR IN READING CHANNEL ' , &
                              'PROGRAM STOPPED'
            Istat = -1
            GOTO 100
         ENDIF
!
         intnum=intnum+1
         LIST%CHAN_REP(intnum) = channel
         LIST%NCHAN_REP=LIST%NCHAN_REP+1
!
         GOTO 20
      ENDIF
 100  RETURN
      END


text:
      program prog
      common /info/ nums(2)
      equivalence (nums(1), n1), (n2, nums(2))
      n1 = 12
      n2 = 34
      write(6, *) nums
      end