Subroutine MoistAirProperties(mode,wbmd,psydat,emode,status)
!**********************************************************************************************
! Psychrometric chart subroutine
!
! modes 1-5:
! These modes take as input patm (in atmospheres), a dry bulb temperature, and one other property:
! wet bulb temperature, relative humidity (fraction), dew point temperature, humidity ratio, or enthalpy, 
! depending on mode. Outputs are humidity ratio (or relative humidity in mode 4), wet bulb temperature, enthalpy
! (or relative humidity in mode 5), mixture density and dry air density.
!
! mode 6:
! Takes as input patm, humidity ratio, and enthalpy, and returns all of the other properties - 
! setting h to hsat at the given w if a saturation condition is encountered
!
! mode 7:
! Takes as input patm, humidity ratio, and enthalpy, and returns all of the other properties - 
! setting w to wsat at the given h if a saturation condition is encountered
!
! mode 8:
! Takes as input patm, relative humidity, and enthalpy, and returns all of the other properties 
!
! The wet bulb temperature is only calculated in wbmode = 1.
!
! Temperatures are in Celsius (iunits=1) or Fahrenheit (iunits=2).
! Enthalpy is in kJ/kg (iunits=1) or BTU/lbm (inuits=2).
! Density is in kg/m^3 (iunits=1) or lbm/ft^3 (iunits=2).
! The reference states for enthalpies are hair=0.0 at 0 deg C and 0 def F and hw(liquid)=0.0 at 0 deg C and 32 deg F
!
! The psydat array contains the moist air properties.
!
!  Last modified:
!  2004-08-11 - JWT - Fixed a bug in MODE 7 where the humidity ratio was being set above the 
!                     inlet humidity ratio for very low temperatures.
!  2005-03-09 - JWT - Small bug with message 118 printed twice to LST, Bug in calculating 
!                     TSLOPE near 0 deg C
!  2006-04-20 - DAA - Fixed bug in printing MODE 2, as IUNIT and ITYPE were being set to zero.
!  2006-05-02 - DEB - Fixed a bug in PSYERR that was causing the UNIT and TYPE number calling
!                     PSYCHROMETRICS to be lost.
!  2009-06-00 - TPM - created MoistAirProperties for version 17 which uses CurUnit and CurType rather than info.
!
!**********************************************************************************************
! Copyright © 2009 Solar Energy Laboratory, University of Wisconsin-Madison. All rights reserved.

!    integer :: mode, wbmd, status, i
!    double precision, dimension(9) :: psydat
!     mode = 4
!     wbmd = 1
!     psydat(1) = 1
!     psydat(2) = 37.78
!     psydat(4) = .45
!     psydat(6) = 0.00523
!
!     
!    call MoistAirProperties(mode,wbmd,psydat,2,status)
!    do i = 1,9
!        write(*,10) psydat(i)
!        10 format(f10.4)
!    end do

Implicit None
    INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
    INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
    INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision

Integer wbmd,wbmode,emode,status,mode,itest,isign_old,isign_new,i_check,iunits
real(DP) psydat(9),ra,paold,patm,tdb,twb,twbold,rh,tdp,w,h,hdp,hmin,psatdb,psat,p,wsat,rhmax,wmax,hmax,balance, &
   dpress,y,tboil,twbnew,wsstar,wnew,err,tslope,pslope,wsslp,wslope,errslp,derrdt,pv,spcvol,rhowa,rhowm,delta_t,h_new,w_temp
save
Data ra/287.055d0/,paold/0.d0/

! Assign local variables:
 wbmode = wbmd
 patm = psydat(1)
 tdb =  psydat(2)
 twb =  psydat(3)
 twbold = -9.999d20
 rh =   psydat(4)
 tdp =  psydat(5)
 w =    psydat(6)
 h =    psydat(7)
 status = 0

! Unit conversions 
 If (iunits == 2) Then
    tdb = (tdb - 32.d0)/1.8d0
    twb = (twb - 32.d0)/1.8d0
    tdp = (tdp - 32.d0)/1.8d0
    h   = (h - 7.687d0)/0.43002d0
 EndIf

