c     kspectrum (http://www.meso-star.com/en_Products.html) - This file is part of kspectrum
c     Copyright (C) 2008-2015 - Méso-Star - Vincent Eymet
c     
c     This file must be used under the terms of the CeCILL license.
c     This source file is licensed as described in the file COPYING, which
c     you should have received as part of this distribution.  The terms
c     are also available at
c     http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
c     
      subroutine density(Rgp,P,Ps,T,iso_abundance,
     &     dens,drho_dP,drho_dT,drho_dx)
      implicit none
      include 'max.inc'
c     
c     Purpose: to compute the density (in molecules/cm³) of the given molecular species
c     
c     Inputs:
c       + Rgp: perfect gas constant [atm.m³/(mol.K)]
c       + P: total pressure [atm]
c       + Ps: partial pressure of the molecular species [atm]
c       + T: temperature [K]
c       + iso_abundance is the abundance of isotope 'iso' for molecule 'mol'
c
c     Outputs:
c       + dens: density of the given species [molecule.cm¯³]
c       + drho_dP: sentivity of density to pressure [molecule.cm¯³.atm¯¹]
c       + drho_dT: sentivity of density to temperature [molecule.cm¯³.K¯¹]
c       + drho_dx: sentivity of density to the concentration [molecule.cm¯³]
c

      include 'parameters.inc'
c     input
      double precision Rgp,P,Ps,T
      double precision iso_abundance
c     output
      double precision dens
      double precision drho_dP
      double precision drho_dT
      double precision drho_dx
c     temp
      double precision x
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine density'

c     WARNING: perfect gas hyp.
c
c     Modification - 2013-03-06
c     The reference line intensiy already takes into account the abundance of isotopologues
c      dens=Nav*Ps/(Rgp*T)*iso_abundance ! [molecule/mol]*[atm]/[atm.m³/(mol.K)]/[K]=[molecule.m¯³]
      dens=Nav*Ps/(Rgp*T) ! [molecule/mol]*[atm]/[atm.m³/(mol.K)]/[K]=[molecule.m¯³]
c     Modification - 2013-03-06
c
      dens=dens*1.0D-6 ! [molecule.cm¯³]

c     sensitivities
c     + to pressure
      if (P.ne.0.0D+0) then
         drho_dP=dens/P
      else
         write(*,*) 'Error from ',label(1:strlen(label)),' :'
         write(*,*) 'P=',P
      endif
c     + to temperature
      if (T.ne.0.0D+0) then
         drho_dT=-dens/T
      else
         write(*,*) 'Error from ',label(1:strlen(label)),' :'
         write(*,*) 'T=',T
      endif
c     + to the concentration of the molecule
      x=Ps/P ! concentration
      if (x.ne.0.0D+0) then
         drho_dx=dens/x
      else
         write(*,*) 'Error from ',label(1:strlen(label)),' :'
         write(*,*) 'x=',x
      endif
      
      return
      end