c

      SUBROUTINE FLARGR
      INCLUDE '/phoenics/d_includ/objnam'
      INCLUDE '/phoenics/d_includ/patcmn'
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoeclos/d_includ/d_earth/parvar'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/grdbfc'
      INCLUDE '/phoenics/d_includ/gracm2'
      COMMON /INDAUX/L0ISL,L0IST,L0SL,L0ST,NSTO,NSOL,L0SLRS,L0TTRS,
     1        L0TTFL,L0TTLS,L0MAXC,L0MAXV,L0MINV,L0RATE,L0MAXI,L0NETT
      INCLUDE '/phoenics/d_includ/parear'
      PARAMETER (NLG=20, NIG=20, NRG=100, NCG=10)
      COMMON/LGRND/LG(NLG)/IGRND/IG(NIG)/RGRND/RG(NRG)/CGRND/CG(NCG)
      LOGICAL LG,  QNE, QEQ, QLE
      CHARACTER*4 CG
      CHARACTER*30 FORM
      LOGICAL ISINZZ, LOG1,LOG2, SLD, RPRINT
      COMMON/GENI/NXNY,IGFL(53),IPRPS,IGFL4(4)
      COMMON /CMOBTP/OBJTYP(20)
      CHARACTER*14 OBJTYP
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      LOGICAL LABSV,LTRES,LTRAD,LPMV,LPPD,LSDEN,LSLEN,LPPM,LSMOK,
     1        LPPDR,LPLOS,LTINS,LSLN2,LPPPM,LCOPM
      LOGICAL G_RESTRT,INPARDOM,LSTART,EQZ,LPAR,LTLINK
      LOGICAL LHRAT, LMH2O, LRELH, LPSAT, LXH2O, LPVAP, LOPTD, FLUID1
      REAL KEXT
C... Regression parameters for Productivity loss
      REAL BB(2,7),FVARIN(150),RESID(150)
      CHARACTER*4 NAM(150), BUFF*1024
      CHARACTER*12 COBNAM,PATNM*8
      CHARACTER*256 CTABLE, QTABLE, STABLE
      SAVE CTABLE, QTABLE, STABLE
      SAVE JPMV,JPPD,JTCL,L0ABSV,L0TCLC,L0TINS
      SAVE CL,RMET,WORK,RMW,FCL,RH,TRAD
      SAVE ACON,BCON,CCON,ACN2,BCN2,CCN2,ROX,YCO,YS,KEXT
      SAVE LPPD,LPMV,LPPDR,LPLOS,LTINS
      SAVE LABSV,LTRES,VOLSUM,RESTIME,IFREQ,LTRAD,LSDEN,
     1     LSLEN,LPPM,LSMOK,LSLN2,LPPPM,LCOPM, LOPTD
      SAVE JTEM1,JTRES,JABSV,JT3,JSDEN,JSLEN,JPPM,JSMOK,JSLN2,JPPPM,
     1                                           JCOPM,JOPTD,JAGE,JTLINK
      SAVE JHRAT, JRELH, JPSAT, JMH2O, JXH2O, JPVAP, JPPDR, JPLOS, JTINS
      SAVE LHRAT, LMH2O, LRELH, LPSAT, LXH2O, LPVAP, LTLINK
      COMMON/SPRAY1/ NSPRAY,L0SPRAY,L0SP_TA,L0SP_RTI,L0SP_TL,L0SP_TL0,
     1     L0SP_TG,L0SP_ABV,L0SP_A,ITURB,NSPRAYT,L0SPRAYT,L0ACTIV
      SAVE G_RESTRT,LPAR
      SAVE NFIR,L0FIR,L0QSOR,LUT, LULOCAL
      SAVE NSMOK,L0SMK,L0SSOR
C... Regression parameters for Productivity loss
      DATA BB / 1.2802070,  -0.15397397,
     1         15.995451,    3.8820297,
     2         31.507402,   25.176447,
     3         11.754937,  -26.641366,
     4          1.4737526,  13.110120,
     5          0.0,        -3.1296854,
     6          0.0,         0.29260920  /
C
      IXL=IABS(IXL)
C... Call HOTBOX to calculate fan operating point or system curve
      CALL HTBXGR
      NAMSUB='FLARGR'
      IF(IGR.EQ.1) GO TO 1
      IF(IGR.EQ.13) GO TO 13
      IF(IGR.EQ.19) GO TO 19
      RETURN
C*****************************************************************
C
C--- GROUP 1. Run title and other preliminaries
C
    1 GO TO (1001,1002,1003),ISC
 1001 CONTINUE
      IF(.NOT.NULLPR) THEN
        CALL WRYT40('Call to spec. Ground FLARGR.F of: 130412')
        CALL WRIT40('Call to spec. Ground FLARGR.F of: 130412')
      ENDIF
C ... Initialise
      OBJTYP(5)='OPENING'
      LPAR=MIMD.AND.NPROC.GT.1
C... Spray settings
      LTLINK=.FALSE.; JTLINK=0
      IF(.NOT.STEADY) THEN
C... loop over allpatches looking for SPRAY
        NSPRAY=0; NSPRAYT=0
        DO IP=1,NUMPAT
          IF(NAMPAT(IP)(1:5).EQ.'SPRAY') THEN
            NSPRAYT=NSPRAYT+1
C... check if link temperature active
            TACTIVE=-999.
            CALL GETSDR(NAMPAT(IP),'T_ACTIV',TACTIVE)
            IF(QNE(TACTIVE,-999.)) THEN
C... increment spray counter
              NSPRAY=NSPRAY+1
            ENDIF
          ENDIF
        ENDDO
        IF(NSPRAY.GT.0) THEN
C... there are some sprays with active link temperatures.
          JTLINK=LBNAME('TLNK'); LTLINK=STORE(JTLINK)
          IF(MYID.EQ.0.AND..NOT.LTLINK) THEN
            SCROLL_BUFFER(1)='WARNING: This run requires STORE(TLNK)'
            SCROLL_BUFFER(2)='in Q1 otherwise will not be able to'
            SCROLL_BUFFER(3)='perform smooth restart'
            DO I=1,3
              CALL WRIT40(SCROLL_BUFFER(I))
            ENDDO
            SCROLL_BUFFER(4)=
     1               'Press Yes to continue anyway or No to abort'
            NLBUF=4
            CALL DISPLAY_MESSAGE(3, IRES)
            IF(IRES.NE.1) THEN
              CALL SET_ERR(580,'Error. No TLNK store for Flair spray',1)
            ENDIF
            NLBUF=0
          ENDIF
          DO IP=1,NUMPAT
            IF(NAMPAT(IP)(1:5).EQ.'SPRAY') THEN
              CALL MAKEPV(30,IP)
            ENDIF
          ENDDO
        ENDIF
      ELSE
        NSPRAY=0; NSPRAYT=0
        G_RESTRT=.FALSE.
      ENDIF
      RETURN
 1002 CONTINUE
C ... Humidity variables
      JHRAT =LBNAME('HRAT') ! humidity ratio (g/kg)
      JRELH =LBNAME('RELH') ! relative humidity
      JPSAT =LBNAME('PSAT') ! water-vapour saturation pressure
      JMH2O =LBNAME('MH2O') ! mass fraction of water vapour
      JXH2O =LBNAME('XH2O') ! water-vapour mole-fraction
      JPVAP =LBNAME('PVAP') ! water vapour partial pressure
      JPPDR =LBNAME('PPDR') ! Draught Rating
      JTINS =LBNAME('TINS') ! Turbulence intensity
      LHRAT =STORE(JHRAT)
      LRELH =STORE(JRELH)
      LPSAT =STORE(JPSAT)
      LMH2O =STORE(JMH2O)
      LXH2O =STORE(JXH2O)
      LPVAP =STORE(JPVAP)
      LPPDR =STORE(JPPDR)
      LTINS =STORE(JTINS)
C
      JTEM1=LBNAME('TEM1') ! temperature
      JTRES=LBNAME('TRES') ! dry resultant temperature
      JT3  =LBNAME('T3')   ! radiant temperature
      JABSV=LBNAME('ABSV') ! absolute velocity
      JSMOK=LBNAME('SMOK') ! smoke concentration
      JSDEN=LBNAME('SDEN') ! smoke density
      JSLEN=LBNAME('SLEN') ! sight length
      JSLN2=LBNAME('SLN2') ! sight length 2
      JPPM =LBNAME('PPM')  ! smoke products parts-per-million
      JPPPM =LBNAME('PPPM') ! smoke particulates parts-per-million
      JCOPM =LBNAME('COPM') ! CO parts-per-million
      JOPTD =LBNAME('OPTD') ! Optical density
      JAGE  =LBNAME('AGE')  ! Age of air
      LABSV=STORE(JABSV)
      LTRES=STORE(JTRES)
      LTRAD=STORE(JT3)
      LSMOK=STORE(JSMOK)
      LSDEN=STORE(JSDEN).AND.LSMOK
      LSLEN=STORE(JSLEN).AND.LSMOK
      LSLN2=STORE(JSLN2).AND.LSMOK
      IF(LSMOK) THEN
        IF(LSLEN) THEN ! constants for sight length
          CALL GETSDR('SLEN','ACON',ACON)
          BCON=-999.0
          CALL GETSDR('SLEN','BCON',BCON)
          CALL GETSDR('SLEN','CCON',CCON)
        ENDIF
        IF(LSLN2) THEN ! constants for sight length 2
          CALL GETSDR('SLN2','ACON',ACN2)
          BCN2=-999.0
          CALL GETSDR('SLN2','BCON',BCN2)
          CALL GETSDR('SLN2','CCON',CCN2)
        ENDIF
        LPPM =STORE(JPPM)
        LPPPM =STORE(JPPPM)
        LCOPM =STORE(JCOPM)
        LOPTD =STORE(JOPTD)
        HFU=25.E6
        CALL GETSDR('FLAIR','HCOMB',HFU) ! Heat of combustion
        ROX=HFU/13.1E6
        CALL GETSDR('FLAIR','STOIC',ROX) ! Stoichiometric ratio
        YS=0.157
        CALL GETSDR('FLAIR','SYIELD',YS) ! Smoke particulate yield
        YCO=0.06
        CALL GETSDR('FLAIR','COYIELD',YCO) ! CO yield
        KEXT=7600.
        CALL GETSDR('FLAIR','KEXT',KEXT) !
        IF(BCON.LE.0) BCON=YS*KEXT/(1.0+ROX)
        IF(BCN2.LE.0) BCN2=YS*KEXT/(1.0+ROX)
      ELSE
        LSLEN=.FALSE.
        LSLN2=.FALSE.
        LPPM=.FALSE.
        LPPPM=.FALSE.
        LCOPM=.FALSE.
        LOPTD =.FALSE.
      ENDIF
      IFREQ=IDISPA
      IF(IFREQ.EQ.0) IFREQ=LSWEEP+1
