C- READ CCP4 density map, calculate min, max, mean, and dev of density
C- Copy map to new version, including density values.                

	CHARACTER*1 BLINE(1024) 
	CHARACTER*36 RECLINE(30)
	EQUIVALENCE (RECLINE(1),BLINE(1))
	real*4 rline(256),TCELL(6),CELL(6)
	EQUIVALENCE (RLINE(1),BLINE(1)),(TCELL(1),RLINE(11))
	INTEGER*4 LINT(256),n,ii
	EQUIVALENCE (LINT(1),BLINE(1))
	INTEGER*2 NX,NY,NZ,NXYZ(3)
	EQUIVALENCE (NX,NXYZ(1)),(NY,NXYZ(2)),(NZ,NXYZ(3))
	CHARACTER*80 TTITLE,TITLE,TSYMTRY,SYMTRY
	EQUIVALENCE (TTITLE,RLINE(57))
	EQUIVALENCE (TSYMTRY,RLINE(1))
	CHARACTER*50 INFILE,OUTFILE,FILENAME
	real RMAP
        CHARACTER*40 FNAME,ASTRING
	CHARACTER*1 OUTSTRING(1200),HSTRING(16),xyzstr(3)
	LOGICAL*4 EFLG
	integer*2 ZLP
	CHARACTER DAT*12,TIM*12
	REAL*8 SUM,SSQU
       	pi=4*ATAN(1.)
	R64K=256*256
 	R32K=128*256
	DATA XYZSTR/'X','Y','Z'/

	write(6,*)'eflg=',eflg
c	DO 2 I=1,16
c2	WRITE(6,*)I,'   ',HSTRING(I)

	WRITE(6,1000)
