subroutine airrhc(sst,qa,rh,qss,qsa,ta,da)
CC known:sst, qa (humidity mixing ratio),  constant rh 
CC calculated values: qss (saturated q for sea surface),
CC             qsa (sat. q for air), ta (air temp), da (air density)

      PARAMETER (airmol=28.97d0, r0=8.3143D+3 )
      DATA PI/3.141592654D0/
      DATA PR/1.01325e5/

c      if(rh.gt.10) rh=rh/100.0
     
CC saturated vapor pressure ess & mixing ratio wss at sea surface
       ess=6.11d+2*dexp((2.5d+6/461.0d0)*(1.0d0/273.15d0-1.0/sst))
       wss=0.622*ess/(PR-ess)

       qss=wss

CC     mixing ratio w & saturated mixing ratio ws
       w=qa
       ws=w/rh
       qsa=qa/rh

       es=ws*PR/(0.622+ws)
CC  for saturated mixing ratio:  ws=0.622*es/(P-es)
       ta1=alog(es/6.11e+2)/(2.5d+6/461.0d0)
       ta=1.0d0/273.15d0-ta1
       ta=1.0/ta

       tva=ta*(1.0+0.6077*w)

       da=PR*airmol/(r0*tva)
c      air's mol number in unit volumn n=PV/RT
c                P in unit Pa.
c      airmol is air's kmol weight,and d is in kg/m**3

      return
      END


      subroutine airtac(sst,qa,ta,qss,qsa,rh,da)
CC known:sst, qa (humidity mixing ratio),  ta (air temp) 
CC calculated values: qss (saturated q for sea surface),
CC       qsa (sat. q for air), rh (relative humidity), da (air density)

      PARAMETER (airmol=28.97d0, r0=8.3143D+3 )
      DATA PI/3.141592654D0/
      DATA PR/1.01325e5/

CC saturated vapor pressure ess & mixing ratio qss at sea surface
       ess=6.11d+2*dexp((2.5d+6/461.0d0)*(1.0d0/273.15d0-1.0/sst))
       qss=0.622*ess/(PR-ess)

CC     mixing ratio w & saturated mixing ratio ws
       w=qa
       es=6.11d+2*dexp((2.5d+6/461.0d0)*(1.0d0/273.15d0-1.0/ta))
       ws=0.622*es/(PR-es)

       rh=w/ws
       qsa=ws

       tva=ta*(1.0+0.6077*w)

       da=PR*airmol/(r0*tva)
c      air's mol number in unit volumn n=PV/RT
c                P in unit Pa.
c      airmol is air's kmol weight,and d is in kg/m**3

      return
      END