      Subroutine rhoofr (r,rhor,nff)
c *** see RH 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)) GoTo 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 GoTo 30
c     Gaussian densities for both n and p matter
      rhor(1) = 0.
      rhl = (r/rhl1p)**2
      If (rhl.gt.150.) GoTo 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.) GoTo 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
