C***********************************************************************
C  FORTRAN77 CODE TO CALCULATE THE ATMOSPHERIC FROST POINT
C  TEMPERATURE WITH RESPECT TO HEXAGONAL ICE
C  FROM 115 K TO 273.16 K
C
C  BY DAN MURPHY <daniel.m.murphy@noaa.gov>
C  AND THOMAS KOOP <thomas.koop@uni-bielefeld.de>
C  Version of March 7, 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  pw = atmospheric water partial pressure in [Pa]
C
C  OUTPUT:
C  tfrost  = atmospheric frost point temperature in [K]
C
C***********************************************************************
C
      PROGRAM tfrost
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C***********************************************************************
C
C   MAIN PROGRAM
C
      OPEN(11,FILE='temp.dat',STATUS='unknown')
 98   FORMAT(1x,E15.7,F11.5,E15.7)
C
C   INPUT
C
      pwstart=3.d0
      pwend=-10.d0
      pwstep=-1.0d0
C
      DO dummy=pwstart,pwend,pwstep
      DO dummy2=1.0d0,0.11d0,-0.1d0

        pw=dummy2*10**dummy
C
        xxx=tfrost(pw)
C
C   OUTPUT
C
        write(11,98)pw,xxx
C
      ENDDO
      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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
