C ---------------------------------------------------------------------
C rhosmp.for
C ---------------------------------------------------------------------
! 27 Jun 2016 - Changed calling sequence from (gravity) to (gravity,coriolis)
! 30 Apr 2010 - Correction
C 25 Mar 2010 - Changed CALL PDENSITY to CALL PDENSITYg
C 25 Mar 2010 - Changed CALL GRAV to CALL GRAVg
C 22 Mar 2010 - Changed subroutine name to RHOSMP
C 22 Mar 2010 - Added ENTRY DENSITY
C 22 Mar 2010 - Added ENTRY DENSITYg
C 22 Mar 2010 - Added calculation of parameters in /DS/
C  6 Aug 2009 - Altered a comment about the units for density.
C  4 May 2009 - Replaced REAL*8 with DOUBLE PRECISION
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 13 Mar 2008 - Corrected length in common block /MM/
C 23 Mar 2007 - Added density, temperature, pressure, and molecular
C               weight common blocks   
C 19 Jul 2006 - Separated idensity (no argument) and idensityg (w. argument)
C  2 Feb 2006 - Deleted common B15 in idensity
C 26 Jan 2006
C 25 Jan 2006 
C 24 Jan 2006 - Different entry points for acoustics and gravity
C 23 Jan 2006
C ---------------------------------------------------------------------

C     HYDRODYNAMIC STABILITY DENSITY MODEL
C     The value for density is in kilograms per cubic centimeter when 
C     pressure is in Newtons per square meter. However, the ratio of density 
C     to its gradients will be correct, independent of the value of pressure,
C     as long as the ratio of pressure to its gradients is correct, 
C     and will give the correct Brunt-Vaisala frequency.             6 Aug 2009

C ---------------------------------------------------------------------
      SUBROUTINE RHOSMP                                           ! 22 Mar 2010
C ---------------------------------------------------------------------

      IMPLICIT double precision (a-h,o-z)

      TYPE :: values_and_derivs_1_2
         REAL(KIND=8) :: value
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: ptmtm,prtm,pthtm,pphtm
         REAL(KIND=8) :: ptmr,prr,pthr,pphr
         REAL(KIND=8) :: ptmth,prth,pthth,pphth
         REAL(KIND=8) :: ptmph,prph,pthph,pphph
      END TYPE 
      TYPE (values_and_derivs_1_2) :: rho
      TYPE (values_and_derivs_1_2) :: pressure
      TYPE (values_and_derivs_1_2) :: temperature
      TYPE (values_and_derivs_1_2) :: molecular_weight
      TYPE (values_and_derivs_1_2) :: potentialrho

      TYPE :: gravity_and_derivs_1
         REAL(KIND=8) :: gs,gsr,gsth,gsph
         REAL(KIND=8) :: ptm,pr,pth,pph
         REAL(KIND=8) :: prtm,prr,prth,prph
         REAL(KIND=8) :: pthtm,pthr,pthth,pthph
         REAL(KIND=8) :: pphtm,pphr,pphth,pphph
      END TYPE 
      TYPE (gravity_and_derivs_1) :: gravity
      TYPE (gravity_and_derivs_1) :: coriolis                      ! 27 Jun 2016
      REAL(KIND=8) :: T,g,gr,gth,gph,mw,p,hscale,rhov

      PARAMETER (NWARSZ=1000)

C     COMMON DECK "WW" INSERTED HERE
      COMMON/WW/dum(10),MAXW,W(NWARSZ)

C     DENSITY     650-674
      EQUIVALENCE (W(650),dmodel),(W(651),dform),(W(652),did)

      REAL(KIND=8) KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)

C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3                      ! 13 Mar 2008
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10

C     COMMON DECK "CC" INSERTED HERE
      REAL(KIND=8) modc
      COMMON/CC/modc(4),cs,pcst,pcsr,pcsth,pcsph

C     COMMON DECK "DS" INSERTED HERE                            ! 27 Mar 2007
      DOUBLE PRECISION ds,modds
      COMMON/DS/MODDS(4),ds,pdst,pdsr,pdsth,pdsph

C     COMMON DECK "TT" INSERTED HERE                            ! 27 Mar 2007
      DOUBLE PRECISION MODT
      COMMON/TT/MODT(4),T,PTT,PTR,PTTH,PTPH