! Check that the total pressure is within the ideal gas range
 If (patm > 5.d0) Call PsyErr (emode,1,status)
 If (patm <= 0.d0) Then
    write(*,20) "total pressure is negative"
    20 format(a)
    Return
 EndIf

! Set w to zero if the value is very small to avoid problems with the correlations
 If (w < 1.d-08) w=0.

! For mode 6, check that the enthalpy is greater than the saturation enthalpy (minimum) for the given
! humidity ratio and that the humidity ratio is greater than o. 
 If (mode ==  6) Then
    If (w < 0.d0) Then
       Call PsyErr(emode,3,status)
       w = 0.d0
       goto 99
    ElseIf (w == 0.d0) Then
       goto 99
    EndIf
    Call DewPt_ASH97(patm,w,tdp,status,emode)
    hdp = tdp + w*(2501.d0 + 1.805d0*tdp)
    hmin = hdp
    If (h < hmin) Then
       Call PsyErr(emode,2,status)
       h = hmin
    EndIf
99  tdp = (h - 2501.d0*w)/(1.d0 + 1.805d0*w)
 EndIf

! For mode 7, check that the humidity ratio is less than the saturation humidity ratio for the given enthalpy
! and that the humidity ratio is greater than 0.
 If (mode == 7) Then
! Set a temporary w as sometimes with very low temps w can be calculated poorly
    w_temp = w
    If (w < 0.d0) Then
       Call PsyErr(emode,3,status)
       w = 0.d0
       goto 991
    ElseIf (w == 0.d0) Then
       goto 991
    EndIf
! Get the temperature and enthalpy at the dew point
    Call DewPt_ASH97(patm,w,tdp,status,emode)
    hdp = tdp + w*(2501.d0 + 1.805d0*tdp)
! Check for enthalpy less than saturation given the humidity ratio
    If (h < hdp) Then
! Start the iterative process by guessing a dry bulb temperature
       Call DewPt_ASH97(patm,w,tdp,status,emode)
	   delta_t=10.d0
	   isign_old=0
! Get the saturation pressure of the water at the current dry-bulb temperature
992    Call Sat_ASH97(tdb,psatdb,status,emode)
! Set the humidity ratio knowing rh=1
       rh = 1.d0
       w = .62198d0*psatdb*rh/(patm-psat*rh)
! Set the enthalpy of the air at this quess
       h_new = tdb + w*(2501.d0 + 1.805d0*tdb)
! Check the enthalpy balance
	   balance = DABS(h_new-h)
	   If (balance <= 0.001d0) goto 991
! Check energy balance
	   If (h_new > h) Then
	      isign_new=1
	   Else
	      isign_new=-1
	   EndIf
	   i_check=isign_new+isign_old
! Adjust the air temperature quess
	   If (i_check == 1) Then
	      tdb=tdb-delta_t
	      isign_old=isign_new
	   ElseIf (i_check == -1) Then
	      tdb=tdb+delta_t
	      isign_old=isign_new
	   ElseIf (i_check == -2) Then
	      tdb=tdb+delta_t
	      isign_old=isign_new
	   ElseIf (i_check == 2) Then
	      tdb=tdb-delta_t
	      isign_old=isign_new
	   Else
	      isign_old=isign_new
	      If (isign_new == 1) Then
	         delta_t=delta_t/10.d0
	         If (delta_t < 0.001) goto 991
	         tdb=tdb-delta_t
	      Else
	         delta_t=delta_t/10.d0
	         If (delta_t < 0.001) goto 991
	         tdb=tdb+delta_t
	      EndIf
	   EndIf
	   goto 992
    EndIf
! Reset the dry bulb temperature from our last guess of w and h
991 tdb = (h - 2501.d0*w)/(1.d0 + 1.805d0*w)
! Check w
    w=DMIN1(w,w_temp)
 EndIf

