C **********************************************************************
C *                                                                    *
C *                                SOFTWARE LICENSING                  *
C *                                                                    *
C *                This program is free software; you can redistribute *
C *                it and/or modify it under the terms of the GNU      *
C *                General Public License as published by the Free     *
C *                Software Foundation, either Version 2 of the        *
C *                license, or (at your option) any later version.     *
C *                                                                    *
C *                This program is distributed in the hope that it     *
C *                will be useful, but without any warranty; without   *
C *                even the implied warranty of merchantability or     *
C *                fitness for a particular purpose. See the GNU       *
C *                General Public License for more details.            *
C *                                                                    *
C *                A copy of the GNU General Public License is         *
C *                available at http://www.gnu.org/copyleft/gpl.html   *
C *                or by writing to the Free Software Foundation, Inc.,*
C *                59 Temple Place - Suite 330, Boston, MA 02111, USA. *
C *                                                                    *
C **********************************************************************

      FUNCTION SADENS(BAROP,AIRTMP,RELHUM) ! NG12022009 
CNG
CNG  ESTIMATE SURFACE AIR DENSITY (kg/m3) BASED ON:
CNG
CNG  BAROP: Barometric Pressure, mbar
CNG  AIRTMP: Surface air temperature, C
CNG  RELHUM: Relative Humidity, %
CNG

      DATA RMOLRAT /0.62197/ ! Molecular weight ratio
      EW=6.1121*(1.0007+3.46E-6*BAROP) *
     .   EXP( (18.677 - AIRTMP/234.50) * (AIRTMP/(257.14 + AIRTMP)) ) !  saturation vapor pressure, new Buck 1996 formula
      QSAT=RMOLRAT*EW/(BAROP-(1.-RMOLRAT)*EW) ! specific humidity at saturation

      Q=(0.01*RELHUM)*QSAT ! vapor pressure (specific humidity of air)
      TT=AIRTMP+273.16
      O61=1./RMOLRAT-1. ! =0.61, moisture correction for temperature
      TV=TT*(1+O61*Q) ! virtual temperature
      SADENS=100*BAROP/(287.04*TV+1.e-30) ! 287.04 is the gas constant for dry air
      SADENS=AMAX1(AMIN1(SADENS,1.6),0.9) ! Make sure they are within "reasonable" range at earth's sea level
      RETURN
      END