CFREQ      IFREQ=IDISPA
CFREQ      IF(IFREQ.EQ.0) IFREQ=LSWEEP+1
      JPMV=LBNAME('PMV')  ! Predicted Mean Vote
      LPMV=JPMV.GT.0
      IF(LPMV) THEN
        JPPD=LBNAME('PPD')  ! Predicted Percentage Dissatified
        LPPD=JPPD.GT.0
        JTCL=LBNAME('TCL')  ! Clothing temperature in deg C
        IF(JTCL.EQ.0) THEN
C... If clothing temperature not STOREd, make 3D store
          CALL GXMAKE0(L0TCLC,NX*NY*NZ,'TCL')
        ENDIF
        JPLOS =LBNAME('PLOS') ! Productivity Loss
        LPLOS =STORE(JPLOS)
c Get SPEDAT settings for PMV parameters (CL,RMET,WORK,RH)
c and deduce others (FCL,RMW)
        CL=-999.
        CALL GETSDR('FLAIR','CLOTHING',CL)
        IF(CL.LT.0.0) THEN
          CALL WRIT40('Invalid value for clothing insulation   ')
          CALL WRYT40('Invalid value for clothing insulation   ')
          CALL SET_ERR(518,
     *          'Error. Invalid value for clothing insulation.',1)
c          CALL WAYOUT(1)
        ENDIF
        IF(CL.LT.0.0078) THEN
          FCL=1.0+1.29*CL
        ELSE
          FCL=1.05+0.645*CL
        ENDIF
        RMET=-999.
        CALL GETSDR('FLAIR','METABOLIC',RMET)
        IF(RMET.LT.0.0) THEN
          CALL WRIT40('Invalid value for metabolic rate        ')
          CALL WRYT40('Invalid value for metabolic rate        ')
          CALL SET_ERR(519,
     *          'Error. Invalid value for metabolic rate.',1)
        ENDIF
        WORK=-999.
        CALL GETSDR('FLAIR','EXTWORK',WORK)
        IF(WORK.LT.0.0) THEN
          CALL WRIT40('Invalid value for external work         ')
          CALL WRYT40('Invalid value for external work         ')
          CALL SET_ERR(520,
     *          'Error. Invalid value for external work.',1)
        ENDIF
        RMW=RMET-WORK
        IF(RMW.LT.0.0) THEN
          CALL WRIT40('Metabolic rate must exceed external work')
          CALL WRYT40('Metabolic rate must exceed external work')
          CALL SET_ERR(521,
     *          'Error. Metabolic rate must exceed external work.',1)
        ENDIF
        RH=-999.
        CALL GETSDR('FLAIR','RELHUMID',RH)
        IF((RH.LT.0.0.OR.RH.GT.100.).AND..NOT.LRELH) THEN
          CALL WRIT40('Invalid value for relative humidity     ')
          CALL WRYT40('Invalid value for relative humidity     ')
          CALL SET_ERR(522,
     *          'Error. Invalid value for relative humidity.',1)
        ENDIF
      ENDIF
C... Radiant temperature
      IF(LTRES.OR.LPMV) THEN
        TRAD=-999.0
        CALL GETSDR('FLAIR','TRAD',TRAD)
        IF(QNE(TRAD,-999.0)) THEN
          LTRAD=.FALSE.
        ELSE
          IF(.NOT.LTRAD) THEN
            CALL WRIT40('No valid setting made for the radiant   ')
            CALL WRIT40('temperature - user-set or Immersol!     ')
            CALL WRYT40('No valid setting made for the radiant   ')
            CALL WRYT40('temperature - user-set or Immersol!     ')
            CALL SET_ERR(523, 'Error. See result file.',1)
          ENDIF
        ENDIF
      ENDIF
C... Draught rating
      IF(LPPDR) THEN
        IF(.NOT.LTINS) CALL GXMAKE0(L0TINS,NX*NY,'TINS')
        IF(SOLVE(KE)) THEN
          ITURB=1
        ELSEIF(STORE(VIST).AND.STORE(LEN1)) THEN
          ITURB=2
        ELSE
          ITURB=0
        ENDIF
      ENDIF
C... Sprinkler settings
      IF(.NOT.STEADY) THEN
        IF(NSPRAY.GT.0) THEN
          WRITE(LUPR1,*) 'Link temperature calculation',
     1                                           ' has been activated '
        ENDIF
      ENDIF
C... If absolute velocity not STOREd, make 3D store
      IF(.NOT.LABSV.AND.(LTRES.OR.LPMV.OR.NSPRAY.GT.0)) THEN
        CALL GXMAKE0(L0ABSV,NX*NY*NZ,'ABSV')
      ENDIF
C
      LUT=87
      CTABLE='conv_table.csv'
C.. Open csv table file and write headings
      IF(STEADY) THEN
        LSTART=FSWEEP.EQ.1
      ELSE
        LSTART=FSTEP.EQ.1
      ENDIF
      ISTO=0
      DO II=1,NPHI
        IF(SOLVE(II).AND.NAME(II).NE.'LTLS') THEN
          ISTO=ISTO+1
          NAM(ISTO)=NAME(II)
        ENDIF
      ENDDO
      IF(MYID.EQ.0) CALL MON_TABL(1,LUT,CTABLE,LSTART,' ') ! Open
      IF(MYID.EQ.0.AND.LSTART) THEN
        IF(STEADY)THEN
          WRITE(BUFF,'(''isweep'',150('',  '',A,5X))')
     1                                                 (NAM(I),I=1,ISTO)
        ELSE
          WRITE(BUFF,'('' istep,    time,   isweep'',
     1                            150('',  '',A,5X))') (NAM(I),I=1,ISTO)
        ENDIF
        CALL MON_TABL(2,LUT,CTABLE,.FALSE.,BUFF) ! Write
      ENDIF
C
C.... prepare heat and mass source tables for FIRE objects
      IF(.NOT.STEADY) THEN
        NFIR0=0; NSMOK0=0; ILIM=ITWO(GD_NUMPAT,NUMPAT,LPAR)
        DO IP=1,ILIM
          IF(LPAR) THEN
            PATNM=GD_NAMPAT(IP)
          ELSE
            PATNM=NAMPAT(IP)
          ENDIF
          IF(PATNM(1:3).EQ.'FIR') THEN
            LL=LENGZZ(PATNM)
            IF(PATNM(LL:LL).EQ.'A') NFIR0=NFIR0+1
            IF(PATNM(LL:LL).NE.'A') NSMOK0=NSMOK0+1
          ENDIF
        ENDDO
        IF(NFIR0.GT.0) THEN
          CALL GXMAKE0(L0FIR,NFIR0,'FIRE')
          CALL GXMAKE0(L0QSOR,NFIR0,'QSOR')
          QTABLE='heat_sources.csv'
          LSTART=FSTEP.EQ.1
          IF(MYID.EQ.0) CALL MON_TABL(1,LUT,QTABLE,LSTART,' ') ! Open
          NFIR=0
          DO IP=1,ILIM
            IF(LPAR) THEN
              PATNM=GD_NAMPAT(IP)
            ELSE
              PATNM=NAMPAT(IP)
            ENDIF
            IF(PATNM(1:3).EQ.'FIR') THEN
              LL=LENGZZ(PATNM)
              IF(PATNM(LL:LL).EQ.'A') THEN
                NFIR=NFIR+1
                F(L0FIR+NFIR)=IP
                IF(NFIR.EQ.NFIR0) EXIT
              ENDIF
            ENDIF
          ENDDO
          IF(LSTART.AND.MYID.EQ.0) THEN
            BUFF=' Time     '
            LL=LENGZZ(BUFF)+5
            DO IFIR=1,NFIR
              IP=NINT(F(L0FIR+IFIR))
              IF(LPAR) THEN
                COBNAM=GD_OBJNAM(IP)
              ELSE
                COBNAM=OBJNAM(IP)
              ENDIF
              BUFF(LL+1:)=', '//COBNAM(1:9)
              LL=LL+11
            ENDDO
            CALL MON_TABL(2,LUT,QTABLE,.FALSE.,BUFF) ! Write
          ENDIF
        ENDIF
        IF(NSMOK0.GT.0) THEN
          CALL GXMAKE0(L0SMK,NSMOK0,'SMK')
          CALL GXMAKE0(L0SSOR,NSMOK0,'SSOR')
          STABLE='smoke_sources.csv'
          LSTART=FSTEP.EQ.1
          IF(MYID.EQ.0) CALL MON_TABL(1,LUT,STABLE,LSTART,' ')
          NSMOK=0
          DO IP=1,ILIM
            IF(LPAR) THEN
              PATNM=GD_NAMPAT(IP)
            ELSE
              PATNM=NAMPAT(IP)
            ENDIF
            IF(PATNM(1:3).EQ.'FIR') THEN
              LL=LENGZZ(PATNM)
              IF(PATNM(LL:LL).NE.'A') THEN
                NSMOK=NSMOK+1
                F(L0SMK+NSMOK)=IP
                IF(NSMOK.EQ.NSMOK0) EXIT
              ENDIF
            ENDIF
          ENDDO
          IF(LSTART.AND.MYID.EQ.0) THEN
            BUFF=' Time     '
            LL=LENGZZ(BUFF)+5
            DO ISMOK=1,NSMOK
              IP=NINT(F(L0SMK+ISMOK))
              IF(LPAR) THEN
                COBNAM=GD_OBJNAM(IP)
              ELSE
                COBNAM=OBJNAM(IP)
              ENDIF
              BUFF(LL+1:)=', '//COBNAM(1:9)
              LL=LL+11
            ENDDO
            CALL MON_TABL(2,LUT,STABLE,.FALSE.,BUFF) ! Write
          ENDIF
        ENDIF
      ENDIF
      RETURN
C
 1003 CONTINUE
      IF(NSPRAY.GT.0) THEN
C... there are some sprays with active link temperatures.
C... reserve storage
        CALL GXMAKE0(L0ACTIV,NSPRAYT,'L0ACTIV')
        CALL GXMAKE0(L0SPRAYT,NSPRAY,'L0SPRAYT')
        CALL GXMAKE0(L0SPRAY,NSPRAY,'L0SPRAY')
        CALL GXMAKE0(L0SP_TA,NSPRAY,'SP_TA')
        CALL GXMAKE0(L0SP_RTI,NSPRAY,'SP_RTI')
        CALL GXMAKE0(L0SP_TL,NSPRAY,'SP_TL')
        CALL GXMAKE0(L0SP_TL0,NSPRAY,'SP_TL0')
        CALL GXMAKE0(L0SP_TG,NSPRAY,'SP_TG')
        CALL GXMAKE0(L0SP_ABV,NSPRAY,'SP_ABV')
        CALL GXMAKE0(L0SP_A,NSPRAY,'SP_A')
