C***********************************************************************
C  FORTRAN77 CODE TO CALCULATE THE VAPOR PRESSURE 
C  OF PURE LIQUID WATER FROM 123 K TO 332 K
C  INCLUDING THE SUPERCOOLED TEMPERATURE RANGE
C
C  BY DAN MURPHY <daniel.m.murphy@noaa.gov>
C  AND 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  MURPHY, D.M., AND KOOP, T.,
C  Review of the vapour pressures of ice and supercooled
C  water for atmospheric applications,
C  Q. J. R. Meteorol. Soc. 131, 1539-1565, 2005.  
C***********************************************************************
C
C  INPUT: 
C  t =Temperature in [Kelvin] 
C  
C  OUTPUT:
C  ph2o0  = vapor pressure over liquid water in [Pa]
C  
C***********************************************************************
C
      PROGRAM p_h2o0
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C***********************************************************************
C
C   MAIN PROGRAM
C
      OPEN(11,FILE='temp.dat',STATUS='unknown')
 98   FORMAT(1x,F10.6,2(E15.7))
C
C   INPUT
C
      tstart=332.0d0
      tend=10.d0
      tstep=-1.0d0 
C
      DO t=tstart,tend,tstep
C
        xxx=ph2o0(t)   
C      
C   OUTPUT
C
        write(11,98)t,xxx,1/t
C
      ENDDO
      close(11)
      END
C
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  LIQUID WATER
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C---------------------------------------
C    VAPOR PRESSURE OF LIQUID WATER 
C---------------------------------------
C
      FUNCTION ph2o0(t)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     Returns -99 for T =< 123K and T >= 332 K.
C
        dum=54.842763d0-6763.22d0/T-4.210d0*dlog(T)+0.000367d0*T
     >    +dtanh(0.0415d0*(T-218.8d0))
     >    *(53.878d0-1331.22d0/T-9.44523d0*dlog(T)+0.014025d0*T)
C
      IF(T.GT.123.d0 .AND. T.LT.332.d0)THEN
        ph2o0=dexp(dum)
      ELSE
        ph2o0=-99.d0
      ENDIF
C
      END
C
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  HEXAGONAL ICE Ih  
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C---------------------------------------
C    VAPOR PRESSURE OF HEXAGONAL ICE Ih
C---------------------------------------
C
      FUNCTION piceIh(t)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     Returns -99 for T =< 110K and T > 273.16 K.
C
      dum=9.550426d0-5723.265d0/T+3.53068d0*dlog(T)-0.00728332d0*T
      IF(T.GT.110.d0 .AND. T.LE.273.16d0)THEN
        piceIh=dexp(dum)
      ELSE
        piceIh=-99.d0
      ENDIF
C
      END
C
C
C---------------------------------------
C    ATMOSPHERIC FROST POINT TEMPERATURE
C    OF HEXAGONAL ICE
C---------------------------------------
C       
      FUNCTION tfrost(pw)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     Returns -99 for Tfrost =< 115K and Tfrost > 273.16K.
C
        dum=(1.814625d0*dlog(pw)+6190.134d0)/(29.120d0-dlog(pw))
C
      IF(dum.GT.115.d0 .AND. dum.LE.273.16d0)THEN
         tfrost=dum
      ELSE
         tfrost=-99.d0
      ENDIF
C
      END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
