      SUBROUTINE RHOOFR (R,RHOR,NFF)
C *** SEE R H LANDAU, PROGRAM LPOTT, 1981
C     RHOR(R) CALULATES UNNORMALIZED NUCLEAR DENSITY IN COORD SPACE
C     IF A LE 16 ,SHELL MODEL(GAUSSIAN) FORM IS USED
C     IF A GT 16, A GENERALIZED FERMI SHAPE IS USED
C          WITH ?C? =ACH, AND ?Z? = ACM
C     N.B. R IS IN MEV-1 AND W IS DIMENSIONLESS
C
      IMPLICIT  REAL*8   (A-H,O-Z)
      REAL*8 IMTIJ,MPI,MN
      DIMENSION NIFTY(20), BNUC(20), BNUCF(20), BN(16), BNF(16), RETIJ
     1(14), IMTIJ(14), RHOR(4)
      COMMON /SEC2/ BNUC,BN,BNUCF,BNF,RETIJ,IMTIJ,HBARC,PI,MPI,MN,NZ,NES
     1,NWAVES,NIFTY,NA
      COMMON /SEC4/ ACHP,ACMP,WSP,ACHN,ACMN,WSN
C
      DATA NDATA/0/,ACHECK/0./
      ANEW = ACHP+ACMP+ACHN+ACMN+WSP+WSN
      IF ((NDATA.NE.0).AND.(ANEW.EQ.ACHECK)) GO TO 10
C     INTIALIZE
      ACHECK = ANEW
      NDATA = 1
      RHL1P = ACHP/HBARC
      RHL1N = ACHN/HBARC
      RHL2P = ACMP/HBARC
      RHL2N = ACMN/HBARC
      AP = (NZ-2.)/6.
      AN = (NA-NZ-2.)/6.
      IF (AP.LT.0.) AP = 0.
      IF (AN.LT.0.) AN = 0.
 10   IF ((NA.GT.16).OR.(ABS(WSP).GT.1.E-4).OR.(ABS(WSN).GT.1.E-4))
     1 GO TO 30
C     GAUSSIAN DENSITIES FOR BOTH N AND P MATTER
      RHOR(1) = 0.
      RHL = (R/RHL1P)**2
      IF (RHL.GT.150.) GO TO 20
      RHOR(1) = (1.+AP*(R/RHL2P)**2)*EXP(-RHL)
 20   IF (NFF.EQ.1) RETURN
      RHOR(2) = 0.
      RHL = (R/RHL1N)**2
      IF (RHL.GT.150.) RETURN
      RHOR(2) = (1.+AN*(R/RHL2N)**2)*EXP(-RHL)
      RETURN
C     GENERALIZED WOOD SAXON
 30   RHOR(1) = 0.
      RHL = (R-RHL1P)/RHL2P
      IF (RHL.GT.150.) GO TO 40
      RHOR(1) = (1.+WSP*(R/RHL1P)**2)/(1.+EXP(RHL))
 40   IF (NFF.EQ.1) RETURN
      RHOR(2) = 0.
      RHOR(3) = 0.
      RHOR(4) = 0.
      RHL = (R-RHL1N)/RHL2N
      IF (RHL.GT.150.) RETURN
      RHOR(2) = (1.+WSN*(R/RHL1N)**2)/(1.+EXP(RHL))
      IF (NFF.EQ.2) RETURN
C     CALC OF FOUR TF OF RHO**2
      RHOR(3) = RHOR(1)**2
      RHOR(4) = RHOR(2)**2
      RETURN
      END