! For mode 8, check the rh and set the humidity ratio
 If (mode == 8) Then
    If (rh < 0.d0) Then
       Call PsyErr(emode,7,status)
	   rh = 0.d0
       w = 0.d0
       goto 995
    ElseIf (rh > 1.d0) Then
       Call PsyErr(emode,8,status)
	   rh = 1.d0
    EndIf

! Start the iterative process by guessing a dry bulb temperature
    tdb=H
	delta_t=10.d0
	isign_old=0
! Get the saturation pressure of the water at the current dry-bulb temperature
996 Call Sat_ASH97(tdb,psatdb,status,emode)
! Set the humidity ratio knowing the rh
    w = .62198d0*psatdb*rh/(patm-psat*rh)
! Set the enthalpy of the air at this quess
    h_new = tdb + w*(2501.d0 + 1.805d0*tdb)
! Check the enthalpy balance
	balance = DABS(h_new-h)
	If (balance <= 0.001d0) goto 995
! Check energy balance
	If (h_new > h) Then
	   isign_new=1
	Else
	   isign_new=-1
	EndIf
	i_check=isign_new+isign_old
! Adjust the air temperature quess
	If (i_check == 1) Then
	   tdb=tdb-delta_t
	   isign_old=isign_new
	ElseIf (i_check == -1) Then
	   tdb=tdb+delta_t
	   isign_old=isign_new
	ElseIf (i_check == -2) Then
	   tdb=tdb+delta_t
	   isign_old=isign_new
	ElseIf (i_check == 2) Then
	   tdb=tdb-delta_t
	   isign_old=isign_new
	Else
	   isign_old=isign_new
	   If (isign_new == 1) Then
	      delta_t=delta_t/10.d0
	      If (delta_t < 0.001) goto 995
	      tdb=tdb-delta_t
	   Else
	      delta_t=delta_t/10.d0
	      If (delta_t < 0.001) goto 995
	      tdb=tdb+delta_t
	   EndIf
	EndIf
	goto 996
995 tdb = (h - 2501.d0*w)/(1.d0 + 1.805d0*w)
 EndIf
! Find saturation pressure of water at wet bulb, dry bulb, or dew point temperature
 Call Sat_ASH97(tdb,psatdb,status,emode)
 Select Case (mode)
 Case (1)
! Check for impossible wet bulb temperatures and correct them if possible.
    If (twb > tdb .and. psatdb >= patm) Then
       write(*,20) "Twb > Tdb and psat >= patm"
       Return
    ElseIf (twb > tdb) Then
       Call PsyErr(emode,4,status)
       twb = tdb
       psat = psatdb
    Else
       Call Sat_ASH97(twb,psat,status,emode)
! Error: if psatwb is greater than patm
       If (psat >= patm) Then
          write(*,20) "psat >= patm"
          Return
       EndIf
    EndIf
 Case (2,4,5,6,7,8)
    psat = psatdb
 Case (3)
! Check for impossible dew point temperatures and correct them if possible
    If (tdp > tdb .and. psatdb >= patm) Then
       write(*,20) "Tdp > Tdb and psat >= patm"
       Return
    ElseIf (tdp > tdb) Then
       Call PsyErr(emode,5,status)
       tdp = tdb
    EndIf
    Call Sat_ASH97(tdp,psat,status,emode)
! Error: if psatdp is greater than patm
    If (psat >= patm) Then
       write(*,20) "psat >= patm"
       Return
    EndIf
 End Select

! mode 1 -- dry bulb and wet bulb supplied
 If (mode == 1) Then
    If (twb <= 0.d0) Then
       p = psat - 5.704d-4*(tdb-twb)*patm
       w = .62198d0 * p/(patm-p)
    Else
       wsat = .62198d0 * psat/(patm-psat)
       w = wsat - (tdb-twb)*(0.24D0 + .441D0*wsat)/(597.31d0 + 0.441D0*tdb - twb)
    EndIf
    If (w < 0.d0) Then
       Call PsyErr(emode,6,status)
       w = 0.0
       h = tdb
       twbold = twb
       wbmode = 1
       psat = psatdb
    Else
