C
C
      SUBROUTINE TCORR (KQ,ID,I)
C
C     -----------------------------------------------------------------
C
C       SUBROUTINE PURPOSES:
C       EVALUATION OF WALL TEMPEARTURE CORRECTION IN WET CASES.
C
C     -----------------------------------------------------------------
C
      INCLUDE 'CMMNW.DAT'
C
      PARAMETER (NT1=45,NT2=55)
      PARAMETER (KMX=NT1*NT2*2, KNX=NMZ*10)
C
      DOUBLE PRECISION BJ00,BJ01,BJ10,BY00,BY01,BY10,RO,R1,H,SKAPPA,
     .                 SKAPA1,STEP
C
      REAL M,LW,BETA
C
      COMMON
     . /BESSEL/ BJ00,BJ01,BJ10,BY00,BY01,BY10
     . /PARAMT/ RO,R1,H,STEP
C
C       NT1 DETERMINES HOW MANY DATA POINTS CAN BE HANDLED FOR TEMP.
C     CORRECTION CALCULATION DUE TO WATER EVAPORATION. NT2 DETERMINES
C     HOW MANY TERMS IN THE INFINITE SERIES OF THE SOLUTION ARE
C     ACTUALLY SUMMED UP. THE CONVERGENCE IS DOMINANTLY CONTROLLED
C     BY THE MAGNITUDE OF R1. ENLARGED R1 GENERALLY REQUIRES LARGER
C     NT2 WHICH IN TURN INCREASES COMPUTING TIME SIGNIFICANTLY. THE
C     CURRENTLY ASSUMED R1 IS GOOD ENOUGH FOR DYNAMIC SIMULATION
C     COVERING A FEW HOURS.
C
      DIMENSION  SKAPPA(NT1,2,NT2),YKAPPA(NT1,2,NT2),H11(NMZ,10),
     .           B1(NMZ,10),NLS(NMZ,10)
      DATA SKAPPA,YKAPPA,H11,B1,NLS
     .     / KMX*0.0,KMX*0.0,KNX*0.0,KNX*0.0,KNX*0 /
C
C     SKAPPA(NT1,2,NT2): EIGENVALUES.
C     YKAPPA(NT1,2,NT2): INTERMEDIATE STORAGE ARRAY.
C     H11(NMZ,10): INTERMEDIATE STORAGE ARRAY.
C     B1(NMZ,10): INTERMEDIATE STORAGE ARRAY.
C     NLS(NMZ,10): INTERVAL RECORD ARRAY.
C
      PI=3.1415926D0
      LL=0
C
C  ** DATA BASE SETUP.
C     DELTAT: TIME INCREMENT, SEC.   TRV: ROCK TEMP. DEGREE F.
C     RO: HYDROLIC RADIUS, FT.     HA: THERMAL DIFFUSIVITY, FT/HR.
C     HC: HEAT TRANSFER COEFF. BTU/(HR*FT**2*F).
C     HK: THERMAL CONDUCTIVEITY, BTU/(HR*FT*F).
C     R1: RADIUS OF INFLUENTIAL ZONE FOR TEMP. CORR. IN WET CASE.
C     NLS: NUMBER OF INTERVALS IN WHICH NO DETAILED CALCULATIONS
C          IN TEMP. CORR. WERE GIVEN. IT IS USED TO DETERMINE WHICH
C          FORMULA TO USE IN THE LATER INTERVALS.
C
      KN=NSEQE(KQ)+ID
      IF (KN.GT.NT1) THEN
         WRITE (8,100) NO(I)
         WTCORR(KQ,ID)=0.0
         WMTR(KQ,ID)=0.0
         BETA=0.0
         RETURN
      ENDIF
      TMSTT=DELTAT/3600.0
C
C     TROCK: ROCK TEMP.   TDA: AIR T. AT STATIONS.  STV: WALL TEMP.
C     HABS: HUMIDITY RATIO OF AIR AT STATIONS.
C
      TRV=TROCK(I)
      TA=TDA(KQ,ID)
      TO=STV(KQ,ID)
      AH=HABS(KQ,ID)
      LW=1050.0
      RO=2.0*A(I)/O(I)
      R1=RO+1.0
      STEP=1.0D0