C... loop round patches again, storing information
        ISPRAY=0; ISPRAYT=0
        DO IP=1,NUMPAT
          IF(NAMPAT(IP)(1:5).EQ.'SPRAY') THEN
            ISPRAYT=ISPRAYT+1; TACTIVE=-999.
            CALL GETSDR(NAMPAT(IP),'T_ACTIV',TACTIVE)
            IF(QNE(TACTIVE,-999.)) THEN
              ISPRAY=ISPRAY+1
              F(L0SPRAYT+ISPRAY)=FLOAT(ISPRAYT)
              F(L0SPRAY+ISPRAY)=FLOAT(IP)
              F(L0SP_TA+ISPRAY)=TACTIVE
              CALL GETSDR(NAMPAT(IP),'RTI',RTI)
              F(L0SP_RTI+ISPRAY)=RTI
            ENDIF
          ENDIF
        ENDDO
        G_RESTRT=QEQ(FIINIT(P1),READFI)
      ENDIF
      RETURN
C*****************************************************************
C
C--- GROUP 13. Boundary conditions and special sources
C                                       Index for Coefficient - CO
C                                       Index for Value       - VAL
   13 IF(ISC.EQ.12) THEN
C------------------- SECTION 12 ------------------- value = GRND
        IF(NPATCH(1:2).EQ.'TM') THEN
C... TM Patches; does not use GXTIM to allow different amplitudes
C....GRND....Linear increase from 0 to Amplitude (carried by RG(40+i))
C... Up to 999 patches allowed
          ILEN=0
          IF(ISINZZ(NPATCH(3:3))) ILEN=1
          IF(ILEN.EQ.1.AND.ISINZZ(NPATCH(4:4))) ILEN=2
          IF(ILEN.EQ.2.AND.ISINZZ(NPATCH(5:5))) ILEN=3
          IF(ILEN.EQ.0) GO TO 999
          FORM='(I1)'
          WRITE(FORM(3:3),'(I1)') ILEN
          READ(NPATCH(3:3+ILEN-1),FORM,ERR=999) IADD
C... Get limits
          IFST=IPATNF(9,IREG)
          ILST=IPATNF(10,IREG)
C... Calculate 'steps'
          RCUT=ISTEP-IFST+0.5
          COEF1=RG(40+IADD)/(ILST-IFST+1)
          SETHS=RCUT*COEF1
C... Set heat sourse as constant over the time-step
          CALL FN1(VAL,SETHS)
        ENDIF
      ELSEIF(ISC.EQ.13) THEN
        IF(INDVAR.EQ.JMH2O) THEN
          RELH=-999.0
          CALL GETSDR(OBJNAM(IPAT),'REL-HUM',RELH)
          IF(QNE(RELH,-999.0)) THEN
            L0VAL=L0F(VAL)
            L0P1=L0F(P1)
            L0TEM1=L0F(JTEM1)
            GWRAT = 18.015/28.96
            DO IX=IXF,IXL
              DO IY=IYF,IYL
                I=(IX-1)*NY+IY
                TIN=F(L0TEM1+I)+TEMP0-273.0   ! in deg C
                PIN=F(L0P1+I)+PRESS0
                PVAP = 1E-2*RELH*PVH2O(TIN)   ! water-vapour pressure
                HRAT = AMAX1(0.0,AMIN1(GWRAT*PVAP/(PIN-PVAP),1.0))
                F(L0VAL+I)=HRAT/(1.0+HRAT)    ! convert to mass fraction
              ENDDO
            ENDDO
          ENDIF
        ENDIF
      ENDIF
      RETURN
C***************************************************************
C
C--- GROUP 19. Special calls to GROUND from EARTH
C
   19 IF(ISC.EQ.1) THEN
C   * ------------------- SECTION 1 ---- Start of time step.
C
C.. Initialise link temperature on 1st time step
        IF(.NOT.STEADY.AND.NSPRAY.GT.0) THEN
          IF(ISTEP.EQ.FSTEP)  THEN
            LOG1=.FALSE.
            IF(G_RESTRT) THEN
C... recover old sprinkler link temperatures from patch-wise store
              DO JSP=1,NSPRAY
                IPT=NINT(F(L0SPRAY+JSP))
                CALL GETPTC(NAMPAT(IPT),GTYPE,JXF,JXL,JYF,JYL,JZF,JZL,
     1                                                          JTF,JTL)
!                L0PAT= L0PVAR(30,IPT, JZF)
!                F(L0SP_TL0+JSP) = F(L0PAT+1)
                IF(LTLINK) THEN
                  L0TLINK=L0F(ANYZ(JTLINK,JZF)); IJ=(JXF-1)*NY+JYF
                  F(L0SP_TL0+JSP) = F(L0TLINK+IJ)
                ENDIF
                IF(F(L0SP_TL0+JSP).GT.F(L0SP_TA+JSP))  THEN
                  WRITE(14,*) '!!!! ',OBJNAM(IPT),' active !!!!'
                  ISP=NINT(F(L0SPRAYT+JSP))
                  F(L0ACTIV+ISP)=1.; LOG1=.TRUE.
                ENDIF
              ENDDO
            ELSE
C... initialise old sprinkler link temperatures to FIINIT
              DO JSP=1,NSPRAY
                F(L0SP_TL0+JSP) = FIINIT(JTEM1)
                IF(F(L0SP_TL0+JSP).GT.F(L0SP_TA+JSP))  THEN
                  IPT=NINT(F(L0SPRAY+JSP))
                  WRITE(14,*) '!!!! ',OBJNAM(IPT),' active !!!!'
                  ISP=NINT(F(L0SPRAYT+JSP))
                  F(L0ACTIV+ISP)=1.; LOG1=.TRUE.
                ENDIF
              ENDDO
            ENDIF
            IF(LOG1) CALL ACTIVATE_SPRAYS
          ENDIF
        ENDIF
      ELSEIF(ISC.EQ.6) THEN
C   * ------------------- SECTION 6 ---- Finish of iz slab.
C ... Calculate absolute velocities and store magnitude in a
C     variable ABSV for analysis of 'dead' circulation zones
C
        LOG1=.FALSE.
        IF(STEADY) THEN
          LOG1=((MOD(ISWEEP+1,IFREQ).EQ.0.AND.CSG1.EQ.'SW') .OR.
     1           MOD(ISWEEP+1,NPRINT).EQ.0) .OR.
     1           ENUFSW.OR.ISWEEP.EQ.LSWEEP-1.OR.ISWEEP.EQ.LSWEEP
        ELSE
          LOG1= (MOD(ISTEP,IFREQ).EQ.0 .OR.
     1           MOD(ISTEP,NTPRIN).EQ.0) .AND.
     1          (ENUFSW.OR.ISWEEP.EQ.LSWEEP-1.OR.ISWEEP.EQ.LSWEEP)
        ENDIF
        IF((LTRES.AND.LOG1).OR.(LPMV.AND.INDVAR.EQ.JTEM1).OR.
     1                                                NSPRAY.GT.0) THEN
          IF(STORE(U1)) L0U1=L0F(U1)
          IF(STORE(V1)) L0V1=L0F(V1)
          IF(STORE(W1)) THEN
            L0W1 =L0F(W1)
            L0W1L=L0F(LOW(W1))
          ENDIF
          IF(LABSV) THEN
            L0AVL=L0F(JABSV)
c.... set liter as a sign that value has been changed
            LITER(JABSV)=1
          ELSE
            L0AVL=L0ABSV+(IZ-1)*NX*NY
          ENDIF
          DO IX=1,NX
            DO IY=1,NY
              I=IY+NY*(IX-1)
C...  Staggered grid velocity components
C...  Average velocities to cell centre before squaring
              IF(NX.GT.1) THEN
                IF(IX.EQ.1) THEN
                  AU1=F(L0U1+I)**2
                ELSEIF(IX.EQ.NX) THEN
                  AU1=F(L0U1+I-NY)**2
                ELSE
                  AU1=(0.5*(F(L0U1+I)+F(L0U1+I-NY)))**2
                ENDIF
              ELSE
                AU1=0.0
              ENDIF
              IF(NY.GT.1) THEN
                IF(IY.EQ.1) THEN
                  BV1=F(L0V1+I)**2
                ELSEIF(IY.EQ.NY) THEN
                  BV1=F(L0V1+I-1)**2
                ELSE
                  BV1=(0.5*(F(L0V1+I)+F(L0V1+I-1)))**2
                ENDIF
              ELSE
                BV1=0.0
              ENDIF
              IF(NZ.GT.1) THEN
                IF(IZ.EQ.1) THEN
                  CW1=F(L0W1+I)**2
                ELSEIF(IZ.EQ.NZ) THEN
                  CW1=F(L0W1L+I)**2
                ELSE
                  CW1=(0.5*(F(L0W1+I)+F(L0W1L+I)))**2
                ENDIF
              ELSE
                CW1=0.0
              ENDIF
C... Set velocity
              F(L0AVL+I)=SQRT(AU1+BV1+CW1)
            ENDDO
          ENDDO
        ENDIF
        IF(LOG1) THEN
C... Calculate overall 'free' volume
          L0VL=L0F(VOL)
          IF(IZ.EQ.1) VOLSUM=0.
          DO IX=1,NX
            DO IY=1,NY
              I=(IX-1)*NY+IY
              IF(.NOT.SLD(I)) THEN
                IF(NPROC.GT.1) THEN
                  IF(INPARDOM(IX,IY,IZ)) VOLSUM=VOLSUM+F(L0VL+I)
                ELSE
                  VOLSUM=VOLSUM+F(L0VL+I)
                ENDIF
              ENDIF
            ENDDO
          ENDDO
        ENDIF
C... Calculate Smoke density in kg/m^3 of mixture
        IF(LSDEN) CALL FN026(JSDEN,DEN1,JSMOK) ! JSDEN=DEN1*JSMOK
C... Calculate Smoke product density in parts/million
        IF(LPPM) CALL FN2(JPPM,JSMOK,0.0,1.E6) ! JPPM=0+1e6*JSMOK
C... Calculate Smoke particulate density in parts/million
        IF(LPPPM) CALL FN2(JPPPM,JSMOK,0.0,1.E6*YS/(1.+ROX))
