C***********************************************************************
C  FORTRAN77 CODE TO CALCULATE THE HOMOGENEOUS ICE NUCLEATION
C  RATE COEFFICIENT FROM LIQUID AQUEOUS SOLUTIONS
C  FROM 170 K TO 240 K.
C
C  BY THOMAS KOOP <thomas.koop@uni-bielefeld.de>
C  Version of March 6, 2006.
C
C  IF YOU USE THIS CODE FOR YOUR SCIENTIFIC WORK, 
C  PLEASE GIVE CREDIT BY CITING THE FOLLOWING PUBLICATION:
C
C  KOOP, T., Luo, B.P., Tsias, A., Peter, T.,
C  Water activity as the determinant for homogeneous 
C  ice nucleation in aqueous solutions,
C  Nature 406, 611-614, 2000.  
C***********************************************************************
C
C  INPUT: 
C  t = temperature of liquid;   range of validity: 170 to 240 K
C  awliq = water activity of liquid;   range of validity: 0.7 to 1.0
C         Note: Calculate the water activity from the ratio of the
C         water vapor pressure of an aerosol droplets and the     
C         water vapor pressure of pure water at the same temperature.
C         DO NOT USE AMBIENT RELATIVE HUMIDITY TO CALCULATE WATER
C         ACTIVITY, SINCE IN STRONG UPDRAFTS THE TWO ARE NOT IDENTICAL!
C
C  OUTPUT:
C  jhom  = homogeneous ice nucleation rate coefficient in [cm^-3 s^-1]
C  
C***********************************************************************
C
      PROGRAM icenucleation_ambientp
      IMPLICIT REAL*8(A-H,O-Z)
C
C***********************************************************************
C
C   MAIN PROGRAM
C
      OPEN(11,FILE='temp.dat',STATUS='unknown')
 98   FORMAT(1x,2(F12.6),1(E15.7))
C
C   INPUT
C
      tstart=245d0
      tend=165.d0
      tstep=-1.0d0 
C      
      print*, ' input water activity of liquid aqueous solution'
      print*, ' (awliq value must be between 0.7 and 1.0):'
      read(5,*)awliq
C
      DO t=tstart,tend,tstep
C
        xxx=xjhom(t,awliq)  
C
C   OUTPUT
C
        write(11,98)awliq,t,xxx
C
      ENDDO
      close(11)
      END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C---------------------------------------
C    ICE NUCLEATION RATE COEFFICIENT 
C---------------------------------------
C
      FUNCTION xjhom(t,awliq)
      IMPLICIT REAL*8(A-H,O-Z)
C
C     Returns -99 for T < 170K and T > 240 K.
C     and 
C     1.e19 for daw > 0.34 and 1.e-15 for daw < 0.26
C
C   The range of validity for the model is 0.26 < daw < 0.34.  
C   For values outside of this range xjhom is set to a very high 
C   value (1.e19) or very low value (1.e-15). These values have no 
C   physical significance, but will normaly yield acceptable values 
C   for use in atmospheric models. 
C
C--------------------
C   calculates the water activity along the ice melting point 
C   curve, awice
C   (Eq. 2 in paper)
        dmu = 210368 +131.438*t -3.32373e+06/t -41729.1*log(t)
C   (Eq. 1 in paper)
        awice=exp(dmu/8.31441/t)
C--------------------
C    calculates the difference between awice and awliq 
C    (Eq. 6 in paper)
        daw=awliq-awice
C--------------------
C    calculates the ice nucleation rate as a function of daw 
C    (eq.7 in paper)
        dum2= -906.688 +8502.28*daw -26924.4*daw**2 +29179.6*daw**3
C--------------------
C    output
C
        IF(T.LT.170.0 .OR. T.GT.240.0)THEN
              xjhom=-99.e0  !OUTSIDE OF VALIDITY RANGE!
        ELSEIF(daw.GT.0.34)THEN
              xjhom=1.e19   !OUTSIDE OF VALIDITY RANGE!
        ELSEIF(daw.LT.0.26)THEN
              xjhom=1.e-15  !OUTSIDE OF VALIDITY RANGE!
        ELSE
              xjhom=10**(dum2)
        ENDIF
C
      END
C
C---------------------------------------