C
C     STEP: LENGTH OF STEPS IN SEARCHING FOR EIGENVALUES.
C     GENERALLY, WHEN (R1-RO)<2, STEP<1.5, WHEN 2<(R1-RO)<10,
C     STEP<0.2, WHEN 10<(R1-RO)<30, STEP<0.05. FOR SAFE SEEK, TAKE
C     SMALLER VALUE OF "STEP".
C
      ALOGR=DLOG(R1/RO)
      HC=BI(I)*HK(I)*O(I)/(2.0*A(I))
      IF (TO.GT.210.0) TO=210.0
C
C     SUB. MTCOEF EVALUATES MASS TRANSFER COEFF. BETA.
C
      CALL MTCOEF (TO,TA,I,BETA)
      TV=TO-10.0
      TU=TO+10.0
      IF (TO.LE.122.0) THEN
         PVM=-1.0682+0.07946*TO-0.00129*TO*TO+1.26E-5*TO**3
         PVB=-1.0682+0.07946*TV-0.00129*TV*TV+1.26E-5*TV**3
         PVC=-1.0682+0.07946*TU-0.00129*TU*TU+1.26E-5*TU**3
      ELSE
         PVM=-33.786+0.97*TO-9.305E-3*TO*TO+3.65E-5*TO**3
         PVB=-33.786+0.97*TV-9.305E-3*TV*TV+3.65E-5*TV**3
         PVC=-33.786+0.97*TU-9.305E-3*TU*TU+3.65E-5*TU**3
      ENDIF
      GAM=0.622*PVM/(101.35-PVM)
      GAV=0.622*PVB/(101.35-PVB)
      GAU=0.622*PVC/(101.35-PVC)
      DRR=(DR+DR*AH*18/28.98)*(TR+460)/((TO+460)*(1+AH))
      SLOPE1=(GAM-GAV)*DRR/10.0
      SLOPE2=(GAU-GAM)*DRR/10.0
      GAA=AH*DRR
      GAW=GAM*DRR
C
C     KEQUIL: 1: THE WALL FILM AND BULK FLOW IS IN A TEMPORARY EQUIL.
C             CONDITION IN MASS TRANSFER. NO WATER EVAPORATION OR
C             CONDENSATION.
C
      KEQUIL=0
      IF (GAA.LT.GAW*PSI(I).AND.TO.LE.CRITW) THEN
         PSII=PSI(I)
         SLOPE=SLOPE1
         NLS(KQ,ID)=NLS(KQ,ID)+1
      ELSE IF (GAA.GT.GAW) THEN
         PSII=1.0
         SLOPE=SLOPE2
         NLS(KQ,ID)=NLS(KQ,ID)+1
      ELSE
         PSII=GAA/GAW
         SLOPE=SLOPE1
         NLS(KQ,ID)=NLS(KQ,ID)+1
         KEQUIL=1
      ENDIF
C
C     NLS: FLAG NUMBER CONDUCTING THE USE OF SOLUTION FORMULAS
C          FOR DIFFERENT TIME INTERVAL.
C
      ITNN=NLS(KQ,ID)
      XI=HC+LW*BETA*PSII*SLOPE
      ETA=LW*BETA*(PSII*GAW-GAA)
      M = ETA/XI
      B = M/TRV
      H = XI/HK(I)
      WN = 1.0+RO*H*ALOGR
      WN1=1.0+RO*H11(KQ,ID)*ALOGR
C
C  ** SOLUTION OF TEMP. CORR.
C
      TERM1 = B/(WN*(1.0+B))
      TERM2 = 0.0
      TERMM=0.0
      DO 40 L=1,NT2