C... Calculate CO density in parts/million
        IF(LCOPM) CALL FN2(JCOPM,JSMOK,0.0,1.E6*YCO/(1.+ROX))
C... Calculate Optical density in 1/m
        IF(LOPTD) CALL FN21(JOPTD,JSMOK,DEN1,0.0,KEXT*YS/(2.3*(1.+ROX)))
C... Calculate sight length based on
C           Sl = max(0, min(C, A/(B*Smoke Density)))
        IF(LSLEN) THEN
          L0SLEN=L0F(JSLEN)
          L0DEN1=L0F(DEN1)
          L0SMOK=L0F(JSMOK)
          DO I=1,NX*NY
            IF(.NOT.SLD(I)) THEN
              F(L0SLEN+I)=
     1  AMAX1(0.0, AMIN1(CCON,ACON/(BCON*F(L0SMOK+I)*F(L0DEN1+I)+TINY)))
            ELSE
              F(L0SLEN+I)=0.0
            ENDIF
          ENDDO
        ENDIF
        IF(LSLN2) THEN
          L0SLEN=L0F(JSLN2)
          L0DEN1=L0F(DEN1)
          L0SMOK=L0F(JSMOK)
          DO I=1,NX*NY
            IF(.NOT.SLD(I)) THEN
              F(L0SLEN+I)=
     1  AMAX1(0.0, AMIN1(CCN2,ACN2/(BCN2*F(L0SMOK+I)*F(L0DEN1+I)+TINY)))
            ELSE
              F(L0SLEN+I)=0.0
            ENDIF
          ENDDO
        ENDIF
C
C .... Humidity parameters
        IF(INDVAR.EQ.JTEM1.AND.LMH2O) THEN
C .. Humidity ratio (g/kg)
          IF(LHRAT.OR.LRELH) THEN
            GMWRAT = 18.015/28.96
            IF(LHRAT) L0HRAT = L0F(JHRAT)
            IF(LRELH) L0RELH = L0F(JRELH)
            L0MH2O = L0F(JMH2O)
            L0TEM1=L0F(JTEM1)
            L0P1 = L0F(P1)
            IF(LPSAT) L0PSAT=L0F(JPSAT)
            IF(LPVAP) L0PVAP=L0F(JPVAP)
            IF(LXH2O) L0XH2O=L0F(JXH2O)
            DO IX=1,NX
              DO IY=1,NY
                I=IY+NY*(IX-1)
                IF(.NOT.SLD(I)) THEN
C .. Humidity ratio (g/kg)
                  YH2O = F(L0MH2O+I)
                  YH2O=AMIN1(1.0,AMAX1(0.0,YH2O))
C                  HRAT = 1.E3*YH2O/(1.-YH2O+1.E-10)
C                  IF(LHRAT) F(L0HRAT+I) = HRAT
                  IF(QEQ(YH2O,1.0)) THEN
c                    WRITE(14,*) 'IX,IY,IZ, HRAT ',IX,IY,IZ, HRAT
                    HRAT=1./TINY
                  ELSE
                    HRAT = YH2O/(1.-YH2O+1.E-10)
                  ENDIF
                  IF(LHRAT) F(L0HRAT+I) = 1.E3*HRAT
                  IF(LRELH) THEN
C .. Water Vapour Saturation Pressure (Pa)
                    TIN=F(L0TEM1+I)+TEMP0-273.0   ! in deg C
                    PVSAT = PVH2O(TIN)
                    IF(LPSAT) F(L0PSAT+I)=PVSAT
C .. Water vapour Pressure (Pa)
C                    HRAT = 1.E-3*HRAT
                    PVAP = HRAT*(PRESS0+F(L0P1+I))/(GMWRAT+HRAT)
                    IF(LPVAP) F(L0PVAP+I)=PVAP
C .. Relative humidity calculation (%)
                    RELH = AMAX1(0.0, AMIN1(100.0, 100.*PVAP/PVSAT))
                    F(L0RELH+I) = RELH
C .. Water vapour mole fraction
                    IF(LXH2O) THEN
                      WMIX = YH2O*18.015 + (1.-YH2O)*28.96
                      F(L0XH2O+I) = YH2O*WMIX/18.015
                    ENDIF
                  ENDIF
                ENDIF
              ENDDO
            ENDDO
          ENDIF
        ENDIF
C... Calculate comfort indices if PMV stored
c Coding based on ISO 7730 BASIC program
c Variables used:
c   Input:
c     CL: thermal insulation of clothing (m2K/W)
c     RMET: metabolic rate (W/m2)
c     WORK: external work (W/m2)
c     RH: relative humidity (%)
c     TAA: air temperature (K)
c     TRA: radiant temperature (K)
c   Used in calculation:
c     FCL: clothing factor (area ratio of clothed to unclothed)
c     TCLA: surface temperature of clothing (K)
c     TCLC: surface temperature of clothing (C)
c     XN: TCLA/100
c     HCF: heat transfer coefficient by forced convection (W/m2K)
c     HC: actual heat transfer coefficient (natural or forced) (W/m2K)
c     PA: partial vapour pressure (Pa)
c     VAR: relative air velocity (m/s)
c   Output:
c     PMV: predicted mean vote (%)
c     PPD: predicted percentage dissatisfied (%) - optional
c     HL1: heat loss diffused through skin (?) - not yet available
c     HL2: heat loss by sweating (?) - not yet available
c     HL3: latent respiration heat loss (?) - not yet available
c     HL4: dry respiration heat loss (?) - not yet available
c     HL5: radiation heat loss (?) - not yet available
c     HL6: convective heat loss (?) - not yet available
c
        IF(INDVAR.EQ.JTEM1) THEN
          IF(LPMV) THEN
            IF(LABSV) THEN
              L0AVL=L0F(JABSV)
            ELSE
              L0AVL=L0ABSV+(IZ-1)*NX*NY
            ENDIF
            L0TEM=L0F(JTEM1)
            IF(LTRAD) THEN
              L0T3=L0F(JT3)
            ENDIF
            IF(JTCL.GT.0) THEN
              L0TCL=L0F(JTCL)
            ELSE
              L0TCL=L0TCLC+(IZ-1)*NX*NY
            ENDIF
            L0PMV=L0F(JPMV)
            IF(LPPD) L0PPD=L0F(JPPD)
            IF(LPLOS) L0PLOS=L0F(JPLOS)
            IF(QEQ(RH,-999.0).AND.LRELH) THEN
              L0RELH=L0F(JRELH)
            ELSE
              L0RELH=0
            ENDIF
C... extra air velocity due to metabolic rate (movement) ISO 7730 Annex D
            RRMET=RMET/58.15 ! Metabolic rate in met
            IF(RRMET.LT.1.0) THEN
              VMETAB=0.0
            ELSE
              VMETAB=0.3*(RRMET-1.0)
            ENDIF
            DO IX=1,NX
              DO IY=1,NY
                I=IY+NY*(IX-1)
                IF(.NOT.SLD(I)) THEN
C Get absolute velocity from store - add movement-induced component
                  VAR=F(L0AVL+I)+VMETAB
                  HCF=12.1*SQRT(VAR)
                  TAA=F(L0TEM+I)+TEMP0
                  TAA=AMIN1(373.0,TAA)
                  IF(LTRAD) THEN
                    TRA=F(L0T3+I)+TEMP0
                  ELSE
                    TRA=TRAD+273.  ! TRAD is in deg C
                  ENDIF
                  TCLA=F(L0TCL+I)+TEMP0
                  TCLC=TCLA-273.0
C Define other terms in the equation
                  IF(L0RELH.EQ.0) THEN
                    RELH=RH
                  ELSE
                    RELH=F(L0RELH+I)
                  ENDIF
                  PA=10.*RELH*EXP(16.6536-4030.183/(TAA-38.))
                  PP1=CL*FCL
                  PP2=3.96*PP1
                  P3=100.*PP1
                  P4=PP1*TAA
                  P5=308.7-0.028*RMW+PP2*(TRA/100.)**4
C Initialise iterative loop
                  XN=TCLA/100.
                  XF=XN
                  N=0
                  EPS=.0000001 ! convergence criterion
C Start of iterations
19611             CONTINUE
                  XF=0.5*(XF+XN)
                  HC=2.38*(ABS(100.*XF-TAA))**0.25
                  IF(HC.LT.HCF) HC=HCF
                  XN=(P5 + P4*HC - PP2*XF**4)/(100.+P3*HC)
                  N=N+1
CFREQ              IF(ABS(XN-XF).GT.EPS.AND.N.LE.150) GO TO 19611
                  IF(ABS(XN-XF).GT.EPS.AND.N.LE.10) GO TO 19611
C End of iterations, so store result
C (Ought to write message about non-convergence, at least on
C  LSWEEP - and probably LSWEEP-1 as well)
                  F(L0TCL+I)=XN*100.-273. ! XN in K, but TCL in deg C
C Calculate the various heat transfer contributions ready for PMV
                  HL1=3.05*0.001*(5733.-6.99*RMW-PA) ! loss through skin
                  IF(RMW.GT.58.15) THEN
                    HL2=0.42*(RMW-58.15)  ! loss by sweating
                  ELSE
                    HL2=0.0
                  ENDIF
                  HL3=1.7*0.00001*RMET*(5867.-PA) ! latent perspiration loss
                  HL4=0.0014*RMET*(34.-(TAA-273.)) ! dry perspiration loss
                  HL5=3.96*FCL*(XN**4-(TRA/100.)**4) ! radiative loss
                  HL6=FCL*HC*(XN*100.-TAA)           ! convective heat loss
                  TS=0.303*EXP(-0.036*RMET)+0.028 ! thermal sensation tran coeff
                  PMV=TS*(RMW-HL1-HL2-HL3-HL4-HL5-HL6) ! Predicted mean vote
                  PMV=AMAX1(-3.0,AMIN1(3.0,PMV)) ! Limit to +/-3.0
                  F(L0PMV+I)=PMV
                  IF(LPPD) THEN
                    PMV2=PMV**2.0
                    F(L0PPD+I)=100.-95.*EXP((-0.03353*PMV2-0.2179)*PMV2)
                    F(L0PPD+I)=AMAX1(5.0,AMIN1(100.0,F(L0PPD+I))) ! Limit
                  ENDIF
C... Productivity loss
                  IF(LPLOS) THEN
                    PMV=AMAX1(-3.0, AMIN1(PMV,3.0)) ! limit to +/ 3
                    II=ITWO(1,2,PMV.LE.0.0) ! cold/hot side
                    PLOS=BB(II,1)+BB(II,2)*PMV+BB(II,3)*PMV**2.0+
     1                   BB(II,4)*PMV**3.0+BB(II,5)*PMV**4.0+
     1                   BB(II,6)*PMV**5.0+BB(II,7)*PMV**6.0
                    F(L0PLOS+I)=AMAX1(0.0,AMIN1(PLOS,100.0))
                  ENDIF
                ENDIF
              ENDDO
            ENDDO
          ENDIF
