************************************************************************
      SUBROUTINE DENS(PATH)
************************************************************************
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)

      PARAMETER (NSPEC=4)

      real*8 path(8)
      COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY
      COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY
      COMMON /QPARAM/RREF,RBTM,RREFP,
     &               PF0(NSPEC),GF0(NSPEC),
     &               XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP,
     &               RRLI,HWLI,
     &               XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2,
     &               RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT,
     &               NPATH,NLPMAX,INDOUT,INDMSH,IEOF,
     &               IYEAR,IMNTH,IDATE
      COMMON /QDENS /XNS   (NSPEC),
     &               DNSDRR(NSPEC),
     &               DNSDTH(NSPEC),
     &               DNSDPH(NSPEC)
      COMMON /QDNSDE/XNDE  (NSPEC),
     &               DNDDRR(NSPEC)
      COMMON /QDNSLI/XNLI  ,DNLDRR
      COMMON /QDNSAO/XNAO  ,DNADRR,DNADTH,DNADPH
      COMMON /QDNSOP/XNOP  ,DNODRR,DNODTH,DNODPH

      call FLPRM(PATH)
      CALL DENSDE(PATH)
      CALL DENSLI(PATH)
      CALL DENSAO(PATH)
      CALL DENSOP(PATH)

      DO 10 IS=1,NSPEC
     
      XNS(IS) = XNDE(IS)*XNLI*XNAO*XNOP

*...  1/N * dn/dr
      DNSDRR(IS) = DNDDRR(IS)+DNLDRR+DNADRR+DNODRR
      DNSDTH(IS) =                   DNADTH+DNODTH
      DNSDPH(IS) = DNADPH+DNODPH     
              
   10 CONTINUE
      RETURN
      END

************************************************************************
      SUBROUTINE DENSAO(PATH)
************************************************************************
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)

      PARAMETER (NSPEC=4)
      real*8 path(8)
      COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY
      COMMON /QPARAM/RREF,RBTM,RREFP,
     &               PF0(NSPEC),GF0(NSPEC),
     &               XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP,
     &               RRLI,HWLI,
     &               XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2,
     &               RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT,
     &               NPATH,NLPMAX,INDOUT,INDMSH,IEOF,
     &               IYEAR,IMNTH,IDATE
      COMMON /QIGRF /YR    ,YT    ,YP    ,
     &               DYRDRR,DYTDRR,DYPDRR,
     &               DYRDTH,DYTDTH,DYPDTH,
     &               DYRDPH,DYTDPH,DYPDPH,
     &               BB    ,
     &               DBBDRR,DBBDTH,DBBDPH
      COMMON /QDNSAO/XNAO  ,DNADRR,DNADTH,DNADPH
      COMMON /QFLPRM/TH0   ,
     &               DT0DRR,DT0DTH,DT0DPH,
     &               PH0   ,
     &               DP0DRR,DP0DTH,DP0DPH,
     &               BB0   ,
     &               DB0DRR,DB0DTH,DB0DPH,
     &               XLV   ,
     &               DLVDRR,DLVDTH,DLVDPH

      XNAO   = 1.0D0
      DNADRR = 0.0D0
      DNADTH = 0.0D0
      DNADPH = 0.0D0

      RETURN
      END

************************************************************************
      SUBROUTINE DENSDE(PATH)
************************************************************************
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)

      PARAMETER (NSPEC=4)
      real*8 path(8)
      COMMON /QPARAM/RREF,RBTM,RREFP,
     &               PF0(NSPEC),GF0(NSPEC),
     &               XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP,
     &               RRLI,HWLI,
     &               XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2,
     &               RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT,
     &               NPATH,NLPMAX,INDOUT,INDMSH,IEOF,
     &               IYEAR,IMNTH,IDATE
      COMMON /QDNSDE/XNDE  (NSPEC),
     &               DNDDRR(NSPEC)

      GPH    = RREF*(PATH(1)-RREF)/PATH(1)
      DGPHDR = (RREF/PATH(1))**2
      ETEXP2 = ETA(2)*EXP(-GPH*SHINV(2))
      ETEXP3 = ETA(3)*EXP(-GPH*SHINV(3))
      ETEXP4 = ETA(4)*EXP(-GPH*SHINV(4))

      SUM1   = ETEXP2
     &        + ETEXP3
     &        + ETEXP4
      SUM2   = ETEXP2*SHINV(2)
     &        + ETEXP3*SHINV(3)
     &        + ETEXP4*SHINV(4)

      XNDE  (1) = dSQRT(SUM1)
      XNDE  (2) = ETEXP2/XNDE(1)
      XNDE  (3) = ETEXP3/XNDE(1)
      XNDE  (4) = ETEXP4/XNDE(1)

      DNDDRR(1) =-DGPHDR*SUM2/SUM1/2.0D0
      DNDDRR(2) =-DGPHDR*SHINV(2)-DNDDRR(1)
      DNDDRR(3) =-DGPHDR*SHINV(3)-DNDDRR(1)
      DNDDRR(4) =-DGPHDR*SHINV(4)-DNDDRR(1)

      RETURN
      END

************************************************************************
      SUBROUTINE DENSLI(PATH)
************************************************************************
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)

      PARAMETER (NSPEC=4)
      real*8 path(8)
      COMMON /QPARAM/RREF,RBTM,RREFP,
     &               PF0(NSPEC),GF0(NSPEC),
     &               XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP,
     &               RRLI,HWLI,
     &               XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2,
     &               RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT,
     &               NPATH,NLPMAX,INDOUT,INDMSH,IEOF,
     &               IYEAR,IMNTH,IDATE
      COMMON /QDNSLI/XNLI, DNLDRR
                 
      IF ((PATH(1)-RRLI).GE.0) THEN
      ALTNRM = (PATH(1)-RRLI)/HWLI
      XNLI   = 1.0D0-dEXP(-ALTNRM**2)
      DNLDRR = dEXP(-ALTNRM**2)*ALTNRM*2.0D0/HWLI/XNLI 
      ELSE
      XNLI=0.0
      DNLDRR=0.0  
      END IF    

      RETURN
      END

************************************************************************
      SUBROUTINE DENSOP(PATH)
************************************************************************
      IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N)

      PARAMETER (NSPEC=4)
      real*8 path(8)
      COMMON /QCONST/PAI,RAD,DEG,RE,CVLCTY
      COMMON /QPATH /FKC,DEL,EPS,GMU,ALF,RAY
      COMMON /QPARAM/RREF,RBTM,RREFP,
     &               PF0(NSPEC),GF0(NSPEC),
     &               XN0,THERM,ETA(NSPEC),SHINV(NSPEC),SHINVP,
     &               RRLI,HWLI,
     &               XLPP,HWPP,ERPP,DRR1,DTH1,DPH1,DRR2,DTH2,DPH2,
     &               RRMIN,RRMAX,ERRR,ERRA,HMIN,DISOUT,
     &               NPATH,NLPMAX,INDOUT,INDMSH,IEOF,
     &               IYEAR,IMNTH,IDATE
      COMMON /QDNSOP/XNOP  ,DNODRR,DNODTH,DNODPH

      XNOP = XN0
      dNODRR = 0.0D0
      DNODTH = 0.0d0
      dnodph = 0.0d0

      RETURN
      END