C
C     SKAPA1: EIGENVALUE OF CHARACTERISTIC EQN.  SKAPPA: STORAGE
C     ARRAY OF EIGENVALUES.
C
            CALL  EIGEN (LL,SKAPA1)
            LL=1
            SKAPPA(KN,1,L) = SKAPA1
            Y1 = (SKAPPA(KN,1,L)*BJ10+H*BJ00)**2-
     .           (SKAPPA(KN,1,L)*SKAPPA(KN,1,L)+H*H)*BJ01*BJ01
            YR1 = Y1/(BJ01*BJ01)
            IF (ITNN.EQ.1) THEN
               YR = 1.0
            ELSE
               Y2 = B1(KQ,ID)*(H11(KQ,ID)-H)/(WN1*(1.0+B)*
     .              SKAPPA(KN,1,L)*SKAPPA(KN,1,L))
               Y2 = Y2+(B1(KQ,ID)-B)*H/((1.0+B)
     .              *SKAPPA(KN,1,L)*SKAPPA(KN,1,L))
               SUMYK = 0.0
               DO 30 J = 1,NT2
                  SUMYK = SUMYK+YKAPPA(KN,2,J)*(H11(KQ,ID)-H)/
     .                    (SKAPPA(KN,2,J)**2-SKAPPA(KN,1,L)**2)
30             CONTINUE
               IF (ITNN.EQ.2) THEN
                  Y2 = Y2-SUMYK*H11(KQ,ID)*B1(KQ,ID)*2.0/
     .                 (RO*(1.0+B))
               ELSE
                  Y2 = Y2+SUMYK*(1.0+B1(KQ,ID))*2.0
     .                 /(RO*(1.0+B))
               ENDIF
               YR = Y2*SKAPPA(KN,1,L)*SKAPPA(KN,1,L)
            ENDIF
            EXP3 = -HA(I)*TMSTT*SKAPPA(KN,1,L)*SKAPPA(KN,1,L)
            IF (ABS(YR).LT.1E-15.AND.YR.GT.0) YR=1E-15
            IF (ABS(YR).LT.1E-15.AND.YR.LT.0) YR=-1E-15
            IF (EXP3.LT.-30.) EXP3=-30.
C           ***************************
            YKAPPA(KN,1,L) = YR*EXP(EXP3)/YR1
            YX=YR/YR1
            VBJ0R = -2.0/(PI*RO)
         TMN = YKAPPA(KN,1,L)*VBJ0R
         IF (ABS(TMN).LT.1E-30.AND.TMN.GT.0) TMN=1E-30
         IF (ABS(TMN).LT.1E-30.AND.TMN.LT.0) TMN=-1E-30
C        **********************************************
         TMN1=YX*VBJ0R
         TERM2 = TERM2+TMN
         TERMM=TERMM-(TMN-TMN1)/(HA(I)*SKAPPA(KN,1,L)*SKAPPA(KN,1,L))
40    CONTINUE
      IF (ITNN.EQ.1) THEN
         VP = TERM1-PI*B*H*(1.0+RO*H*ALOGR)*TERM2
     .        /(WN*(1.0+B))
         WM=TERM1*TMSTT-PI*B*H*(1.0+RO*H*ALOGR)*TERMM/
     .      (WN*(1.0+B))
      ELSE
         VP = TERM1+PI*TERM2
         WM=TERM1*TMSTT+PI*TERMM
      ENDIF
C
C     WTCORR: TEMP. CORR. OF WALL AT STATIONS.
C     WMTR: MASS TRANSFER BETWEEN WALL FILM AND BULK FLOW AT
C           STATIONS IN THE PRESENT INTERVAL, LBM/FT2.
C
      WTCORR(KQ,ID) = VP*(TRV+M)-M
       IF (WTCORR(KQ,ID).LT.0) WTCORR(KQ,ID)=0.
C      ****************************************
      WMTR1=(WM*(TRV+M)-M*TMSTT)*PSII*SLOPE
      WMTR2=(PSII*GAW-GAA)*TMSTT
      WMTR(KQ,ID)=(WMTR1+WMTR2)*BETA*DRR
      IF (KEQUIL.EQ.1) WMTR(KQ,ID)=0.0
C
C  ** SAVE HISTORICAL DATA FOR THE NEXT INTERVAL.
C
      B1(KQ,ID) = B
      H11(KQ,ID) = H
      DO 70 J = 1,NT2
         SKAPPA(KN,2,J) = SKAPPA(KN,1,J)
         YKAPPA(KN,2,J) = YKAPPA(KN,1,J)
70    CONTINUE
C
100   FORMAT (/,' * ATTN * TEMP. CORR. DATA FOR AIRWAY ',I5,' EXCEEDED',
     .       ' THE ARRAY',/,' CAPACITY. IT IS ASSUMED DRY NOW.')
C
      RETURN
      END

