PROGRAM FORGRID

!PURPOSE: Compute & write a Fourier map over the entire unit cell 

	INCLUDE		'../INCLDS/COPYRIGT.FOR'

!PSEUDOCODE:

!CALLING ARGUMENTS:

!INCLUDE STATEMENTS:

	INCLUDE 	'../INCLDS/HEADSCOM.FOR'

	INCLUDE		'../INCLDS/SPGCOMI.FOR'

	INCLUDE		'../INCLDS/FMPDATCM.FOR'

!LOCAL VARIABLES:

	INTEGER*4	IUOUT		!Output file
	CHARACTER*68	EXPNAM		! Experiment name
	CHARACTER*66	TITLE(2)	!Titles
	CHARACTER*12	KEYVAL		!ISAM key
	REAL*4		ABC(3),ANGLES(3)!Lattice parameters
	INTEGER*4	IPHAS		!Phase no.
	REAL*4		TPOS(3)		!Rho test point
	REAL*4		FMULT		!Rho multiplier
	CHARACTER*68	TEXT		!ISAM file read/write buffer
C	INTEGER*4	NPHAS(9)	!Phase flags
      CHARACTER*4   MAPTYP              !Map type
      INTEGER*4     NOXYZT(3)           !No. of steps in map
      INTEGER*4     IRHO                !Array pointer
      INTEGER*4     IDUM(3)             !Dummy array

	integer jx,jy,jz
	character*50 filout
!SUBROUTINES CALLED: 
	
!FUNCTION DEFINITIONS: 

	INTEGER*4	GETRHO		!=0 if rho obtained, =1 if error
	INTEGER*4	READEXP		!ISAM file read function

!DATA STATEMENTS:

!CODE:
 
	CALL STRTRN('FORGRID','SHARED','TERM',IUEXP,IULST,IUTRM)
	CALL PROGNAM(IULST,'FORGRID',
	1	'Fourier map to grid convertion program')
	INQUIRE (IUEXP,NAME=EXPNAM)
	NEND = INDEX(EXPNAM,'.EXP')-1

C get the phases in use
C	ISAM = READEXP(IUEXP,' EXPR NPHAS ',TEXT)
C	READ(TEXT,'(9I5)') NPHAS
C
C get the phase number used for map computations -- currently there can only 
C be one phase used for all current map computations. However, there can
C be a problem if old map files are laying around computed with other phases 
C -- this problem can happen within all GSAS Fourier map programs as well.
	TEXT = ' ' 
	ISAM = READEXP(IUEXP,'  FOUR CDAT1',TEXT)
	READ(TEXT,'(11X,I2)') IPHAS
	IF (IPHAS .EQ. 0) THEN
	   WRITE (*,*) 'No Fourier Phase was found. '//
	1	'Has a Map been computed?'
	   stop 'no map'
	ENDIF
	WRITE(KEYVAL,'(A,I1,A)') 'CRS',IPHAS,' '
	ISAM = READEXP(IUEXP,KEYVAL(1:6)//'ABC   ',TEXT)
	READ (TEXT,'(3F10.0)') ABC
	ISAM = READEXP(IUEXP,KEYVAL(1:6)//'ANGLES',TEXT)
	READ (TEXT,'(3F10.0)') ANGLES

C read in the Fourier map
	CALL OPNFMP(IUEXP,'OLD',MAPTYP,IUFMP)
	READ (IUFMP,REC=2) NXYZI,IDUM,IDUM,MSECT,NRHO,
     1      RMAX,RMIN,SUMRHO
	CALL GETVM(IRHO,4*NRHO)
	CALL REDFMP(IUEXP,IUFMP,%val(IRHO),IPHAS,NOXYZT,RMIN,RMAX,
     1      FMULT,TITLE,MAPTYP)
	CALL FREEVM(IRHO)
        
	CALL GETUNIT(IUOUT)
	open(unit=IUOUT,file=EXPNAM(:NEND)//'.grd',status='unknown')
	inquire(unit=IUOUT,name=filout)
	write (*,*) 'Writing file: ',filout

	write(IUOUT,'(A)') Title(1)
	write(IUOUT,'(3f10.5,3f10.4)') ABC,ANGLES
	write (IUOUT,'(3i4)') NXYZI

	do jx = 1,NXYZI(1)
	   tpos(1) = (jx-1.)/NXYZI(1)
	   do jy = 1,NXYZI(2)
	      tpos(2) = (jy-1.)/NXYZI(2)
	      do jz = 1,NXYZI(3)
		 tpos(3) = (jz-1.)/NXYZI(3)
		 IERR = GETRHO(TPOS,RHO)
		 write (IUOUT,*) RHO
!		 write (IUOUT,*) TPOS,RHO
	      ENDDO
	   ENDDO
	ENDDO
	close(unit=IUOUT)
	END