C     COMMON DECK "MM" INSERTED HERE                            ! 27 Mar 2007
      DOUBLE PRECISION m,modm
      COMMON/MM/modm(4),m,pmt,pmr,pmth,pmph
     1 ,XO2,XN2,XO,XN,XH2O,XO3                     ! 5 Mar 2007 ! 13 Mar 2008

C     COMMON DECK "PP" INSERTED HERE                            ! 27 Mar 2007
      DOUBLE PRECISION modp
      COMMON/PP/modp(4),p,ppt,ppr,ppth,ppph
      EQUIVALENCE(w(1),earthr)

C ---------------------------------------------------------------------

      ENTRY DENSITYg(rho,potentialrho)
      ENTRY DENSITY(rho,potentialrho)

c      print *, ' before tempg'
      CALL tempg(temperature) ! 24 January 2006
c      print *, ' after tempg'
      T = temperature%value

c      print *, ' before grav'
      CALL gravg(gravity,coriolis)                   ! 25 Mar 2010 ! 27 Jun 2016
c      print *, ' after grav'
      g = gravity%gs
      gr = gravity%gsr
      gth = gravity%gsth
      gph = gravity%gsph

      CALL molwtg(molecular_weight) ! 24 January 2006
      mw = molecular_weight%value

      CALL presg(pressure) ! 24 January 2006
      p = pressure%value

      hscale = -p/pressure%pr

      rho%value = (p*mw)/(rgas*T)
      ds = rho%value                                              ! 22 Mar 2010
      rhov = rho%value

      rho%ptm = rhov * ((pressure%ptm/p) - (temperature%ptm/T) +
     & (molecular_weight%ptm/mw))
      pdst = rho%ptm                                              ! 22 Mar 2010
      rho%pr = rhov * ((pressure%pr/p) - (temperature%pr/T) +
     & (molecular_weight%pr/mw))
      pdsr = rho%pr                                               ! 22 Mar 2010
      rho%pth = rhov * ((pressure%pth/p) - (temperature%pth/T) +
     & (molecular_weight%pth/mw))
      pdsth = rho%pth                                             ! 22 Mar 2010
      rho%pph = rhov * ((pressure%pph/p) - (temperature%pph/T) +
     & (molecular_weight%pph/mw))
      pdsph = rho%pph                                             ! 22 Mar 2010