C... Draught rating (as per ISO 7730)
          IF(LPPDR) THEN
            L0PPDR=L0F(JPPDR)
            IF(LTINS) L0TINS=L0F(JTINS)
            IF(LABSV) THEN
              L0AVL=L0F(JABSV)
            ELSE
              L0AVL=L0ABSV+(IZ-1)*NX*NY
            ENDIF
            L0TEM=L0F(JTEM1)
            IF(ITURB.EQ.1) THEN
              L0KE=L0F(KE)
            ELSEIF(ITURB.EQ.2) THEN
              L0ENUT=L0F(VIST)
              L0EL1=L0F(LEN1)
            ENDIF
            VMIN=VARMIN(JPPDR); VMAX=VARMAX(JPPDR)
            DO I=1,NX*NY
              IF(.NOT.SLD(I)) THEN
                VABS=AMAX1(0.05,F(L0AVL+I))
                IF(ITURB.EQ.0) THEN   ! laminar
                  TINS=0.0
                ELSE
                  IF(ITURB.EQ.1) THEN ! KE solved
                    SQRTKE=SQRT(F(L0KE+I))
                  ELSE                ! ENUT & LEN1 stored
                    SQRTKE=F(L0ENUT+I)/(CMU*F(L0EL1+I))
                  ENDIF
                  TINS=100.0*SQRTKE/VABS
                ENDIF
                PPDR=(34.0-F(L0TEM+I))*((VABS-0.05)**0.62)*
     1                                        (0.37*VABS*TINS+3.14)
                F(L0PPDR+I)=AMAX1(VMIN,AMIN1(VMAX,PPDR))
              ELSE
                F(L0PPDR+I)=0.0
                TINS=0.0
              ENDIF
              IF(LTINS) F(L0TINS+I)=TINS
            ENDDO
          ENDIF
        ENDIF
      ELSEIF(ISC.EQ.7) THEN
C   * ------------------- SECTION 7 ---- Finish of sweep.
C...Set print control variable RPRINT for csv residual table
        RPRINT=.FALSE.
C...Resultant dry temperature TRES
        LOG1=.FALSE.
        IF(STEADY) THEN
C'####      LOG1=((MOD(ISWEEP+1,IFREQ).EQ.0.AND.CSG1.EQ.'SW') .OR.
          LOG1=((MOD(ISWEEP,IFREQ).EQ.0.AND.CSG1.EQ.'SW') .OR.
     1           MOD(ISWEEP,NPRINT).EQ.0) .OR.
     1           ENUFSW.OR.ISWEEP.EQ.LSWEEP
          RPRINT=((MOD(ISWEEP,IFREQ).EQ.0.AND.CSG1.EQ.'SW') .OR.
     1           MOD(ISWEEP,NPRINT).EQ.0) .OR.
     1           (MOD(ISWEEP,NPLT).EQ.0) .OR.
     1           ENUFSW.OR.ISWEEP.EQ.LSWEEP
        ELSE
          LOG1= (MOD(ISTEP,IFREQ).EQ.0 .OR.
     1           MOD(ISTEP,NTPRIN).EQ.0) .AND.
     1          (ENUFSW.OR.ISWEEP.EQ.LSWEEP)
          RPRINT=(MOD(ISTEP,IFREQ).EQ.0 .OR.
     1            MOD(ISTEP,NTPRIN).EQ.0) .AND.
     1           (ENUFSW.OR.ISWEEP.EQ.LSWEEP)
        ENDIF
        IF(NULLPR) RPRINT=.FALSE.
        IF(LTRES.AND.LOG1) THEN
          DO IZZ=1,NZ
            IF(LTRAD) THEN
              L0TRD=L0F(ANYZ(JT3,IZZ))
            ENDIF
            CALL SUB2(L0TEM,L0F(ANYZ(JTEM1,IZZ)),
     1                  L0TRES,L0F(ANYZ(JTRES,IZZ)))
            IF(LABSV) THEN
              L0AVL=L0F(ANYZ(JABSV,IZZ))
            ELSE
              L0AVL=L0ABSV+(IZZ-1)*NX*NY
            ENDIF
            DO IX=1,NX
              DO IY=1,NY
                I=IY+NY*(IX-1)
                IF(LTRAD) THEN
                  TRA=F(L0TRD+I)
                ELSE
                  TRA=TRAD
                ENDIF
C... Calculate 'dry resultant temperature'
                F(L0AVL+I)=ABS(F(L0AVL+I))
                F(L0TRES+I)=(TRA+F(L0TEM+I)
     1                 *SQRT(10.*F(L0AVL+I)))/(1.+SQRT(10.*F(L0AVL+I)))
              ENDDO
            ENDDO
          ENDDO
        ENDIF
C
        IF (LPAR) THEN
          CALL LGSUM_OR(LOG1)
          CALL LGSUM_OR(RPRINT)
        ENDIF
        IF(LOG1.OR.RPRINT) THEN
C... Calculate overall 'residence' time ie. 'free' volume
C    divided by the volumetric flow-rate. Also average inlet temperature
          FMASIN=0.0
          FMASIN2=0.0
          FVARIN=TINY
C... Calculate smoke, energy and momentum inflows for csv residual table
C... Calculate overall mass-inflow
          ILIM=ITWO(GD_NUMPAT,NUMPAT,LPAR)
          DO I=1,ILIM
            IF(LPAR) THEN
              IR=GD_INDPAT(I,1)
            ELSE
              IR=I
            ENDIF
            IF(LPAR) THEN
              CALL PGETCV(I,R1,GCO,GVAL)
            ELSE
              CALL GETCOV(NAMPAT(IR),R1,GCO,GVAL)
            ENDIF
            IF(QEQ(GVAL,-999.)) CYCLE
            IF(IR.LT.0) THEN
              SORCE=0.0
            ELSE
              CALL GETSO(IR,R1,SORCE)
            ENDIF
            IF(LPAR) CALL GLSUM(SORCE)
            IF(SORCE.GT.0.0) THEN ! mass inflow
              FMASIN=FMASIN+SORCE  ! sum mass
              DO IVAR=1,NPHI
                IF(SOLVE(IVAR)) THEN
                  IF(FLUID1(IVAR)) THEN
                    IF(IR.LT.0) THEN
                      SOR=0.0
                    ELSE
                      CALL GETSO(IR,IVAR,SOR)
                    ENDIF
                    IF(LPAR) CALL GLSUM(SOR)
                    FVARIN(IVAR)=FVARIN(IVAR)+ABS(SOR)
                  ENDIF
                ENDIF
              ENDDO
            ENDIF
            IF(.NOT.ONEPHS) THEN
              IF(IR.LT.0) THEN
                SORCE=0.0
              ELSE
                CALL GETSO(IR,R2,SORCE)
              ENDIF
              IF(LPAR) CALL GLSUM(SORCE)
              IF(SORCE.GT.0.0) THEN
                FMASIN2=FMASIN2+SORCE
                DO IVAR=1,NPHI
                  IF(SOLVE(IVAR)) THEN
                    IF(.NOT.FLUID1(IVAR)) THEN
                      IF(IR.LT.0) THEN
                        SOR=0.0
                      ELSE
                        CALL GETSO(IR,IVAR,SOR)
                      ENDIF
                      IF(LPAR) CALL GLSUM(SOR)
                      FVARIN(IVAR)=FVARIN(IVAR)+ABS(SOR)
                    ENDIF
                  ENDIF
                ENDDO
              ENDIF
            ENDIF
          ENDDO
C... Maximum momentum inflow taken as reference flux for csv residual table
          FMOMIN=AMAX1(FVARIN(3),FVARIN(5),FVARIN(7))
          IF(.NOT.ONEPHS) THEN
            FMOMIN2=AMAX1(FVARIN(4),FVARIN(6),FVARIN(8))
          ENDIF
