      SUBROUTINE PSY(IPSY,TDB,TWB,TDP,RH,W,H)
      LOGICAL SI
C
C  THIS ROUTINE EVALUATES PSYCHOMETRIC CONDITIONS FOR AIR AT 1 ATM
C  REFER TO ASHRAE 1977 FUNDAMENTALS, CHAPTER 5
C   BY F.E.JAKOB, BATTELLE (ABOUT 1984)
C  IPSY IS A KEY INDICATOR TO THE PROGRAM.
C  GIVEN ANY TWO PSYCHOMETRIC VARIABLES, THE OTHER 4 ARE EVALUATED.
C  THE TABLE BELOW SPECIFIES THE VALUE OF IPSY REQUIRED FOR EACH
C  PAIR COMBINATION.  NOTE THAT TDB IS ALWAYS REQUIRED.
C
C  IF(IPSY.GT.ZERO) THEN THE INCOMING DATA IS ASSUMED TO BE IN SI UNITS
C                   AND THEY ARE CONVERTED IN ORDER TO BE COMPATIBLE
C                   WITH THE ASHRAE EQUATIONS PROGRAMMED HERE.
C                   THE VARIABLES ARE CONVERTED BACK TO SI ON OUTPUT.
C
C  IF(IPSY.LT.ZERO) THEN THE INCOMING DATA IS ASSUMED TO BE
C                   IN ENGLISH UNITS. NO CONVERSIONS ARE MADE
C
C     IPSY   DESCRIPTION
C     ****   ************************
C      +-1    TDB AND TWB SPECIFIED
C      +-2    TDB AND TDP SPECIFIED
C      +-3    TDB AND RH  SPECIFIED
C      +-4    TDB AND  W  SPECIFIED
C      +-5    TDB AND  H  SPECIFIED
C
C  UNITS : 
C        VARIABLE  UNITS
C        ********  **********************************************
C        ********  (IPSY.LT.ZERO)************(IPSY.GT.ZERO)******
C        ********  **********************************************
C          TDB     DEG F                     DEG C
C          TWB     DEG F                     DEG C
C          TDP     DEG F                     DEG C
C           RH     PERCENT (%)               PERCENT (%)
C            W     (LBM H20)/(LBM AIR)       (KG H20)/(KG AIR)
C            H     (BTU)/(LBM OF DRY AIR)    KJ/(KG OF DRY AIR)
C
C        PRESSURE  ATMOSPHERES  (29.92 INCHES OF HG PER ATM)
C
C
C
C  RELATIVE HUMIDITY (EQ 25): 
      FRH(W,WS,PWS)=(W/WS)/(1.0-(1.0-W/WS)*PWS/P)*100.
C
C  PARTIAL VAPOR PRESSURE (NOT SATURATION VP) (EQ 36): 
      FPW(W)=P*W/(0.62198+W)
C
C  HUMIDITY RATIO (EQ 32 SOLVED FOR W): 
      FW32(TDB,H)=(H-0.240*TDB)/(1061.0+0.444*TDB)
C
C  ENTHALPY (EQ 32): 
      FH(TDB,W)=0.240*TDB+W*(1061.0+0.444*TDB)
C
C  DEW POINT (PW REQUIRED IN INCHES OF HG) (EQ 37): 
      FTDP(PW)=((ALOG(PW)*1.8893)+30.5790)*ALOG(PW)+79.047
C
C  TEMPERATURE CONVERSIONS
      FTOC(T) = (T - 32.0) / 1.8
      CTOF(T) = T * 1.8 + 32.0