! According to ASHRAE 97
	   h = 1.006d0*tdb + w*(2501.d0 + 1.805d0*tdb)
       goto 100
    EndIf
 EndIf

! mode 2 -- dry bulb and relative humidity supplied
 If (mode == 2) Then 
    If (rh < 0.d0) Then
       Call PsyErr(emode,7,status)
       rh = 0.d0
    ElseIf (psat >= patm) Then
       rhmax = patm/psat
       If (rh >= (.99d0*rhmax)) Then
	      write(*,20) "fog"
          !do not stop if fog occurs
          rh = 0.99d0
       EndIf
    ElseIf (rh > 1.d0) Then
       rh = 1.d0
    EndIf
    w = .62198d0 * psat*rh/(patm-psat*rh)
 EndIf

! mode 3 -- dry bulb and dew point supplied  
 If (mode == 3) w = .62198d0 * psat/(patm-psat)

! Find enthalpy for modes 2-4
 If ((mode==2).or.(mode==3).or.(mode==4)) Then 
    If (psatdb < patm) Then
       wmax = .62198d0 * psatdb/(patm-psatdb)
       If (w > wmax) Then
          Call PsyErr(emode,9,status)
          w = wmax
       EndIf
    EndIf
    If (w < 0.d0) Then
       Call PsyErr(emode,3,status)
       w = 0.d0
    EndIf
    h = tdb + w*(2501.d0 + 1.805d0*tdb)
 EndIf

! mode 5 -- dry bulb and enthalpy supplied 
 If (mode==5) Then
    If (psatdb < patm) Then
       wmax = .62198d0 * psat/(patm-psat)
       hmax = tdb + wmax*(2501.d0 + 1.805d0*tdb)
       If (h > hmax) Then
          Call PsyErr(emode,10,status)
          h = hmax
       EndIf
    EndIf
    hmin = tdb
    If (h < hmin) Then
       Call PsyErr(emode,11,status)
       h = hmin
    EndIf
    w = (h-tdb)/(2501.d0+1.805d0*tdb)
 EndIf

! Find wet bulb temperatures if wbmode equals 1 
 If (wbmode /= 1) Then
    twb = tdb
    goto 100
 EndIf
 dpress = DABS(1.d0-patm)

! The following correlation is for 1 atmospheres total pressure. If outside the correlation range, the 
! correlation is used for the initial guess in the iterative method.
 If (h > 0.d0 .and. h < 2000.d0) Then
    y = DLOG(h*.43002d0+7.687d0)
    If (h <= 9.473d0) Then
       twb=-17.4422d0+1.9356d0*y+.7556d0*y**2+.5406d0*y**3
    ElseIf (h > 9.473d0) Then
       twb=-.6008d0-22.04556d0*y+11.4356d0*y**2-.97667d0*y**3
    EndIf
 Else
    twb = 9.99999d25
 EndIf

! Use a Newton's iterative method to find the wet bulb temperature
 If (dpress>=0.d0 .or. h<=0.d0 .or. h>275.d0) Then
    itest=0
    If (DABS(patm-paold) > 1.0d-10) Then
       Call Boil(patm,tboil)
       paold = patm
    EndIf
! Inital guess
    twbnew = DMIN1(twb,(tboil-0.1d0),tdb)
70  If (twbnew >= (tboil-0.09d0)) twbnew=tboil-0.1d0
    Call Sat_ASH97(twbnew,psat,status,emode)
    wsstar=.62198d0*psat/(patm-psat)
    If (mode == 5) w = (h-tdb)/(2501.d0+1.805d0*tdb)
    wnew=((2501.d0-2.381d0*twbnew)*wsstar-(tdb-twbnew))/(2501.d0+1.805d0*tdb-4.186d0*twbnew)
    err = w - wnew
    If (DABS(err) <= (.001d0*w)) goto 75
    If (w == 0.d0) Then
       If (DABS(err) <= .0001d0) goto 75
    EndIf
    itest = itest + 1
    If (itest >= 25) goto 75