C... Compute reference momentum outflow if reference inflow is zero
          IF(QLE(FMOMIN,1.E-10)) THEN
            DO IVAR=1,NPHI
              IF(FLUID1(IVAR)) FVARIN(IVAR)=TINY
            ENDDO
            FMASIN=0.0
            DO I=1,ILIM
              IF(LPAR) THEN
                IR=GD_INDPAT(I,1)
              ELSE
                IR=I
              ENDIF
              IF(LPAR) THEN
                CALL PGETCV(I,R1,GCO,GVAL)
              ELSE
                CALL GETCOV(NAMPAT(IR),R1,GCO,GVAL)
              ENDIF
              IF(QEQ(GCO,-999.0)) CYCLE
              IF(IR.LT.0) THEN
                SORCE=0.0
              ELSE
                CALL GETSO(IR,R1,SORCE)
              ENDIF
              IF(LPAR) CALL GLSUM(SORCE)
              IF(SORCE.LT.0) THEN ! mass outflow
                FMASIN=FMASIN+ABS(SORCE)  ! sum mass
                DO IVAR=1,NPHI
                  IF(SOLVE(IVAR)) THEN
                    IF(FLUID1(IVAR)) THEN
                      IF(IR.LT.0) THEN
                        SOR=0.0
                      ELSE
                        CALL GETSO(IR,IVAR,SOR)
                      ENDIF
                      IF(LPAR) CALL GLSUM(SOR)
                      FVARIN(IVAR)=FVARIN(IVAR)+ABS(SOR)
                    ENDIF
                  ENDIF
                ENDDO
              ENDIF
            ENDDO
            FMOMIN = AMAX1(FVARIN(3),FVARIN(5),FVARIN(7))
          ENDIF
          IF(.NOT.ONEPHS) THEN
            IF(QLE(FMOMIN2,1.E-10)) THEN
              DO IVAR=1,NPHI
                IF(.NOT.FLUID1(IVAR)) FVARIN(IVAR)=TINY
              ENDDO
              FMASIN2=0.0
              DO I=1,ILIM
                IF(LPAR) THEN
                  IR=GD_INDPAT(I,1)
                ELSE
                  IR=I
                ENDIF
                IF(LPAR) THEN
                  CALL PGETCV(I,R2,GCO,GVAL)
                ELSE
                  CALL GETCOV(NAMPAT(IR),R2,GCO,GVAL)
                ENDIF
                IF(QNE(GCO,-999.0)) THEN
                  IF(IR.LT..0) THEN
                    SORCE=0.0
                  ELSE
                    CALL GETSO(IR,R2,SORCE)
                  ENDIF
                  IF(LPAR) CALL GLSUM(SORCE)
                  IF(SORCE.LT.0) THEN ! mass outflow
                    FMASIN2=FMASIN2+ABS(SORCE)  ! sum mass
                    DO IVAR=1,NPHI
                      IF(SOLVE(IVAR)) THEN
                        IF(.NOT.FLUID1(IVAR)) THEN
                          IF(IR.LT.0) THEN
                            SOR=0.0
                          ELSE
                           CALL GETSO(IR,IVAR,SOR)
                          ENDIF
                          IF(LPAR) CALL GLSUM(SOR)
                          FVARIN(IVAR)=FVARIN(IVAR)+ABS(SOR)
                        ENDIF
                      ENDIF
                    ENDDO
                  ENDIF
                ENDIF
              ENDDO
              FMOMIN2= AMAX1(FVARIN(4),FVARIN(6),FVARIN(8))
            ENDIF
          ENDIF
          IF(JTEM1.GT.0) TBAR=FVARIN(JTEM1)
          IF(SOLVE(KE)) THEN
            CALL GETSOR('KESOURCE',KE,GKESO)
            IF(LPAR) CALL GLSUM(GKESO)
            FVARIN(KE)=AMAX1(FVARIN(KE),ABS(GKESO))
          ENDIF
          IF(SOLVE(EP)) THEN
            CALL GETSOR('KESOURCE',EP,GKESO)
            IF(LPAR) CALL GLSUM(GKESO)
            FVARIN(EP)=AMAX1(FVARIN(EP),ABS(GKESO))
          ENDIF
          IF(JAGE.GT.0) THEN
            IF(SOLVE(JAGE)) THEN
              CALL GETSOR('AGE',JAGE,GAGESO)
              IF(LPAR) CALL GLSUM(GAGESO)
              FVARIN(JAGE)=AMAX1(TINY,GAGESO)
            ENDIF
          ENDIF
          IF(NSMOK.GT.0.AND.JSMOK.GT.0) THEN
            FVARIN(JSMOK)=TINY
            IF(LPAR) THEN
              DO I=1,NSMOK
                IPG=NINT(F(L0SMK+I))
                IP=GD_INDPAT(IPG,1)
                IF(IP.LT.0) THEN
                  SSOR=0.0
                ELSE
                  CALL GETSO(IP,JSMOK,SSOR)
                ENDIF
                CALL GLSUM(SSOR)
                FVARIN(JSMOK)=FVARIN(JSMOK)+SSOR
              ENDDO
            ELSE
              DO I=1,NSMOK
                IP=NINT(F(L0SMK+I))
                CALL GETSO(IP,JSMOK,SSOR)
                FVARIN(JSMOK)=FVARIN(JSMOK)+SSOR
              ENDDO
            ENDIF
          ENDIF
C
          RHOIN=1.189 ! default inlet density
          IF(CP1.GT.0.AND.JTEM1.GT.0) THEN ! specific heat is constant
            TBAR=TBAR/(CP1*FMASIN+TINY) ! Tbar = Qin/(Cp*massin) - in K
            IF(QEQ(RHO1,GRND5))THEN ! check for Ideal Gas Law
              RHOIN=PRESS0*RHO1B/(TBAR+TINY) ! get average inlet density
            ENDIF
          ENDIF
C          IF (MIMD.AND.NPROC.GT.1) CALL DGSUM2(FMASIN,VOLSUM)
          IF(LPAR) CALL GLSUM(VOLSUM)
          VFLRT=AMAX1(TINY,FMASIN/RHOIN)
          RESTIME=VOLSUM/(VFLRT+TINY)
          IF(LOG1.AND.MYID.EQ.0.AND..NOT.NULLPR) THEN
            CALL WRITST
            IF(STEADY) THEN
              CALL WRIT1I('Sweep  ',ISWEEP)
            ELSE
              CALL WRIT2I('Sweep  ',ISWEEP,'; Step ', ISTEP)
            ENDIF
            CALL WRITBL
            CALL WRIT40('Overall residence time calculated as    ')
            CALL WRIT40('free volume/volum.flow-rate (in seconds)')
            IF(.NOT.STEADY) THEN
              WRITE(LUPR1,*)'for the time period ending at ',
     1                TIM,' seconds'
            ENDIF
            CALL WRITBL
            CALL WRIT1R('RES.TIME',RESTIME)
            CALL WRITBL
            CALL WRIT40('Ventilation rate in air changes per hour')
            CALL WRIT1R('ACH     ',3600./RESTIME)
            CALL WRITBL
            CALL WRIT40('The total free volume in the domain is (m^3)')
            CALL WRIT1R('VOLUME  ',VOLSUM)
            CALL WRITBL
            CALL WRIT40
     1             ('The total volumetric flow rate is (m^3/s, m^3/hr)')
            CALL WRIT2R('VFLRT s ',VFLRT,';VFLRT h',VFLRT*3600.)
            CALL WRITBL
            CALL WRITST
          ENDIF
C... Print normalised residuals to csv residual table
        ENDIF
        IF(RPRINT) THEN
          ISTO=0
          DO IVAR=1,NPHI
            IF(SOLVE(IVAR)) THEN
              ISTO=ISTO+1
              RNORM=ABS(FVARIN(IVAR))
              IF(IVAR.GE.3.AND.IVAR.LE.8) THEN
                IF(FLUID1(IVAR)) THEN
                  RNORM=FMOMIN+TINY
                ELSE
                  RNORM=FMOMIN2+TINY
                ENDIF
              ENDIF
              RESR=RESREF(IVAR); TOTR=TOTRES(IVAR)
              IF(LPAR) TOTR=TOTR*NPROC
              RESID(ISTO)=100.*RESR*TOTR/(RNORM+TINY)
            ENDIF
          ENDDO
C
          IF(MYID.EQ.0) THEN
            IF(STEADY) THEN
C Steady residual table
              WRITE(BUFF,'(I6,150('','',1PE11.3))')
     1                                       ISWEEP,(RESID(I),I=1,ISTO)
            ELSE
C Transient residual table
              WRITE(BUFF,'(I6,'','',1PE11.3,'','',I6,
     1        150('','',1PE11.3))') ISTEP,TIM,ISWEEP,(RESID(I),I=1,ISTO)
            ENDIF
            CALL MON_TABL(2,LUT,CTABLE,.FALSE.,BUFF)
C Write Reference inflow fluxes to the RESULT file
            CALL WRITBL
            IF(STEADY) THEN
              WRITE(LUPR1,1978)ISWEEP
            ELSE
              WRITE(LUPR1,1973)ISTEP,TIM,ISWEEP
            ENDIF
            IF(ONEPHS) THEN
              WRITE(LUPR1,'('' P1       '',1PE11.3,'' kg/s'',/,
     1                    '' U1,V1,W1 '',1PE11.3,'' N'')')FMASIN,FMOMIN
            ELSE
              WRITE(LUPR1,'('' R1       '',1PE11.3,'' kg/s'',/,
     1                    '' U1,V1,W1 '',1pe11.3,'' N'')')FMASIN,FMOMIN
              WRITE(LUPR1,'('' R2       '',1pe11.3,'' kg/s'',/,
     1                   '' U2,V2,W2 '',1pe11.3,'' N'')')FMASIN2,FMOMIN2
            ENDIF
            DO IVAR=12,NPHI
              IF(SOLVE(IVAR)) THEN
                IF(IVAR.EQ.JTEM1) THEN
                  WRITE(LUPR1,'(1X,A9,1PE11.3,'' W'')') NAME(IVAR),
     1                                                      FVARIN(IVAR)
                ELSEIF(IVAR.EQ.JSMOK) THEN
                  WRITE(LUPR1,'(1X,A9,1PE11.3,'' kg/s'')') NAME(IVAR),
     1                                                      FVARIN(IVAR)
                ELSE
                  WRITE(LUPR1,'(1X,A9,1PE11.3)') NAME(IVAR),FVARIN(IVAR)
                ENDIF
              ENDIF
            ENDDO
            CALL WRITBL
            CALL WRITST
          ENDIF
        ENDIF
1971    FORMAT(1X,I5,',',1PE11.3,',',I5,6(',',1PE11.3))
1972    FORMAT(1X,I5,',',1PE11.3,',',I5,5(',',1PE11.3))
1973    FORMAT(1X,'Inflow fluxes for normalised residual table at',
     1    ' ISTEP = ',I5,/' Time = ',1PE11.3,'s and ISWEEP = ',I5)
1974    FORMAT(1X,' P1',9x,'= ',1PE11.3,' kg/s'/1X,' U1, V1, W1 = ',
     1    1PE11.3,' N '/1X,' TEM1',7X,'= ',1PE11.3,' W')
1975    FORMAT(1X,' SMOK',7X,'= ',1PE11.3,' kg/s')
1976    FORMAT(1X,I5,6(',',1PE11.3))
1977    FORMAT(1X,I5,5(',',1PE11.3))
1978    FORMAT(1X,'Inflow fluxes for normalised residual table at',
     1    ' ISWEEP = ',I5)
C
C... Sprinkler-head activation
C
        IF(.NOT.STEADY.AND.NSPRAY.GT.0.AND.LTLINK) THEN
          DO JSP=1,NSPRAY
            IPT=NINT(F(L0SPRAY+JSP))
            CALL GETPTC(NAMPAT(IPT),GTYPE,JXF,JXL,JYF,JYL,JZF,JZL,
     1                                                          JTF,JTL)
            IXS=(JXF+JXL)/2; IYS=(JYF+JYL)/2; IZS=(JZF+JZL)/2
            L0TEM1 =ANYZ(JTEM1,IZS)
            IF(LABSV) THEN
              L0AVL=ANYZ(JABSV,IZS)
            ELSE
              L0AVL=-(L0ABSV+(IZS-1)*NX*NY)
            ENDIF
            IJSP=(IXS-1)*NY+IYS
            L0TEM1=IABS(L0TEM1); L0AVL=IABS(L0AVL)
            F(L0SP_TG+JSP)=F(L0TEM1+IJSP)
            F(L0SP_ABV+JSP)=F(L0AVL+IJSP)
C.. compute link temperature
            F(L0SP_A+JSP)= DT*SQRT(F(L0SP_ABV+JSP))/F(L0SP_RTI+JSP)
            F(L0SP_TL+JSP)=(F(L0SP_TL0+JSP)+F(L0SP_A+JSP)*
     &                     F(L0SP_TG+JSP))/(1.+F(L0SP_A+JSP))
            DO IZZ=JZF,JZL ! store link temp into TLNK over patch
              L0TLINK=L0F(ANYZ(JTLINK,IZZ))
              DO IX=JXF,JXL
                DO IY=JYF,JYL
                  IJ=(IX-1)*NY+IY
                  F(L0TLINK+IJ)=F(L0SP_TL+JSP)
                ENDDO
              ENDDO
            ENDDO
          ENDDO
        ENDIF
      ELSEIF(ISC.EQ.8) THEN