c      AA = pressure%pr/p
c      BB = temperature%pr/T
c      print*,'Pressure Term = ',AA
c      print*,'Temperature Term = ',BB
c      print*,'pressure = ',p
c      print*,'pressure gradient = ',pressure%pr
c      print*,'temperature = ',temperature%value
c      print*,'temperature gradient = ',temperature%pr
c      print*,'rho%pr (rho) = ',rho%pr
c      PRINT *,'temperature%ptm = ',temperature%ptm
c      PRINT *,'temperature%pr = ',temperature%pr
c      PRINT *,'temperature%ptm = ',temperature%pth
c      PRINT *,'temperature%pph = ',temperature%pph

      rho%ptmtm = rhov * (((1.d0/rhov)*rho%ptm)**2 +
     & (1.d0/p)*pressure%ptmtm - ((1.d0/p)*pressure%ptm)**2 -
     & (1.d0/T)*temperature%ptmtm + ((1.d0/T)*temperature%ptm)**2 +
     & (1.d0/mw)*molecular_weight%ptmtm -
     & ((1.d0/mw)*molecular_weight%ptm)**2)
      rho%prtm = rhov * (((1.d0/(rhov**2))*rho%pr*rho%ptm) +
     & (1.d0/p)*pressure%prtm-((1.d0/(p**2))*pressure%pr*pressure%ptm)-
     & (1.d0/T)*temperature%prtm +
     & ((1.d0/(T**2))*temperature%pr*temperature%ptm) +
     & (1.d0/mw)*molecular_weight%prtm -
     & ((1.d0/(mw**2))*molecular_weight%pr*molecular_weight%ptm))
      rho%pthtm = rhov * (((1.d0/(rhov**2))*rho%pth*rho%ptm) +
     & (1.d0/p)*pressure%pthtm -
     & ((1.d0/(p**2))*pressure%pth*pressure%ptm) -
     & (1.d0/T)*temperature%pthtm +
     & ((1.d0/(T**2))*temperature%pth*temperature%ptm) +
     & (1.d0/mw)*molecular_weight%pthtm -
     & ((1.d0/(mw**2))*molecular_weight%pth*molecular_weight%ptm))
      rho%pphtm = rhov * (((1.d0/(rhov**2))*rho%pph*rho%ptm) +
     & (1.d0/p)*pressure%pphtm -
     & ((1.d0/(p**2))*pressure%pph*pressure%ptm) -
     & (1.d0/T)*temperature%pphtm +
     & ((1.d0/(T**2))*temperature%pph*temperature%ptm) +
     & (1.d0/mw)*molecular_weight%pphtm -
     & ((1.d0/(mw**2))*molecular_weight%pph*molecular_weight%ptm))

      rho%ptmr = rhov * ((rho%ptm*rho%pr/(rhov**2)) +
     & (pressure%ptmr/p) - (pressure%ptm*pressure%pr/(p**2)) -
     & (temperature%ptmr/T) + (temperature%ptm*temperature%pr/(T**2)) +
     & (molecular_weight%ptmr/mw) -
     & (molecular_weight%ptm*molecular_weight%pr/(mw**2)))
      rho%prr = rhov * ((rho%pr/rhov)**2 + (pressure%prr/p) -
     & (pressure%pr/p)**2 - (temperature%prr/T) +
     & (temperature%pr/T)**2 + (molecular_weight%prr/mw) +
     & (molecular_weight%pr/mw)**2)
      rho%pthr = rhov * ((rho%pth*rho%pr/(rhov**2)) +
     & (pressure%pthr/p) - (pressure%pth*pressure%pr/(p**2)) -
     & (temperature%pthr/T) + (temperature%pth*temperature%pr/(T**2)) +
     & (molecular_weight%pthr/mw) -
     & (molecular_weight%pth*molecular_weight%pr/(mw**2)))
      rho%pphr = rhov * (((1.d0/(rhov**2))*rho%pph*rho%pr) +
     & (1.d0/p)*pressure%pphr-((1.d0/(p**2))*pressure%pph*pressure%pr)-
     & (1.d0/T)*temperature%pphr +
     & ((1.d0/(T**2))*temperature%pph*temperature%pr) +
     & (1.d0/mw)*molecular_weight%pphr -
     & ((1.d0/(mw**2))*molecular_weight%pph*molecular_weight%pr))

c      somethingpr1 = pressure%pr/p
c      somethingpr2 = pressure%prr/p - (pressure%pr/p)**2
c      PRINT *,'somethingpr1 = ',somethingpr1
c      PRINT *,'somethingpr2 = ',somethingpr2

c      term1 = (rho%pr/rhov)**2
c      term2 = (pressure%prr/p)
c      term3 = (pressure%pr/p)**2 
c      term4 = (temperature%prr/T)
c      term5 = (temperature%pr/T)**2
c      PRINT *,'term1 = ',term1
c      PRINT *,'term2 = ',term2
c      PRINT *,'term3 = ',term3
c      PRINT *,'term4 = ',term4
c      PRINT *,'term5 = ',term5