1000	FORMAT(//' CCP4 DENSITY MAP density fix.')                       
5	WRITE(6,1100)
1100	FORMAT(/'$Input filename (CR to quit- do not give version!):  ')
	READ(5,1200) INFILE
1200	FORMAT(A)
	IF (INFILE(:5).EQ.'     ') STOP
	  INQUIRE(FILE=INFILE, NAME=FILENAME,EXIST=EFLG)
	WRITE(6,*)'INFILE=',INFILE,'EFLG=',EFLG,'FILENAME=',FILENAME
	  IF (.NOT.EFLG) GOTO 5
        OPEN (UNIT=2,FILE=INFILE,STATUS='OLD',FORM='UNFORMATTED',
     .   READONLY)
C1400	write(6,*)'Enter window (pixels; xmin,xmax,ymin,ymax;',
C     .             ' 0 0 0 0 for default):'
C	read(5,*) i1,i2,j1,j2
C	WRITE(6,*)'Enter density limits (min,max) or 0,0 to autoscale'
C	WRITE(6,*)'Density autoscaling requires two passes through the file.'
C	READ(5,*)X0,XM
C	  WRITE(6,1600)
	
	READ(2)(RLINE(I),I=1,256)
C  These dimensions are actual dimens present in map- maybe < full unit cell
	NFAST=LINT(1)
	NMED=LINT(2)
	NSLOW=LINT(3)
	maptype=LINT(4)
	nsfast=lint(5)
	nsmed=lint(6)
	nsslow=lint(7)
   	DO 290 I=1,6
290	CELL(I)=TCELL(I)
C   These dimensions are sampling dimensions, corresp to full unit cell:
c	correspond to true x,y,z axes, regardles of order in map.
        NX=LINT(8)
	NY=LINT(9)
	NZ=LINT(10)
	iaxf=lint(17)
	iaxm=lint(18)
	iaxsl=lint(19)
C  MFAST, MED, SLOW ARE FULL UNIT CELL IN THOSE DIRECTIONS
	mfast=nxyz(iaxf)
	mmed=nxyz(iaxm)
	mslow=nxyz(iaxsl)
	olddmin=rline(20)
	olddmax=rline(21)
	olddmean=rline(22)
	nsg=lint(23)
	nsymchar=lint(24)
	oldddev=RLINE(55)
	ntitle=lint(56)
	TITLE=TTITLE
	WRITE(6,*)'NFAST, MEDIUM, SLOW=',NFAST,NMED,NSLOW
	WRITE(6,*)'MAP TYPE=',MAPTYPE
	WRITE(6,*)'START (FAST, MED, SLOW)=',NSFAST,NSMED,NSSLOW
	WRITE(6,*)'CELL: ',(CELL(I),I=1,6)
	WRITE(6,*)'SAMPLING(NX, NY, NZ)=',NX, NY, NZ
	WRITE(6,*)'AXIS PERMUTAION (FAST, MEDIUM, SLOW)=',IAXF,
     & IAXM,IAXSL
        WRITE(6,*)'old DMIN, DMAX, DMEAN, DDEV=',OLDDMIN,OLDDMAX,OLDDMEAN,OLDDDEV
	WRITE(6,*)'SPACE GROUP #',NSG
	WRITE(6,*)'# CHAR OF SYM INFO:',NSYMCHAR
	WRITE(6,*)'NUMBER OF TITLES:',NTITLE
	WRITE(6,*) TITLE

	READ(2)(RLINE(I),I=1,256)
	SYMTRY=TSYMTRY
	WRITE(6,*)SYMTRY
	IF((NFAST.GT.500).OR.(NMED.GT.500)) STOP 
     &	'DIMENSIONED FOR 500X500 SECTIONS ONLY!'
	
        II=1+NSYMCHAR/4
C  	!!!IF NSYMCHAR NOT A MULTIPLE OF 4, SHOULD ROUND UP, NOT DOWN?	
	SUM=0
	SSQU=0
	DMIN=1E31
	DMAX=-1E31
	NREC=2
	DO 300 LAYER=1,nslow
	DO 300 J=1,NMED
	DO 300 K=1,NFAST
	IF (II.GT.256) THEN
	  READ(2,END=900)(RLINE(IJ),IJ=1,256)
	  NREC=NREC+1
	  II=1
	  ENDIF

	RMAP     =RLINE(II)
	II=II+1
	if (rmap.lt.dmin) dmin=rmap
	if (rmap.gt.dmax) dmax=rmap
	sum=sum+rmap
	ssqu=ssqu+rmap*rmap
300	CONTINUE

380	CLOSE(UNIT=2)

	NN=(NFAST*NMED*NSLOW)  
	DMEAN=SUM/NN
	DDEV=SQRT(SSQU/NN -DMEAN*DMEAN)
	WRITE(6,*)' PARAM       OLD VALUE        ACTUAL VALUE'
	WRITE(6,*)'DMIN  ',OLDDMIN,DMIN
	WRITE(6,*)'DMAX  ',OLDDMAX,DMAX
	WRITE(6,*)'DMEAN  ',OLDDMEAN,DMEAN
	WRITE(6,*)'DDEV  ',OLDDDEV,DDEV

C	rewind(unit=2)
	CLOSE(UNIT=2)
        OPEN (UNIT=2,FILE=INFILE,STATUS='OLD',FORM='UNFORMATTED',
     .   READONLY)
        OPEN (UNIT=3,FILE=infile,STATUS='NEW',FORM='UNFORMATTED',
     & RECL=256, RECORDTYPE='FIXED', CARRIAGECONTROL='NONE')

	READ(2)(RLINE(I),I=1,256)
	rline(20)=DMIN
	rline(21)=DMAX
	rline(22)=DMEAN
	RLINE(55)=DDEV
        WRITE(3)(RLINE(IJ),IJ=1,256)
	DO 390 K=2,NREC
	READ(2)(RLINE(I),I=1,256)
390     WRITE(3)(RLINE(IJ),IJ=1,256)

	STOP 'NORMAL END'
900	STOP 'END OF FILE BEFORE END OF MAP'
	END