! Find the slope of the error function
    tslope = 0.9999d0*(twbnew+273.15)-273.15
    If (twbnew == 0.d0) tslope = -.005d0
    Call Sat_ASH97(tslope,pslope,status,emode)
    wsslp=.62198d0*pslope/(patm-pslope)
    wslope=((2501.d0-2.381d0*tslope)*wsslp-(tdb-tslope))/(2501.d0+1.805d0*tdb-4.186d0*tslope)
    errslp = w - wslope
	If (ABS(errslp)<=0.0001d0) goto 75  !exit iteration if the slope of the error is small
    derrdt = (errslp - err)/(tslope - twbnew)
    twbnew = twbnew - err/derrdt
    goto 70
75  twb = twbnew
 EndIf
 If (twb < twbold) twb = twbold

100 Continue

! Find relative humidity, dew point, mixture density, and dry air density 
 pv  = patm*w/(.62198d0+w)
 If (mode /= 2) rh = pv/psatdb
 If ((mode == 7) .and. (rh > 1.d0)) Then
    rh =1.d0
	w = 0.62198d0*psatdb/(patm-psatdb)
 EndIf
 If (mode /= 3 .and. pv > 0.d0) Then
    Call DewPt_ASH97(patm,w,tdp,status,emode)
 ElseIf (mode /= 3 .and. pv <= 0) Then
! For dry air there is no dew point temperature
    tdp=-9.99999d25
 EndIf
 spcvol = ra*(tdb+273.15d0)/(patm*101325d0)*(1.d0+1.6078d0*w)
 rhowa = 1.d0/spcvol
 rhowm = rhowa*(1.d0+w)

! Convert outputs to appropriate units 
 If (iunits == 2) Then
    h = h*0.43002d0+7.68d0
    tdb = 1.8d0*tdb + 32.d0
    twb = 1.8d0*twb + 32.d0
    If (tdp > -9.99999d24) tdp = 1.8d0*tdp + 32.d0
    rhowm = rhowm/16.02d0
    rhowa = rhowa/16.02d0
 EndIf

! Set outputs
 If ((mode == 6).or.(mode == 7).or.(mode == 8)) psydat(2) = tdb
 psydat(3) = twb
 psydat(4) = DMIN1(1.,(DMAX1(0.,rh)))
 psydat(5) = tdp
 psydat(6) = w
 psydat(7) = h
 psydat(8) = rhowm
 psydat(9) = rhowa
 
 Return
End

Subroutine Sat_ASH97(tin,psat,status,emode)
!***************************************************************************
! Subroutine for finding saturation pressure or water at a given temperature
!***************************************************************************
Implicit none
Integer emode,status
Double Precision tin,psat,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,t
! The following correlation for the saturation pressure of water vapor (in Pascals) as a function of
! temperature is taken from the 1997 ASHRAE Fundamentals Handbook (SI).
Data c1/-5.674359d3/,c2/6.3925247d0/,c3/-9.677843d-3/,c4/6.2215701d-7/,c5/2.0747825d-9/,c6/-9.484024d-13/
Data c7/4.1635019d0/,c8/-5.8002206d3/,c9/1.3914993d0/,c10/-4.8640239d-2/,c11/4.1764768d-5/,c12/-1.4452093d-8/,c13/6.5459673d0/
 
 t = tin + 273.15d0
 If (t <= 0.d0) Then
    write(*,20) "Fatal error in saturation subroutine"
    20 format(a)
    Return
 EndIf
! Saturation pressure over ice (-100 to 0 C) 
 If (t < 273.15d0) Then
    psat=DEXP(c1/t+c2+c3*t+c4*t**2+c5*t**3+c6*t**4+c7*DLOG(t))
! Saturation pressure over liquid water (0 to 200 C)
 ElseIf (t >= 273.15d0) Then
    psat=DEXP(c8/t+c9+c10*t+c11*t**2+c12*t**3+c13*DLOG(t))
 EndIf