c      height = r-earthr
c      PRINT *,'rhov = ',rhov
c      PRINT *,'rho%pr = ',rho%pr
c      PRINT *,'rho%prr = ',rho%prr
c      PRINT *,'pressure%prr = ',pressure%prr
c      PRINT *,'pressure%pr = ',pressure%pr
c      PRINT *,'height = ',height
c      PRINT *,'p = ',p
c      PRINT *,'temperature%prr = ',temperature%prr
c      PRINT *,'temperature%pr = ',temperature%pr
c      PRINT *,'T = ',T
c      PRINT *,'molecular_weight%prr = ',molecular_weight%prr
c      PRINT *,'molecular_weight%pr = ',molecular_weight%pr
c      PRINT *,'mw = ',mw

      rho%ptmth = rhov * (((1.d0/(rhov**2))*rho%ptm*rho%pth) +
     & (1.d0/p)*pressure%ptmth -
     & ((1.d0/(p**2))*pressure%ptm*pressure%pth)-
     & (1.d0/T)*temperature%ptmth +
     & ((1.d0/(T**2))*temperature%ptm*temperature%pth) +
     & (1.d0/mw)*molecular_weight%ptmth -
     & ((1.d0/(mw**2))*molecular_weight%ptm*molecular_weight%pth))
      rho%prth = rhov * (((1.d0/(rhov**2))*rho%pr*rho%pth) +
     & (1.d0/p)*pressure%prth-((1.d0/(p**2))*pressure%pr*pressure%pth)-
     & (1.d0/T)*temperature%prth +
     & ((1.d0/(T**2))*temperature%pr*temperature%pth) +
     & (1.d0/mw)*molecular_weight%prth -
     & ((1.d0/(mw**2))*molecular_weight%pr*molecular_weight%pth))
      rho%pthth = rhov * (((1.d0/rhov)*rho%pth)**2 +
     & (1.d0/p)*pressure%pthth - ((1.d0/p)*pressure%pth)**2 -
     & (1.d0/T)*temperature%pthth + ((1.d0/T)*temperature%pth)**2 +
     & (1.d0/mw)*molecular_weight%pthth -
     & ((1.d0/mw)*molecular_weight%pth)**2)
      rho%pphth = rhov * (((1.d0/(rhov**2))*rho%pph*rho%pth) +
     & (1.d0/p)*pressure%pphth -
     & ((1.d0/(p**2))*pressure%pph*pressure%pth) -
     & (1.d0/T)*temperature%pphth +
     & ((1.d0/(T**2))*temperature%pph*temperature%pth) +
     & (1.d0/mw)*molecular_weight%pphth -
     & ((1.d0/(mw**2))*molecular_weight%pph*molecular_weight%pth))

      rho%ptmph = rhov * (((1.d0/(rhov**2))*rho%ptm*rho%pph) +
     & (1.d0/p)*pressure%ptmph -
     & ((1.d0/(p**2))*pressure%ptm*pressure%pph) -
     & (1.d0/T)*temperature%ptmph +
     & ((1.d0/(T**2))*temperature%ptm*temperature%pph) +
     & (1.d0/mw)*molecular_weight%ptmph -
     & ((1.d0/(mw**2))*molecular_weight%ptm*molecular_weight%pph))
      rho%prph = rhov * (((1.d0/(rhov**2))*rho%pr*rho%pph) +
     & (1.d0/p)*pressure%prph-((1.d0/(p**2))*pressure%pr*pressure%pph)-
     & (1.d0/T)*temperature%prph +
     & ((1.d0/(T**2))*temperature%pr*temperature%pph) +
     & (1.d0/mw)*molecular_weight%prph -
     & ((1.d0/(mw**2))*molecular_weight%pr*molecular_weight%pph))
      rho%pthph = rhov * (((1.d0/(rhov**2))*rho%pth*rho%pph) +
     & (1.d0/p)*pressure%pthph -
     & ((1.d0/(p**2))*pressure%pth*pressure%pph) -
     & (1.d0/T)*temperature%pthph +
     & ((1.d0/(T**2))*temperature%pth*temperature%pph) +
     & (1.d0/mw)*molecular_weight%pthph -
     & ((1.d0/(mw**2))*molecular_weight%pth*molecular_weight%pph))
      rho%pphph = rhov * (((1.d0/rhov)*rho%pph)**2 +
     & (1.d0/p)*pressure%pphph - ((1.d0/p)*pressure%pph)**2 -
     & (1.d0/T)*temperature%pphph + ((1.d0/T)*temperature%pph)**2 +
     & (1.d0/mw)*molecular_weight%pphph -
     & ((1.d0/mw)*molecular_weight%pph)**2)

c      tempor = rho%pr/rhov
c      PRINT *,'tempor = ',tempor
c      tempor2 = rho%prr/rhov - (rho%pr/rhov)**2
c      tempor2 = - (temperature%prr/T) +
c     & (temperature%pr/T)**2 + (molecular_weight%prr/mw) +
c     & (molecular_weight%pr/mw)**2
c      PRINT *,'tempor2 = ',tempor2

c      PRINT *,'molecular_weight%pr = ',molecular_weight%pr
c      PRINT *,'molecular_weight%prr = ',molecular_weight%prr

c      PRINT *,'temperature%pr = ',temperature%pr
c      PRINT *,'temperature%prr = ',temperature%prr

