C***********************************************************************
C  FORTRAN77 CODE TO CALCULATE THE DENSITY OF
C  MULTICOMPONENT AQUEOUS SOLUTIONS CONTAINING
C  H+/NH4+/SO42–/NO3–/H2O
C
C  BY MICHA SEMMLER
C  AND THOMAS KOOP <thomas.koop@uni-bielefeld.de>
C  Version of March 31, 2006.
C
C  IF YOU USE THIS CODE FOR YOUR SCIENTIFIC WORK,
C  PLEASE GIVE CREDIT BY CITING THE FOLLOWING PUBLICATION:
C
C  SEMMLER, M, LUO, B.P., AND KOOP, T.,
C  Densities of liquid H+/NH4+/SO42–/NO3–/H2O solutions
C  at tropospheric temperatures,
C  Atmos. Environment, 40, 467-483, 2006.
C***********************************************************************
C
C  INPUT:
C  temp = temperature in [K]
C  wt2 = mass fraction of sulfuric acid (H2SO4)
C  wt3 = mass fraction of ammonium sulfate ((NH4)2SO4)
C  wt4 = mass fraction of ammonium nitrate (NH4NO3)
C  wt5 = mass fraction of nitric acid (HNO3)
C     Note: all mass fraction values are 0.0 < wt < 1.0
C
C  VALIDITY RANGE:
C      253K  =< temp =< 298K
C      0.1   =< wt2+wt3+wt3+wt4  =<  0.6
C
C  THE MODEL IS ASSSUMED TO PROVIDE REASONABLE RESULTS OVER
C      180K   =< temp =< 330K with
C      0.1   =< wt2+wt3+wt3+wt4  =<  0.6
C
C  AND SHOULD NOT BE USED OUTSIDE THIS RANGE OF VALUES. SEE
C  THE ABOVE PAPER FOR A MORE DETAILED DISCUSSION.
C
C  OUTPUT:
C  rho  = density of liquid aqueous solution in [g cm^-3]
C
C***********************************************************************
C
      PROGRAM density_loop
      IMPLICIT REAL*8 (A-H,O-Z)
C
C***********************************************************************
C
C   MAIN PROGRAM
C
      OPEN(11,FILE='temp.dat',STATUS='unknown')
 98   FORMAT(1x,6(F12.5))
C
C   INPUT
C
       wt2 = 0.0
       wt3 = 0.0
       wt4 = 0.0
       wt5 = 0.0
C
      tempstart=340.d0
      tempend=170.d0
      tempstep=-5.d0
C
      wtstart=0.00d0
      wtend=0.9d0
      wtstep=0.1d0
C
      DO temp=tempstart,tempend,tempstep
      DO wt2=wtstart,wtend,wtstep
C
        dens=rho(wt2,wt3,wt4,wt5,temp)
C
C   OUTPUT
C
        write(11,98)temp,wt2,wt3,wt4,wt5,dens
C
      ENDDO
      ENDDO
      close(11)
      END
C
C***********************************************************************
c
c
      function rho(wts0,wtn0,wta0,wtb0,temp)
c
C
C     Returns -99 for T < 180K and T > 330 K
C     and for total solute mass fraction wlim < 0.1 and wlim > 0.6
C
       implicit real*8 (a-h,o-z)
       parameter (nP=9,nS=5)
       real*8 x(9,5),wm(6),z(3)
       data z/0.178852E-01,0.181472E-02,0.000000E-00/
       data wm/18.015,98.078,132.140,80.043,63.012,115.109/
       data (x(i,1), i=1,nP)
     *      /0.362249E-01, 0.546891E-02, 0.873468E-02,
     *       0.344556E-01,-0.650688E-01, 0.417174E-01,
     *       0.136245E-01,-0.391920E-01, 0.252329E-01/
       data (x(i,2), i=1,nP)
     *      /0.558085E-01, 0.406200E-01,-0.110797E-01,
     *       0.135837E-01,-0.158075E-01, 0.000000E+00,
     *       0.000000E+00, 0.000000E+00, 0.000000E+00/
       data (x(i,3), i=1,nP)
     *      /0.503779E-01,-0.691069E-02, 0.113024E-01,
     *      -0.287432E-01, 0.155419E+00,-0.157049E+00,
     *       0.161094E+00,-0.585880E+00, 0.541843E+00/
       data (x(i,4), i=1,nP)
     *      /0.287384E-01,-0.136204E-02, 0.159091E-01,
     *       0.150333E-01, 0.313911E-03, 0.000000E+00,
     *       0.000000E+00, 0.000000E+00, 0.000000E+00/
       data (x(i,5), i=1,nP)
     *      /0.476267E-01, 0.282931E-01,-0.890139E-02,
     *       0.343616E-01,-0.458649E-01, 0.000000E+00,
     *       0.000000E+00, 0.000000E+00, 0.000000E+00/