C   * ------------------- SECTION 8 ---- Finish of time step.
C
C... For transient problems calculate percentage of air
C    changed during the time-step (divide tstep/res time?)
        LOG1=.FALSE.
        IF(.NOT.STEADY) THEN
          LOG1= (MOD(ISTEP,IFREQ).EQ.0 .OR.
     1           MOD(ISTEP,NTPRIN).EQ.0) .AND.
     1          (ENUFSW.OR.ISWEEP.EQ.LSWEEP-1.OR.ISWEEP.EQ.LSWEEP)
        ENDIF
        IF(LOG1.AND.MYID.EQ.0.AND..NOT.NULLPR) THEN
          AIRRT=DT/RESTIME*100
          CALL WRITBL
          CALL WRIT40('Percentage of air exchanged during this ')
          CALL WRIT40('time period was                         ')
          CALL WRIT1R('  AIR % ',AIRRT)
          CALL WRITBL
          CALL WRITST
        ENDIF
C... Sprinkler-head activation
C
        IF(.NOT.STEADY.AND.NSPRAY.GT.0) THEN
          CALL WRITST
          CALL WRIT1I('ISTEP ',ISTEP)
          CALL WRIT1R('TIME  ',TIM)
          LOG2=.FALSE.
          DO JSP=1,NSPRAY
            IPT=NINT(F(L0SPRAY+JSP))
            WRITE(14,*) ' Spray-head link data for: ', OBJNAM(IPT)
            SP_TAC = F(L0SP_TA+JSP)+TEMP0-273.
            SP_TLC = F(L0SP_TL+JSP)+TEMP0-273.
            SP_TGC = F(L0SP_TG+JSP)+TEMP0-273.
            WRITE(14,1985) 'activation temperature = ',SP_TAC,'degC'
            WRITE(14,'(1X,A25,1PE11.3,1X,A9)')
     1        'response time index    = ',F(L0SP_RTI+JSP),'(m.s)^0.5'
            WRITE(14,1985) 'link temperature       = ',SP_TLC,'degC'
            WRITE(14,1985) 'gas temperature        = ',SP_TGC,'degC '
            WRITE(14,1985) 'flow velocity          = ',F(L0SP_ABV+JSP),
     1                                               ' m/s '
            WRITE(14,1985) 'link coefficient       = ',F(L0SP_A+JSP)
            SP_TL0C=F(L0SP_TL0+JSP)+TEMP0-273.
            WRITE(14,1985) 'old link temperature   = ',SP_TL0C,'degC'
 1985       FORMAT(1X,A25,1PE11.3,1X,A4)
            IF(F(L0SP_TL+JSP).GT.F(L0SP_TA+JSP))  THEN
              WRITE(14,*) '!!!! ',OBJNAM(IPT),' active !!!!'
              ISP=NINT(F(L0SPRAYT+JSP))
              IF(EQZ(F(L0ACTIV+ISP))) LOG2=.TRUE. ! a spray has been activated
              F(L0ACTIV+ISP)=1.                   ! set active spray flag
            ENDIF
            CALL WRITBL
C... save current value into old store for next step
            F(L0SP_TL0+JSP)=F(L0SP_TL+JSP)
!C... also save into patch-wise store, so it will be written to PHI
!            L0PAT= L0PVAR(30,IPT, JZF)
!            F(L0PAT+1)=F(L0SP_TL+JSP)
          ENDDO
C... a spray was activated so modify SPIN file for GENTRA
          IF(LOG2) CALL ACTIVATE_SPRAYS
C
          NTB=20 ! Number of sprays per TLINK file
          IF(ISTEP.EQ.FSTEP) THEN
            LSTART=.TRUE.
            IS1=1;IS2=MIN(NSPRAY,NTB);ITAB=0
            DO WHILE (IS1.LE.NSPRAY)
              ITAB=ITAB+1; FORM='tlinkiiiiii.csv'
              WRITE(FORM(6:11),'(I6.0)') ITAB
              LL=LENGZZ(FORM);CALL REMSPC(FORM,LL)
              CALL MON_TABL(1,LUT,FORM,LSTART,' ') ! Open
              IS1=IS1+NTB; IS2=MIN(NSPRAY,IS2+NTB)
            ENDDO
            IF(LSTART) THEN
              IS1=1;IS2=MIN(NSPRAY,NTB);ITAB=0
              DO WHILE (IS1.LE.NSPRAY)
                FORM='(1X,"TIME",8X,iiiiii(",",A13))'
                WRITE(FORM(15:20),'(I6.0)') IS2-IS1+1
                LL=LENGZZ(FORM);CALL REMSPC(FORM,LL)
                WRITE(BUFF,FORM)(OBJNAM(NINT(F(L0SPRAY+IP))),IP=IS1,IS2)
                ITAB=ITAB+1; FORM='tlinkiiiiii.csv'
                WRITE(FORM(6:11),'(I6)') ITAB
                LL=LENGZZ(FORM);CALL REMSPC(FORM,LL)
                CALL MON_TABL(2,LUT,FORM,.FALSE.,BUFF) ! Write header
                IS1=IS1+NTB; IS2=MIN(NSPRAY,IS2+NTB)
              ENDDO
            ENDIF
          ENDIF
          IS1=1;IS2=MIN(NSPRAY,NTB);ITAB=0
          DO WHILE (IS1.LE.NSPRAY)
            FORM='(1PE13.6,iiiiii(",",1PE13.6))'
            WRITE(FORM(10:15),'(I6.0)') IS2-IS1+1
            LL=LENGZZ(FORM);CALL REMSPC(FORM,LL)
            WRITE(BUFF,FORM) TIM,(F(L0SP_TL+JSP),JSP=IS1,IS2)
            ITAB=ITAB+1; FORM='tlinkiiiiii.csv'
            WRITE(FORM(6:11),'(I6)') ITAB
            LL=LENGZZ(FORM);CALL REMSPC(FORM,LL)
            CALL MON_TABL(2,LULOCAL,FORM,.FALSE.,BUFF) ! Write data
            IS1=IS1+NTB; IS2=MIN(NSPRAY,IS2+NTB)
          ENDDO
        ENDIF
        CALL WRITST
C.... Write fire heat and mass sources to file
        IF(LOG1) THEN
          IF(LPAR) THEN
            DO I=1,NFIR
              IPG=NINT(F(L0FIR+I))
              IP=GD_INDPAT(IPG,1)
              IF(IP.LT.0) THEN
                FSOR=0.0
              ELSE
                CALL GETSO(IP,JTEM1,FSOR)
              ENDIF
              CALL GLSUM(FSOR)
              F(L0QSOR+I)=FSOR
            ENDDO
            IF(MYID.EQ.0.AND.NFIR.GT.0) THEN
              WRITE(BUFF,'(1PE10.3,99('','',1PE10.3))')
     1                                 TIM-0.5*DT,(F(L0QSOR+I),I=1,NFIR)
              CALL MON_TABL(2,LUT,QTABLE,.FALSE.,BUFF)
            ENDIF
            DO I=1,NSMOK
              IPG=NINT(F(L0SMK+I))
              IP=GD_INDPAT(IPG,1)
              IF(IP.LT.0) THEN
                SSOR=0.0
              ELSE
                CALL GETSO(IP,JSMOK,SSOR)
              ENDIF
              CALL GLSUM(SSOR)
              F(L0SSOR+I)=SSOR
            ENDDO
            IF(MYID.EQ.0.AND.NSMOK.GT.0) THEN
              WRITE(BUFF,'(1PE10.3,99('','',1PE10.3))')
     1                                TIM-0.5*DT,(F(L0SSOR+I),I=1,NSMOK)
              CALL MON_TABL(2,LUT,STABLE,.FALSE.,BUFF)
            ENDIF
          ELSE
            DO I=1,NFIR
              IP=NINT(F(L0FIR+I))
              CALL GETSO(IP,JTEM1,F(L0QSOR+I))
            ENDDO
            IF(NFIR.GT.0) THEN
              WRITE(BUFF,'(1PE10.3,99('','',1PE10.3))')
     1                                 TIM-0.5*DT,(F(L0QSOR+I),I=1,NFIR)
              CALL MON_TABL(2,LUT,QTABLE,.FALSE.,BUFF)
            ENDIF
            DO I=1,NSMOK
              IP=NINT(F(L0SMK+I))
              CALL GETSO(IP,JSMOK,F(L0SSOR+I))
            ENDDO
            IF(NSMOK.GT.0) THEN
              WRITE(BUFF,'(1PE10.3,99('','',1PE10.3))')
     1                                TIM-0.5*DT,(F(L0SSOR+I),I=1,NSMOK)
              CALL MON_TABL(2,LUT,STABLE,.FALSE.,BUFF)
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      RETURN
C... Error exit
  999 CALL WRIT40('Error in number denoting a TM patch.    ')
      CALL WRIT40('Check Q1 for TMx, where 0 < x <999      ')
      CALL SET_ERR(524, 'Error. See result file.',1)
c      CALL WAYOUT(1)
      END
C**********************************************************************
      FUNCTION PVH2O(TS)
C   *****************************************************************
C   * Purpose       : This function determines the saturation       *
C   *                 pressure of the water vapour in Pascals.      *
C   * Ref           : Arden-Buck Correlation.                       *
C   *                 Buck,A.L.,"New equations for computing vapor  *
C   *                 pressure and enhancement factor", J.Appl.     *
C   *                 Meteorol., Vol.20, p1527-1532,(1981).         *
C   * Validity range: -80 Deg C to 100 Deg C.                       *
C   * Author        : M.R.MALIN                                     *
C   * Date revised  : 19/01/09                                      *
C   * Q.A. Engineer : M. R. MALIN                                   *
C   *---------------------------------------------------------------*
C   * Input arguments :-                                            *
C   *   TS     -    Real , Saturation temperature (Deg C)           *
C   *****************************************************************
       T       = AMAX1(-80.0,TS)
       PVH2O = 611.21*EXP(T*(18.678-T/234.5)/(257.14+T))
       END
C**********************************************************************
      FUNCTION PVH2O_old(TS)