c      PRINT *,'pressure%pr = ',pressure%pr
c      PRINT *,'pressure%prr = ',pressure%prr

c      PRINT *,'rho%ptm = ',rho%ptm
c      PRINT *,'rho%pr = ',rho%pr
c      PRINT *,'rho%pth = ',rho%pth
c      PRINT *,'rho%pph = ',rho%pph
c      PRINT *,'rho%ptmtm = ',rho%ptmtm
c      PRINT *,'rho%prtm = ',rho%prtm
c      PRINT *,'rho%pthtm = ',rho%pthtm
c      PRINT *,'rho%pphtm = ',rho%pphtm
c      PRINT *,'rho%ptmr = ',rho%ptmr
c      PRINT *,'rho%prr = ',rho%prr
c      PRINT *,'rho%pthr = ',rho%pthr
c      PRINT *,'rho%pphr = ',rho%pphr
c      PRINT *,'rho%ptmth = ',rho%ptmth
c      PRINT *,'rho%prth = ',rho%prth
c      PRINT *,'rho%pthth = ',rho%pthth
c      PRINT *,'rho%pphth = ',rho%pphth
c      PRINT *,'rho%ptmph = ',rho%ptmph
c      PRINT *,'rho%prph = ',rho%prph
c      PRINT *,'rho%pthph = ',rho%pthph
c      PRINT *,'rho%pphph = ',rho%pphph

      potentialrho%value = rho%value
      potentialrho%ptm = rho%ptm - pressure%ptm/cs
      potentialrho%pr = rho%pr - pressure%pr/cs
      potentialrho%pth = rho%pth - pressure%pth/cs
      potentialrho%pph = rho%pph - pressure%pph/cs

c      print*,'cs = ',cs

      potentialrho%ptmtm = potentialrho%ptm**2/rhov + rho%ptmtm -
     & rho%ptm**2/rhov - pressure%ptmtm/cs +
     & (pressure%ptm*pcst)/(cs**2) + (pressure%ptm*rho%ptm)/(rhov*cs)
      potentialrho%prtm = (potentialrho%pr*rho%ptm)/rhov + rho%prtm -
     & (rho%pr*rho%ptm)/rhov - (pressure%pr*pressure%ptm)/cs +
     & (pressure%pr*pcst)/(cs**2) + (pressure%pr*rho%ptm)/(rhov*cs)
      potentialrho%pthtm = (potentialrho%pth*rho%ptm)/rhov +
     & rho%pthtm - (rho%pth*rho%ptm)/rhov -
     & (pressure%pth*pressure%ptm)/cs + (pressure%pth*pcst)/(cs**2) +
     & (pressure%pth*rho%ptm)/(rhov*cs)
      potentialrho%pphtm = (potentialrho%pph*rho%ptm)/rhov +
     & rho%pphtm - (rho%pph*rho%ptm)/rhov -
     & (pressure%pph*pressure%ptm)/cs + (pressure%pph*pcst)/(cs**2) +
     & (pressure%pph*rho%ptm)/(rhov*cs)

      potentialrho%ptmr = (potentialrho%ptm*rho%pr)/rhov + rho%ptmr -
     & (rho%ptm*rho%pr)/rhov - (pressure%ptm*pressure%pr)/cs +
     & (pressure%ptm*pcsr)/(cs**2) + (pressure%ptm*rho%pr)/(rhov*cs)
      potentialrho%prr = potentialrho%pr*rho%pr/rhov + rho%prr -
     & rho%pr**2/rhov - pressure%prr/cs + (pressure%pr*pcsr)/(cs**2) +
     & (pressure%pr*rho%pr)/(rhov*cs)
      potentialrho%pthr = (potentialrho%pth*rho%pr)/rhov +
     & rho%pthr - (rho%pth*rho%pr)/rhov -
     & (pressure%pth*pressure%pr)/cs + (pressure%pth*pcsr)/(cs**2) +
     & (pressure%pth*rho%pr)/(rhov*cs)
      potentialrho%pphr = (potentialrho%pph*rho%pr)/rhov +
     & rho%pphr - (rho%pph*rho%pr)/rhov -
     & (pressure%pph*pressure%pr)/cs + (pressure%pph*pcsr)/(cs**2) +
     & (pressure%pph*rho%pr)/(rhov*cs)

      potentialrho%ptmth = (potentialrho%ptm*rho%pth)/rhov +
     & rho%ptmth - (rho%ptm*rho%pth)/rhov -
     & (pressure%ptm*pressure%pth)/cs + (pressure%ptm*pcsth)/(cs**2) +
     & (pressure%ptm*rho%pth)/(rhov*cs)
      potentialrho%prth = (potentialrho%pr*rho%pth)/rhov + rho%prth -
     & (rho%pr*rho%pth)/rhov - (pressure%pr*pressure%pth)/cs +
     & (pressure%pr*pcsth)/(cs**2) + (pressure%pr*rho%pth)/(rhov*cs)
      potentialrho%pthth = potentialrho%pth*rho%pth/rhov + rho%pthth -
     & rho%pth**2/rhov - pressure%pthth/cs +
     & (pressure%pth*pcsth)/(cs**2) + (pressure%pth*rho%pth)/(rhov*cs)
      potentialrho%pphth = (potentialrho%pph*rho%pth)/rhov +
     & rho%pphth - (rho%pph*rho%pth)/rhov -
     & (pressure%pph*pressure%pth)/cs + (pressure%pph*pcsth)/(cs**2) +
     & (pressure%pph*rho%pth)/(rhov*cs)

      potentialrho%ptmph = (potentialrho%ptm*rho%pph)/rhov +
     & rho%ptmph - (rho%ptm*rho%pph)/rhov -
     & (pressure%ptm*pressure%pph)/cs + (pressure%ptm*pcsph)/(cs**2) +
     & (pressure%ptm*rho%pph)/(rhov*cs)
      potentialrho%prph = (potentialrho%pr*rho%pph)/rhov + rho%prph -
     & (rho%pr*rho%pph)/rhov - (pressure%pr*pressure%pph)/cs +
     & (pressure%pr*pcsph)/(cs**2) + (pressure%pr*rho%pph)/(rhov*cs)
      potentialrho%pthph = (potentialrho%pth*rho%pph)/rhov +
     & rho%pthph - (rho%pth*rho%pph)/rhov -
     & (pressure%pth*pressure%pph)/cs + (pressure%pth*pcsph)/(cs**2) +
     & (pressure%pth*rho%pph)/(rhov*cs)
      potentialrho%pphph = potentialrho%pph*rho%pph/rhov + rho%pphph -
     & rho%pph**2/rhov - pressure%pphph/cs +
     & (pressure%pph*pcsph)/(cs**2) + (pressure%pph*rho%pph)/(rhov*cs)