! Temperature out of the range used for the correlation 
 If (t < 173.15d0 .or. t > 473.15d0) Then
     Call PsyErr(emode,12,status)
 EndIf
! Convert pressure from Pascals to atmospheres
 psat = psat/101325.d0
 Return
End 

Subroutine DewPt_ASH97(patm,w,tdp,status,emode)
!***************************************************************************
! Subroutine for finding the dew point temperature given the humidity ratio. The correlation is from the
! 1997 ASHRAE Fundamentals Handbook. The dew point temperature is in degrees C and atmospheric pressure in kilopascals.
!***************************************************************************
Implicit none
Integer emode,status
Double Precision patm,w,tdp,pv,y     
 
 pv  = 101.3d0*patm*w/(.62198d0+w)
! If pv is negative, you're going to get a math error below. If this becomes a problem, we may want to change it...
! however, check if the restriction is correct first.
 y = DLOG(pv) !y=ALOG(MAX(pv,0.))
 tdp =6.54d0+14.526d0*y+0.7389d0*y*y+0.09486d0*y*y*y+0.4569d0*pv**0.1984d0
 If (tdp < 0.d0) tdp = 6.09d0 + 12.608d0*y + 0.4959d0*y*y
 If (tdp > 93.d0 .or. tdp < -60.d0) Call PsyErr(emode,13,status)
 Return
End

Subroutine Boil(patm,tboil)
!****************************************************************************
! Subroutine for finding the boiling temperature of water given the total pressure. A Newton's method is used with the
! saturated water vapor pressure correlation from the 1985 ASHRAE Fundamentals Handbook.
!****************************************************************************
Implicit none
Integer itest
Double Precision patm,tboil,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,pboil,t1,zz,p1,err,derrdt,told

Data c1/-5674.5359d0/,c2/6.3925247d0/,c3/-0.9677843d-2/,c4/0.62215701d-6/,c5/0.20747825d-8/,c6/-0.9484024d-12/
Data c7/4.1635019d0/,c8/-5800.2206d0/,c9/1.3914993d0/,c10/-0.048640239d0/,c11/0.41764768d-4/,c12/-0.14452093d-7/,c13/6.5459673d0/
 
 pboil = patm*101325.d0
 itest = 0
 If (pboil >= 611.21d0) Then
! Using ASHRAE correlation for 0 to 200 C. First guess tboil = 100 C.
    t1 = 100.d0 + 273.15d0
10  zz = c8/t1+c9+c10*t1+c11*t1**2+c12*t1**3+c13*DLOG(t1)
    p1 = DEXP(zz)
    err = pboil - p1
    If (DABS(err) <= (.01d0*pboil)) Then
       tboil = t1 - 273.15d0
       Return
    EndIf
    itest = itest + 1
    If (itest > 100) Then
       write(*,20) "Fatal error in boiling subroutine"
       20 format(a)
       Return
    EndIf
    derrdt = -(p1)*(-c8/t1**2+c10+2*c11*t1+3*c12*t1**2+c13/t1)
    told = t1
    t1 = t1 - err/derrdt
    goto 10
 Else
! Using ASHRAE correlation for -100 to 0 C. First guess tboil = 0 C.
    t1 = 273.15d0
110 zz = c1/t1+c2+c3*t1+c4*t1**2+c5*t1**3+c6*t1**4+c7*DLOG(t1)
    p1 = DEXP(zz)
    err = pboil - p1
    If (ABS(err) <= (.01d0*pboil)) Then
       tboil = t1 - 273.15d0
       Return
    EndIf
    itest = itest +1
    derrdt = -p1*(-c1/t1**2+c3+2*c4*t1+3*c5*t1**2+4*c6*t1**3+c7/t1)
    told = t1
    t1 = t1 - err/derrdt
    goto 110
 EndIf
End
      