C   *****************************************************************
C   * Purpose      : This function determines the saturation        *
C   *                pressure of the water vapour.                  *
C   * Ref          : ASHRAE Correlation - 0 to 100 DEG c            *
C   * Author       : M.R.MALIN                                      *
C   * Date revised : 18/02/04                                       *
C   * Q.A. Engineer: M. R. MALIN                                    *
C   *---------------------------------------------------------------*
C   * Input arguements :-                                           *
C   *   TS     -    Real , Saturation temperature (Deg C)           *
C   *****************************************************************
      LOGICAL QLE
      IF(QLE(TS,-50.0)) THEN
        PVH2O=6.5
      ELSE
        T1  = AMIN1(TS,500.0)
C... Convert Tsat from degC to degK
        T   = T1+273.15
        TC  = 647.286
        C1  = -7.419242
        C2  = 2.9721E-1
        C3  = -1.155286E-1
        C4  = 8.685635E-3
        C5  = 1.094098E-3
        C6  = -4.39993E-3
        C7  = 2.520658E-3
        C8  = -5.218684E-4
        GA  = 0.01
        TP  = 338.15
        DT1 = GA*(T-TP)
        DT2 = DT1*DT1
        DT3 = DT2*DT1
        DT4 = DT3*DT1
        DT5 = DT4*DT1
        DT6 = DT5*DT1
        DT7 = DT6*DT1
        SUM = C1 + C2*DT1 + C3*DT2 + C4*DT3 + C5*DT4 + C6*DT5
     1           + C7*DT6 + C8*DT7
        TERM= SUM*(TC/T-1.0)
        PC    = 22.089E6
        PVH2O = EXP(TERM)*PC
      ENDIF
      pvh2o_old=pvh2o
      END
C**********************************************************************
      SUBROUTINE MON_TABL(IGO,LUT,TABLE,LSTART,BUFF)
      COMMON /LDATA/ LDAT1(31),NULLPR,LDAT2(52)
      LOGICAL LDAT1,NULLPR,LDAT2
      LOGICAL LSTART
      CHARACTER*(*) TABLE,BUFF,BUFF2(5)*80
  111 CONTINUE
      IF(IGO.EQ.1) THEN
        IF(LSTART) THEN ! create new file at start of run
          OPEN(LUT,FILE=TABLE,STATUS='NEW',IOSTAT=IOS)
          IF(IOS.NE.0) THEN ! file existed, so delete and start again
            OPEN(LUT,FILE=TABLE,STATUS='UNKNOWN',IOSTAT=IOS)
            CLOSE(LUT,STATUS='DELETE',IOSTAT=IOS)
            IF(IOS.NE.0) GO TO 700
            OPEN(LUT,FILE=TABLE,STATUS='NEW',IOSTAT=IOS)
          ENDIF
        ELSE ! restart run - try to use existing file
          OPEN(LUT,FILE=TABLE,STATUS='OLD',IOSTAT=IOS)
          IF(IOS.NE.0) THEN ! old file does not exist - make new one
            OPEN(LUT,FILE=TABLE,STATUS='UNKNOWN',IOSTAT=IOS)
            CLOSE(LUT,STATUS='DELETE',IOSTAT=IOS)
            IF(IOS.NE.0) GO TO 700
            OPEN(LUT,FILE=TABLE,STATUS='NEW',IOSTAT=IOS)
            LSTART=.TRUE.
          ENDIF
        ENDIF
      ELSE
        LL=LENGZZ(BUFF)
        IF(BUFF(LL:LL).EQ.',') LL=LL-1 ! strip off any trailing ,
        OPEN(LUT,FILE=TABLE,STATUS='OLD',ACCESS='APPEND',IOSTAT=IOS)
        IF(IOS.EQ.0) THEN ! file exists
          WRITE(LUT,'(1X,A)',IOSTAT=IOS) BUFF(1:LL)
          IF(IOS.NE.0) GO TO 700
        ELSE              ! file does not exist, create one
          OPEN(LUT,FILE=TABLE,STATUS='UNKNOWN',IOSTAT=IOS)
          CLOSE(LUT,STATUS='DELETE',IOSTAT=IOS)
          IF(IOS.NE.0) GO TO 700
          OPEN(LUT,FILE=TABLE,STATUS='NEW',IOSTAT=IOS)
          IF(IOS.GT.0) THEN
            CALL WRIT40('Failed to open file '//TABLE)
            CALL WRIT40(BUFF(1:LL))
          ELSE
            WRITE(LUT,'(1X,A)',IOSTAT=IOS) BUFF(1:LL)
            IF(IOS.NE.0) GO TO 700
          ENDIF
        ENDIF
      ENDIF
      CLOSE(LUT,STATUS='KEEP',IOSTAT=IOS)
      RETURN
  700 continue
      BUFF2(1)='Error opening or writing file '//TABLE
      BUFF2(3)='Is it already opened in another program?'
      BUFF2(4)=
     1 'If so, close the other program and click OK to try again,'
      BUFF2(5)='or Cancel to abort the run'
      CALL IOEMZZ(IOS,BUFF2(2)) ! get system error message
      CALL SET_MOUSE_FLAGS
      IACT=1 ! initialise for NULLPR=T
      IF(.NOT.NULLPR) CALL ERRMSG(BUFF2,5,222,IACT) ! display error dialog
      IF(IACT.EQ.0 )THEN
        CLOSE(LUT,IOSTAT=IOS)
        GO TO 111 ! go to top and try again
      ELSE
        CALL SET_ERR(548, 'Program terminated by user.',1)
      ENDIF
      END
C**********************************************************************
      SUBROUTINE ACTIVATE_SPRAYS
C.... This subroutine manipulates the GENTRA SPIN input file in order to reset
C     track start times if the spray activation criterion has been satisfied.
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/tracmn'
      INCLUDE '/phoenics/d_includ/pltcfile'
      COMMON/SPRAY1/ NSPRAY,L0SPRAY,L0SP_TA,L0SP_RTI,L0SP_TL,L0SP_TL0,
     1     L0SP_TG,L0SP_ABV,L0SP_A,ITURB,NSPRAYT,L0SPRAYT,L0ACTIV
      CHARACTER LINE*200,WD*20,BUFF(6)*80
      DIMENSION WD(20),NCH(20),TEMPRL(20)
      LOGICAL ERROR,QEQ
C
      LUINP=LUNIT(37); LUOUT=137
      NMFIL(37)=GINFIL
      CALL OPENZZ(37,IERR) ! open SPIN file
      IF(IERR.NE.0) GO TO 1000 ! problem opening it
      OPEN(LUOUT,FILE='spin_tmp',STATUS='UNKNOWN',IOSTAT=IOS) ! open scratch
      IF(IOS.NE.0) GO TO 1000 ! problem opening it
      CLOSE(LUOUT,STATUS='DELETE',IOSTAT=IOS) ! delete it
      IF(IOS.NE.0) GO TO 1000 ! problem opening it
      OPEN(LUOUT,FILE='spin_tmp',STATUS='NEW',IOSTAT=IOS) ! open new copy
      IF(IOS.NE.0) GO TO 1000 ! problem opening it
      DO ISP=1,NSPRAYT ! loop over all sprays
        LINE=' '; IPORT=0
 434    READ(LUINP,'(A)',END=431) LINE
        IF(LINE(1:).EQ.' '.OR.INDEX(LINE,'*').NE.0) THEN ! found header or blank
          LINLEN=LENGZZ(LINE)
          IF(INDEX(LINE,'ports').NE.0) THEN ! found number of ports
            CALL SPLTZZ(LINE, WD, NWDS, NCH, LINLEN)
            NPORTS=RRDGGG(WD(NWDS),NCH(NWDS),ERROR) ! read from last item on lin
            IF(ERROR) GO TO 1001
          ELSEIF(INDEX(LINE,'size').NE.0) THEN ! found number of size ranges
            CALL SPLTZZ(LINE, WD, NWDS, NCH, LINLEN)
            NSIZES=RRDGGG(WD(NWDS),NCH(NWDS),ERROR) ! read from last item on lin
            IF(ERROR) GO TO 1001
          ENDIF
          WRITE(LUOUT,'(A)') LINE(1:LINLEN) ! copy to scratch
          GOTO 434
        ENDIF
        BACKSPACE LUINP ! have reached start of data
        DO IPORT=1, NPORTS*NSIZES ! loop over all droplets for this spray
 435      READ(LUINP,'(A)',END=431) LINE
          IF(LINE(1:).EQ.' '.OR.INDEX(LINE,'*').NE.0) GOTO 435 ! skip comment or
          LINLN1=LENGZZ(LINE)
          IF(NINT(F(L0ACTIV+ISP)).EQ.1) THEN ! current spray is active
            CALL UCASZZ (LINE,LINLN1)
            LINLEN=LENGZZ(LINE)
            CALL SPLTZZ(LINE, WD, NWDS, NCH, LINLEN)
            TSTRT=RRDGGG(WD(NWDS-1),NCH(NWDS-1),ERROR) ! get start time for inje
            IF(ERROR) GO TO 1002
            IF(QEQ(TSTRT,1.0E+10)) THEN ! 1.E10 flags auto-start
              WRITE(WD(NWDS-1),'(1PE9.3)') TIM-1.001*DT ! overwrite start with c
              TEND=RRDGGG(WD(NWDS),NCH(NWDS),ERROR)
              IF(ERROR) GO TO 1003
              WRITE(WD(NWDS),'(1PE9.3)') TEND+TIM ! overwrite end with duration
            ENDIF
            WRITE(LUOUT,'(2X,12(A11))') (WD(IWD),IWD=1,NWDS) ! copy to scratch
          ELSE
            WRITE(LUOUT,'(A)') LINE(1:LINLN1)
          ENDIF
        ENDDO
      ENDDO
      CLOSE(LUINP,STATUS='DELETE',IOSTAT=IOS1) ! delete existing SPIN file
      CLOSE(LUOUT,STATUS='KEEP',IOSTAT=IOS2)   ! close scratch file
      IF(IOS1+IOS2.NE.0) GO TO 1000            ! problem closing or deleting
      CALL COPY_FILE('spin_tmp',GINFIL,1)      ! copy scratch to SPIN
      CALL DEL_FILE('spin_tmp',1)              ! delete scratch
      IF(IOS.NE.0) GO TO 1000
      RETURN
  431 CONTINUE
      CALL SET_ERR(576,'Error. Incomplete GENTRA SPIN file',1)
 1000 CONTINUE
      CALL SET_ERR(575,'Error. Cannot rewrite GENTRA SPIN file',1)
 1001 CONTINUE
      CALL SET_ERR(577,'Error. Cannot read NPORTS from SPIN file',1)
 1002 CONTINUE
      CALL SET_ERR(578,'Error. Cannot read start time from SPIN file',1)
 1003 CONTINUE
      CALL SET_ERR(579,'Error. Cannot read end time from SPIN file',1)
      END
c