c      testing1 = potentialrho%pr/rhov
c      testing2 = potentialrho%prr/rhov -
c     & (rho%pr*potentialrho%pr)/(rhov**2)

c      PRINT *,'testing1 = ',testing1
c      PRINT *,'testing2 = ',testing2
c      PRINT *,' = '

      CALL PDENSITYg(rho,potentialrho)                            ! 25 Mar 2010
      RETURN

C ---------------------------------------------------------------------

      ENTRY DENSITYa

c      print*,'Beginning of DENSITYa'
      CALL temp
c      print*,'After call to temp'
      CALL molwt
c      print*,'After call to molwtg'
      CALL pres
c      print*,'After call to pres'

c      print*,'p = ',p
c      print*,'m = ',m
c      print*,'T = ',T
c      print*,'rgas = ',rgas

c      hscale = -p/ppr

      ds = (p*m)/(rgas*T)                                        ! 22 Mar 2010
c      print*,'ds = ',ds

      pdst = ds * ((ppt/p) - (pTt/T) + (pmt/m))                  ! 22 Mar 2010
      pdsr = ds * ((ppr/p) - (pTr/T) + (pmr/m))                  ! 22 Mar 2010
      pdsth = ds * ((ppth/p) - (pTth/T) + (pmth/m))              ! 22 Mar 2010
      pdsph = ds * ((ppph/p) - (pTph/T) + (pmph/m))              ! 22 Mar 2010

c      print*,'End of DENSITYa'

      RETURN

C ---------------------------------------------------------------------
      END SUBROUTINE RHOSMP
C ---------------------------------------------------------------------

C ---------------------------------------------------------------------
      SUBROUTINE setdensity
C ---------------------------------------------------------------------
C     This subroutine is called once.

      IMPLICIT double precision (a-h,o-z)