Subroutine PsyErr(emode,enmb,status)
!***************************************************************************
! Subroutine for handling warning messages
!
! emode is 0, 1, OR 2.  
!      0 - Do not print warning messages, return warning number (status)
!      1 - Print only one warning for the entire simulation
!      2 - Print one warning each timestep
! Warnings are printed for each warning condition and each unit (up to iumax) calling psych
!***************************************************************************
      
Double Precision Time,wtime,Timestep
Integer enmb,status,emode,line
Character (len=12) IUMAXStr
Character (len=300) PsychErr
 
! iumax is the maximum number of units numbers that will be stored. If more than iumax units call psych with the
! same warning condition then the unit numbers are no longer used in determining if warnings have already been printed.
! imax is the number of possible warning conditions (enmb must be between 1 and iwmax)

! Check for allowable emode (0,1,2) 
 If (emode < 0 .or. emode > 2) Then
    write(*,20) "Fatal error reported by psychrometrics; emode is out of range"
    20 format(a)
    Return 
 EndIf
! Set status 
 status = enmb
! Transfer according to emode
 line = emode + 1
 If (emode == 0) Return
 
! Warning conditions
 Select Case (enmb) 
 Case (1) 
    PsychErr = TRIM(PsychErr)//'a total pressure greater than 5 atmospheres. Beware that the idea gas relations used in the psychrometrics routine are not accurate at such pressures.'
 Case (2) 
 	PsychErr = TRIM(PsychErr)//'an enthalpy less than possible for the given humidity ratio. The air was assumed to be saturated at the given humidity ratio and the enthalpy was set for this condition.'
 Case (3) 
 	PsychErr = TRIM(PsychErr)//'a humidity ratio less than 0.0. Dry air was assumed and the humidity ratio was set to 0.0.'
 Case (4) 
 	PsychErr = TRIM(PsychErr)//'a wet bulb temperature greater than the drybulb temperature.  Saturated air was assumed and the wet bulb temperature set to the dry bulb temperature.'
 Case (5) 
 	PsychErr = TRIM(PsychErr)//'a dew point temperature greater than the dry bulb temperature.  The dew point was set equal to the dry bulb temperature.'
 Case (6) 
 	PsychErr = TRIM(PsychErr)//'a wet bulb temperature below that for dry air at the given dry bulb temperature. Dry air was assumed and a new wet bulb temperature found.'
 Case (7) 
 	PsychErr = TRIM(PsychErr)//'a relative humidity less than 0.0. The relative humidity was set to 0.0.'
 Case (8) 
 	PsychErr = TRIM(PsychErr)//'a relative humidity greater than 1.0. The relative humidity was set to 1.0.'
 Case (9) 
 	PsychErr = TRIM(PsychErr)//'a humidity ratio above the saturation humidity ratio for the given dry bulb temperature. The air was assumed saturated at the given dry bulb temperature and the humidity ratio set for this condition.'
 Case (10) 
 	PsychErr = TRIM(PsychErr)//'an enthalpy greater than that for saturated air at the given dry bulb temperature. The air was assumed saturated and the enthalpy set for this condition.'
 Case (11) 
  	PsychErr = TRIM(PsychErr)//'an enthalpy less then that for dry air at the given dry bulb temperature. Dry air was assumed and the enthalpy set for this condition.'
 Case (12) 
 	PsychErr = TRIM(PsychErr)//'conditions outside the intended range of the water vapor saturation pressure correlation. The correlation, from the 1985 ASHRAE fundamentals handbook, is for temperature range of -100 C to 200 C.'
 Case (13) 
 	PsychErr = TRIM(PsychErr)//'conditions outside the intended range of the dew point temperature correlation. The correlation from the 1981 ASHRAE fundamentals handbook is for the temperature range of -60 C to 70 C.'
 Case (14) 
 	PsychErr = TRIM(PsychErr)//'a humidity ratio less than possible for the given enthalpy. The air was assumed to be saturated at the given enthalpy and the humidity ratio set for this condition.'
 End Select

 Return
End 