C
C
      P=1.0
      SI = .FALSE.
      IF(IPSY.GT.0) SI = .TRUE.
      IF(.NOT.SI) GO TO 2
      TDB = CTOF(TDB)
      IF( IABS(IPSY).EQ.1 ) TWB = CTOF(TWB)
      IF( IABS(IPSY).EQ.2 ) TDP = CTOF(TDP)
      IF( IABS(IPSY).EQ.5 ) H = H / 2.3278
    2 IF(IABS(IPSY).GE.1.AND.IABS(IPSY).LE.5) GO TO 5
      WRITE(60,1000) IPSY
 1000 FORMAT(1H0,'!!!!!! TROUBLE!, IPSY=',I5,' IN SUBROUTINE PSY')
      RETURN
    5 GO TO (10,20,30,40,50) IABS(IPSY)
C
   10 CONTINUE
      IF(TWB .GT. TDB)THEN
         WRITE(60,1010) TWB,TDB
         TWB=TDB-.001
         ENDIF
      PWSS=FPWS(TWB)
      WSS=FW22(PWSS)
      W=FW35(TDB,TWB,WSS)
      PWS=FPWS(TDB)
      WS=FW22(PWS)
      RH=FRH(W,WS,PWS)
      H=FH(TDB,W)
      PW=FPW(W)
      TDP=FTDP(PW*29.92)
      GO TO 100
C
   20 CONTINUE
      IF(TDP .GT. TDB)THEN
         WRITE(60,1020) TDP,TDB
         TDP=TDB-.001
         ENDIF
      PW=FPWS(TDP)
      W=FW22(PW)
      PWS=FPWS(TDB)
      WS=FW22(PWS)
      RH=FRH(W,WS,PWS)
      H=FH(TDB,W)
      TWB=FTWB(TDB,RH,W)
      GO TO 100
C
   30 CONTINUE
      IF( RH.LT.-0.001.OR. RH.GT.100.001)THEN
         IF(RH.LT.-0.001)THEN
           WRITE(60,1030)RH,TDB
           RH=0.+.001
           ENDIF
         IF(RH.GT.100.001)THEN
           WRITE(60,1040)RH,TDB
           RH=99.999
           ENDIF
         ENDIF
      PWS=FPWS(TDB)
      PW=PWS*RH/100.0
      W=FW22(PW)
      WS=FW22(PWS)
      H=FH(TDB,W)
      TDP=FTDP(PW*29.92)
      TWB=FTWB(TDB,RH,W)
      GO TO 100
C
   40 CONTINUE
      PWS=FPWS(TDB)
      PW=FPW(W)
      RH=PW/PWS*100.0
      IF(RH .LE. 100.001) GO TO 45
      RH=99.999
      GO TO 30
   45 CONTINUE
      WS=FW22(PWS)
      H=FH(TDB,W)
      TDP=FTDP(PW*29.92)
      TWB=FTWB(TDB,RH,W)
      GO TO 100
C
   50 CONTINUE
      PWS=FPWS(TDB)
      W=FW32(TDB,H)
      IF(W .LT. 0.00001) W=0.00001
      PW=FPW(W)
      RH=PW/PWS*100.0
      IF(RH .LE. 100.001) GO TO 55
      RH=99.999
      GO TO 30
   55 CONTINUE
      WS=FW22(PWS)
      TDP=FTDP(PW*29.92)
      TWB=FTWB(TDB,RH,W)
      GO TO 100
C
  100 CONTINUE
      IF(.NOT.SI) RETURN
      TDB = FTOC(TDB)
      TWB = FTOC(TWB)
      TDP = FTOC(TDP)
      H = H * 2.3271
      RETURN
 1010 FORMAT(1H ,'TROUBLE IN PSY: TWB = ',G13.3,' > TDB= ',G13.3)
 1020 FORMAT(1H ,'TROUBLE IN PSY: TDP = ',G13.3,' > TDB= ',G13.3)
 1030 FORMAT(1H ,'TROUBLE IN PSY: RH = ',G13.3,' < 0., TDB= ',G13.3)
 1040 FORMAT(1H ,'TROUBLE IN PSY: RH = ',G13.3,' > 100., TDB= ',G13.3)
      END
      FUNCTION FTWB(TDB,RH,W)