C     COMMON DECK "B15" INSERTED HERE
      INTEGER PRMX,PRNTBL,PRITBL,PRFRMTB,IDSPR(10)
      COMMON/B15/PRMX,PRNTBL(10),PRITBL(10),PRFRMTB(10),PRGP(11)
      EQUIVALENCE (PRGP,IDSPR),(ANP,PRGP(11))
      INTEGER PFMX,PFNTBL(10),PFITBL(10),PFFRMTB(10)
      DATA PFMX/1/
      DATA PFNTBL/1,11,8*0/
      DATA PFITBL/1,9*0/
      DATA PFFRMTB/1,9*0/
      DATA RECOGP/1.d0/

c      print*,'Beginning of setdensity'

      PRMX=PFMX
      CALL imove(PRNTBL,PFNTBL,10)  
      CALL imove(PRITBL,PFITBL,10)  
      CALL imove(PRFRMTB,PFFRMTB,10)

      CALL setpdensity                                          ! 24 Jan 2006
      RETURN                        

C ---------------------------------------------------------------------
      END SUBROUTINE setdensity
C ---------------------------------------------------------------------


C ---------------------------------------------------------------------
      SUBROUTINE idensity
C ---------------------------------------------------------------------
C     This subroutine is called once for each new W array.

      IMPLICIT double precision (a-h,o-z)

      TYPE :: model
         CHARACTER :: model_name
         REAL(KIND=8) :: data_id
         CHARACTER :: pert_name
         REAL(KIND=8) :: pert_id
      END TYPE model
      TYPE (model) :: moddg

      PARAMETER (NWARSZ=1000)
C     COMMON DECK "WW" INSERTED HERE
      COMMON/WW/dum(10),MAXW,W(NWARSZ)

C     DENSITY     650-674
      EQUIVALENCE (W(650),dmodel),(W(651),dform),(W(652),did)

C     COMMON DECK "DS" INSERTED HERE
      REAL(KIND=8) ds,modds
      COMMON/DS/MODDS(4),ds,pdst,pdsr,pdsth,pdsph

C     COMMON DECK "B15" INSERTED HERE
c      INTEGER PRMX,PRNTBL,PRITBL,PRFRMTB,IDSPR(10)
c      COMMON/B15/PRMX,PRNTBL(10),PRITBL(10),PRFRMTB(10),PRGP(11)
c      EQUIVALENCE (PRGP,IDSPR),(ANP,PRGP(11))
c      INTEGER PFMX,PFNTBL(10),PFITBL(10),PFFRMTB(10)
c      DATA PFMX/1/
c      DATA PFNTBL/1,11,8*0/
c      DATA PFITBL/1,9*0/
c      DATA PFFRMTB/1,9*0/
c      DATA RECOGP/1.d0/

c      print*,'Beginning of idensity'

      IF (dmodel /= 1.d0) THEN
           PRINT *,' Model check number mismatch: '
           STOP ' W(650) should be 1 if you are using model RHOSMP.'
      ENDIF

      modds(1) = 'RHOSMP'
      modds(2) = did                                            ! 19 Jul 2006

c      print*,'Before call ipdensity'
      CALL ipdensity                                            ! 19 Jul 2006
c      print*,'After call ipdensity'

      RETURN                                                    ! 19 Jul 2006

C ---------------------------------------------------------------------

      ENTRY idensityg(moddg)                                    ! 19 Jul 2006
      IF (dmodel /= 1.d0) THEN
           PRINT *,' Model check number mismatch: '
           STOP ' W(650) should be 1 if you are using model RHOSMP.'
      ENDIF

      moddg%model_name = 'RHOSMP'
      moddg%data_id = did
      moddg%pert_name = ''
      moddg%pert_id = 0.d0

!      modds(1) = moddg%model_name                                ! 30 Apr 2010
!      modds(2) = moddg%data_id                                   ! 30 Apr 2010
      modds(1) = 'RHOSMP'                                         ! 30 Apr 2010
      modds(2) = did                                ! 19 Jul 2006 ! 30 Apr 2010

      CALL ipdensityg(moddg)                                    ! 19 Jul 2006

      RETURN
C ---------------------------------------------------------------------
      END SUBROUTINE idensity
C ---------------------------------------------------------------------