c
c     Temperature is converted to reduced temperature relative to 273.15K
c
       t=(temp-273.15)/273.15
c
c
       w20=wts0
       w30=wtn0
c
       w4=wta0
       w5=wtb0
c
       w=w20+w30+w4+w5
       w1=1.-w
c
        xm10=1d3/wm(1)
        xm20=1d3/wm(2)*w20/w1
        xm30=1d3/wm(3)*w30/w1
c
        xm6 = min(xm20,xm30)
        xm2 = xm20-xm6
        xm3 = xm30-xm6
c
c      Attention: ONE AmmoSulf + ONE SulfAcid -> TWO BiSulf !!
        xm6 = 2d0*xm6
c
        w2=xm2*wm(2)*w1/1d3
        w3=xm3*wm(3)*w1/1d3
        w6=xm6*wm(6)*w1/1d3
c
        v1=z(1)+z(2)*t+z(3)*t**2
c
        v20=x(1,1)+x(2,1)*w+x(3,1)*w**2
        v21=x(4,1)+x(5,1)*w+x(6,1)*w**2
        v22=x(7,1)+x(8,1)*w+x(9,1)*w**2
c
        v2=v20 + v21*t + v22*t**2
c
        v30=x(1,2)+x(2,2)*w+x(3,2)*w**2
        v31=x(4,2)+x(5,2)*w+x(6,2)*w**2
        v32=x(7,2)+x(8,2)*w+x(9,2)*w**2
c
        v3=v30 + v31*t + v32*t**2
c
        v40=x(1,3)+x(2,3)*w+x(3,3)*w**2
        v41=x(4,3)+x(5,3)*w+x(6,3)*w**2
        v42=x(7,3)+x(8,3)*w+x(9,3)*w**2
c
        v4=v40 + v41*t + v42*t**2
c
        v50=x(1,4)+x(2,4)*w+x(3,4)*w**2
        v51=x(4,4)+x(5,4)*w+x(6,4)*w**2
        v52=x(7,4)+x(8,4)*w+x(9,4)*w**2
c
        v5=v50 + v51*t + v52*t**2
c
        v60=x(1,5)+x(2,5)*w+x(3,5)*w**2
        v61=x(4,5)+x(5,5)*w+x(6,5)*w**2
        v62=x(7,5)+x(8,5)*w+x(9,5)*w**2
c
        v6=v60 + v61*t + v62*t**2
c
        vmcal=w1*v1/wm(1) + w2*v2/wm(2) +
     *        w3*v3/wm(3) + w4*v4/wm(4) +
     *        w5*v5/wm(5) + w6*v6/wm(6)
c
c    Now set limits for range of validity
c
      wlim=wts0+wtn0+wta0+wtb0
      tlim=t*273.15+273.15
c
      IF(Tlim.LT.180.d0 .OR. Tlim.GT.330.d0)THEN
            rho=-99.d0
      ELSEIF(wlim.LT.0.1d0 .OR. wlim.GT.0.6d0)THEN
            rho=-99.d0
      ELSE
            rho=1/vmcal/1000.d0
      ENDIF
c
       return
      end
C***********************************************************************