C
C  THIS ROUTINE ITERATIVLY SOLVES THE
C  WET BULB TEMPERATURE EQUATION
C
C     VARIABLE  UNITS
C     ********  ************************
C       TDB     DEG F
C        RH     PERCENT
C         W     (LBM H20)/LBM AIR)
C
      DATA LIMIT /20/ , XL100 /0.1/
C
C  EMPIRICAL APPROXIMATION FOR WET BULB TEMPERATURE
C  (USED FOR INITIAL GUESS)
      FTWBG(TDB,RH)=9.64*EXP(-0.05*RH)+
     + (0.15+0.18*ALOG(RH))*TDB
C
C  WET BULB TEMPERATURE (EQ 35 SOLVED FOR T*=TWB): 
      FTWBN(TDB,W,WSS)=(1093.0*(W-WSS)+
     + (0.444*W+0.240)*TDB)/
     + (W+0.240-0.556*WSS)
C
C  HUMIDITY RATIO (EQ 23): 
      FW(PWS)=0.62198*PWS/(P-PWS)
C-------------------------------------------------------------
C
      IF(ABS(RH - 100.).LE.0.01) GO TO 200
      IF(RH.LT.-0.001 .OR. RH.GT.100.001 )THEN
        IF(RH.LT.-.001)THEN
          RH=0.001
        ELSE IF (RH.GT.100.01)THEN
          RH=99.999
          ENDIF
        ENDIF
      P=1.0
      TWBO=TDB+5.0
      TWBNO=FTWBG(TWBO,RH)
      TWB=TDB
      TWBN=FTWBG(TWB,RH)
C
      DO 100 I100=1,LIMIT
      SLOPE=(TWBN-TWBNO)/(TWB-TWBO)
      XINT=TWBN-SLOPE*TWB
      TWBO=TWB
      TWBNO=TWBN
      TWB=XINT/(1.0-SLOPE)
      PWSS=FPWS(TWB)
      WSS=FW(PWSS)
      TWBN=FTWBN(TDB,W,WSS)
  100 IF(ABS(TWB-TWBN).LE.XL100) GO TO 110
C
      IF(I100.GE.LIMIT) WRITE(60,1000) I100,TDB,RH,W,
     + TWBO,TWBNO,TWB,TWBN,SLOPE,XINT
      STOP 'FTWB CONVERGANCE PROBLEM'
C
  110 FTWB=TWBN
      RETURN
  200 CONTINUE
      FTWB=TDB
      RETURN
 1000 FORMAT(1H0, I10, 3G10.4, (/1X,2G10.4) )
 1010 FORMAT(1H0,'TROUBLE IN FTWB',3G10.3)
      END
      FUNCTION FPWS(TF)
C
C  EVALUATE THE SATURATION PRESSURE OF
C  WATER VAPOR.
C
C  TF IS EXPECTED IN DEG F
C
C
C  KEENAN, KEYES, HILL, MOORE FORMULA (EQ 5)
C
      DIMENSION F(8)
      DATA F/-741.9242,-29.72100,-11.55286,
     + -0.8685635,0.1094098,0.439993,
     + 0.2520658,0.05218684/
C
      TC=(TF-32.0)/1.8
      TK=TC + 273.15
      SUM=F(1)
C
      DO 10 I10=2,8
   10 SUM=F(I10)*(0.65-0.01*TC)**(I10-1)+SUM
      FPWS=217.99*EXP(0.01/TK*(374.136-TC)*SUM)
      RETURN
      END
      FUNCTION FW22(PW)
C  HUMIDITY RATIO (EQ 22): 
      DATA P/1./
      FW22=0.62198*PW/(P-PW)
      RETURN
      END
      FUNCTION FW35(TDB,TWB,WSS)
C  HUMIDITY RATIO (EQ 35): 
      FW35=((1093.0-0.556*TWB)*WSS
     + -0.240*(TDB-TWB))/(1093.0+0.444*TDB-TWB)
      RETURN
      END
      
C
