c 

      SUBROUTINE FLARGR
      INCLUDE 'objnam'
      INCLUDE 'patcmn'
      INCLUDE 'farray'
      USE human_thermal_comfort_pet_mod
      INCLUDE 'd_earth/parvar'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      INCLUDE 'grdbfc'
      INCLUDE 'gracm2'
      COMMON/F0/IF01(24),KZZW,IF05(280)
      COMMON /INDAUX/L0ISL,L0IST,L0SL,L0ST,NSTO,NSOL,L0SLRS,L0TTRS,
     1        L0TTFL,L0TTLS,L0MAXC,L0MAXV,L0MINV,L0RATE,L0MAXI,L0NETT
      INCLUDE 'parear'
      PARAMETER (NLG=20, NIG=20, NRG=100, NCG=10)
      COMMON/LGRND/LG(NLG)/IGRND/IG(NIG)/RGRND/RG(NRG)/CGRND/CG(NCG)
      COMMON /FACES/L0FACE,L0FACZ
      COMMON /MDPTLN/NMDPTE,L0PTLN,LOBJBX
      COMMON /L0VRPAT/ L0IVRP
      COMMON /GEODMN0/ I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
     1                I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
      LOGICAL LG,  QNE, QEQ, QLE
      CHARACTER*4 CG
      CHARACTER*30 FORM
      LOGICAL ISINZZ, LOG1,LOG2, SLD, RPRINT,GRN
      COMMON/GENI/NXNY,IG1(8),NFM,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,LTAPP,LTAPR,ALLCELLS
      LOGICAL G_RESTRT,INPARDOM,LSTART,EQZ,LPAR,LTLINK,ISACTIVECELLS
      LOGICAL LHRAT, LMH2O, LRELH, LPSAT, LXH2O, LPVAP, LOPTD, FLUID1
      LOGICAL LUTCI,LPET,LTSIB
      COMMON/GSLIPL/C1_V,DFMODL,GSLIPX,GSLIPY,GSLIPZ,S_UP,S_DOWN,S_VERT
      LOGICAL C1_V,DFMODL,GSLIPX,GSLIPY,GSLIPZ,S_UP,S_DOWN,S_VERT
      LOGICAL LSLDIR,NEZ
      REAL KEXT
C... Regression parameters for Productivity loss
      REAL BB(2,7),FVARIN(150),RESID(150),VALS(20),MOUT
      INTEGER IDUM(NUMPAT)
      INTEGER, ALLOCATABLE :: IRMPAT(:),IRMSOR(:)
      REAL, ALLOCATABLE :: HVABS(:)
      CHARACTER*4 NAM(150), BUFF*1024
      CHARACTER*12 COBNAM,PATNM*8
      CHARACTER*256 CTABLE, QTABLE, STABLE,CH13*13,CHGR_13_R*13
      DOUBLE PRECISION UTCI_APPROX
C... SAVE variables required each visit
      SAVE CTABLE, QTABLE, STABLE
      SAVE JPMV,JPPD,JTCL,L0ABSV,L0TCLC,L0TINS,JTAPP,JTAPR,JLIT,JRM
      SAVE CL,RMET,WORK,RMW,FCL,RH,TRAD,QRAD
      SAVE ACON,BCON,CCON,ACN2,BCN2,CCN2,ROX,YCO,YS,KEXT
      SAVE LPPD,LPMV,LPPDR,LPLOS,LTINS,LTAPP,LTAPR
      SAVE LABSV,LTRES,VOLSUM,RESTIME,IFREQ,LTRAD,LSDEN,
     1     LSLEN,LPPM,LSMOK,LSLN2,LPPPM,LCOPM,LOPTD,ALLCELLS,LUTCI,LPET,
     1     LTSIB
      SAVE JTEM1,JTRES,JABSV,JT3,JSDEN,JSLEN,JPPM,JSMOK,JSLN2,JPPPM,
     1                    JCOPM,JOPTD,JAGE,JAGE2,JTLINK,JUTCI,JPET,JTSIB
      SAVE JHRAT, JRELH, JPSAT, JMH2O, JXH2O, JPVAP, JPPDR, JPLOS, JTINS
      SAVE LHRAT, LMH2O, LRELH, LPSAT, LXH2O, LPVAP, LTLINK, JAEE,JACE
      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,NROOM
      SAVE NSMOK,L0SMK,L0SSOR
      SAVE IRMPAT,IRMSOR
      SAVE HVABS,AGESUM
!C### MRM 05.06.15 Aerosol dispersion/deposition model
!      SAVE DFMODL
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
      IF(USP) RETURN
      IXL=IABS(IXL)
C... Call HOTBOX to calculate fan operating point or system curve
      CALL HTBXGR
      NAMSUB='FLARGR'
      IF(IGR==1) GO TO 1
      IF(IGR==8) GO TO 8
      IF(IGR==13) GO TO 13
      IF(IGR==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: 301122')
        CALL WRIT40('Call to spec. Ground FLARGR.F of: 301122')
      ENDIF
C ... Initialise
      OBJTYP(5)='OPENING'
      LPAR=MIMD.AND.NPROC>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)=='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>0) THEN
C... there are some sprays with active link temperatures.
          JTLINK=LBNAME('TLNK'); LTLINK=STORE(JTLINK)
          IF(MYID==0.AND..NOT.LTLINK.AND..NOT.NULLPR) 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/=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)=='SPRAY') THEN
              CALL MAKEPV(30,IP)
            ENDIF
          ENDDO
        ENDIF
      ELSE
        NSPRAY=0; NSPRAYT=0
        G_RESTRT=.FALSE.
      ENDIF
C ... Aerosol dispersion & deposition model - initialisation I
      DFMODL=.FALSE.
      CALL GETSDL('DFLUX','DFMODL',DFMODL)
      IF(DFMODL) CALL GXDFM
C
      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
      JUTCI =LBNAME('UTCI') ! Universal Thermal Climate Index
      JPET  =LBNAME('PET')  ! Physiological Equivalent Temperature
      JTSIB =LBNAME('TSIB') ! Thermal Sensetion Index (BEAM)
      LHRAT =STORE(JHRAT)
      LRELH =STORE(JRELH)
      LPSAT =STORE(JPSAT)
      LMH2O =STORE(JMH2O)
      LXH2O =STORE(JXH2O)
      LPVAP =STORE(JPVAP)
      LPPDR =STORE(JPPDR)
      LTINS =STORE(JTINS)
      LUTCI =STORE(JUTCI)
      LPET  =STORE(JPET)
      LTSIB =STORE(JTSIB)
C
      JTEM1=LBNAME('TEM1') ! temperature
      JTRES=LBNAME('TRES') ! dry resultant temperature
      JT3  =LBNAME('T3')   ! radiant temperature
      JABSV=LBNAME('VABS') ! 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
      JAGE2 =LBNAME('AGE2') ! Age of air in each room
      JAEE  =LBNAME('AEE')  ! Air Exchange Effectiveness
      JACE  =LBNAME('ACE')  ! Air Exchange Effectiveness
      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(JAGE>0.AND.DEN1<=0) THEN
        CALL WRIT40('Calculation of Air Exchange Effectiveness')
        CALL WRIT40('requires storage of density (STORE(DEN1)in Q1)')
        CALL WRYT40('Calculation of Air Exchange Effectiveness')
        CALL WRYT40('requires storage of density (STORE(DEN1)in Q1)')
        CALL SET_ERR(524, 'Error. See result file.',1)
      ENDIF
      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<=0) BCON=YS*KEXT/(1.0+ROX)
        IF(BCN2<=0) BCN2=YS*KEXT/(1.0+ROX)
      ELSE
        LSLEN=.FALSE.; LSLN2=.FALSE.; LPPM=.FALSE.
        LPPPM=.FALSE.; LCOPM=.FALSE.; LOPTD =.FALSE.
      ENDIF
C
      IFREQ=IDISPA
      IF(IFREQ==0) IFREQ=LSWEEP+1
C
      JPMV=LBNAME('PMV')  ! Predicted Mean Vote
      LPMV=STORE(JPMV)
      JTAPP=LBNAME('TAPP')  ! Apparent Temperature
      LTAPP=STORE(JTAPP)
      JTAPR=LBNAME('TAPR')  ! Apparent Temperature with radiation
      LTAPR=STORE(JTAPR)
      IF(LPMV) THEN
        JPPD=LBNAME('PPD')  ! Predicted Percentage Dissatified
        LPPD=JPPD>0
        JTCL=LBNAME('TCL')  ! Clothing temperature in deg C
        IF(JTCL==0) THEN
C... If clothing temperature not STOREd, make 3D store
          CALL GXMAKE0(L0TCLC,NXNY*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<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<0.078) THEN
          FCL=1.0+1.29*CL
        ELSE
          FCL=1.05+0.645*CL
        ENDIF
        RMET=-999.
        CALL GETSDR('FLAIR','METABOLIC',RMET)
        IF(RMET<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<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<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
      ENDIF
      IF(LPMV.OR.LUTCI.OR.LTAPR.OR.LTAPP) THEN
        RH=-999.
        CALL GETSDR('FLAIR','RELHUMID',RH)
        IF((RH<0.0.OR.RH>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.OR.LUTCI.OR.LPET) 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
        IF(LTAPR.OR.LTSIB) THEN
          QRAD=75.0; CALL GETSDR('FLAIR','TAPR-QRAD',QRAD)
          JLIT=LBNAME('LIT')
          ALLCELLS=.FALSE.; CALL GETSDL('SUNLIGHT','ALLCELLS',ALLCELLS)
        ENDIF
C... Draught rating
      IF(LPPDR.OR.LTINS) THEN
        IF(.NOT.LTINS) CALL GXMAKE0(L0TINS,NXNY,'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>0) THEN
          WRITE(LUPR1,*) 'Link temperature calculation',
     1                                           ' has been activated '
        ENDIF
      ENDIF
C
C... Marker for rooms
      JRM=LBNAME('ROOM')
C... Find mass source patches inside rooms
      IF(JRM>0) THEN
        NROOM=0
        ILIM=ITWO(GD_NUMPAT,NUMPAT,LPAR)
        DO IPG=1,ILIM ! loop over patches looking for rooms
          IF(LPAR) THEN
            PATNM=GD_NAMPAT(IPG)(2:)
          ELSE
            PATNM=NAMPAT(IPG)
          ENDIF
          IF(PATNM(1:4)=='ROOM') THEN ! found a room
            NROOM=NROOM+1; IDUM(NROOM)=IPG ! store its patch number
          ENDIF
        ENDDO ! end loop over patches
        IF(NROOM>0) THEN ! rooms exist
          ALLOCATE(IRMPAT(NROOM),STAT=IERR) ! stores patch number of room IRM
          ALLOCATE(HVABS(NROOM),STAT=IERR)  ! stores height above floor for average velocity
          IRMPAT=IDUM(1:NROOM); IDUM=0; I0=1 ! copy from IDUM, set in loop over patches above
          HVABS(1:NROOM)=1.2 ! default height for area-averaged VABS is 1.2m
          DO IRM=1,NROOM    ! loop over rooms to find internal sources
            IPG=IRMPAT(IRM) ! get the global patch number for this room
            IF(LPAR) THEN
              IPR=GD_INDPAT(IPG,1)   !  Local Index of Patch
            ELSE
              IPR=IPG
            ENDIF
            IF(IPR<0) CYCLE ! not on this processor
            NSO=0 ! initialise to no internal sources
            CALL GETPAT(IPR,IR,TYP,IX1,IX2,IY1,IY2,IZ1,IZ2,IT1,IT2) ! get limits
            CALL GETSDR(OBJNAM(IPR),'HVABS',HVABS(IRM)) ! get VABS height for room
            DO IS=1,NUMPAT  ! now search for mass sources inside the rrom
              IF(IS==IPR) CYCLE
              CALL GETCV(IS,P1,GCO,GVAL) ! get CO & VAL for P1
              IF(GCO>0.0.OR.GRN(GCO)) THEN ! there is a source
                CALL GETPAT(IS,IR,TYP,I1,I2,J1,J2,K1,K2,IT1,IT2) ! get limits
                IF(I1>=IX1.AND.I2<=IX2.AND.J1>=IY1.AND.J2<=IY2.AND.
     1                                     K1>=IZ1.AND.K2<=IZ2) THEN
                  NSO=NSO+1; IDUM(NSO+I0)=IS ! increment counter and save patch number
                ENDIF
              ENDIF
            ENDDO
            IDUM(I0)=NSO; I0=I0+NSO+1 ! save number of sources and increment offset
          ENDDO ! end loop over rooms
          ALLOCATE(IRMSOR(I0),STAT=IERR) ! save source patch number for all rooms
          IRMSOR=IDUM(1:I0) ! copy from IDUM set in loop above
        ENDIF
      ENDIF
C... If absolute velocity not STOREd, make 3D store
      IF(.NOT.LABSV.AND.(LTRES.OR.LPMV.OR.NSPRAY>0.OR.LPPDR.OR.
     1                LTAPP.OR.LTAPR.OR.LTINS.OR.NROOM>0.OR.LUTCI)) THEN
        CALL GXMAKE0(L0ABSV,NXNY*NZ,'ABSV')
      ENDIF
C... Normalised Residual Tabular Printout
      LUT=87; CTABLE='conv_table.csv'
C.. Open csv table file and write headings
      IF(STEADY) THEN
        LSTART=FSWEEP==1
      ELSE
        LSTART=FSTEP==1
      ENDIF
      ISTO=0
      DO II=1,NPHI
        IF(SOLVE(II).AND.NAME(II)/='LTLS') THEN
          ISTO=ISTO+1; NAM(ISTO)=NAME(II)
        ENDIF
      ENDDO
      IF(MYID==0) CALL MON_TABL(1,LUT,CTABLE,LSTART,' ') ! Open
      IF(MYID==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)=='FIR') THEN
            LL=LENGZZ(PATNM)
            IF(PATNM(LL:LL)=='A') NFIR0=NFIR0+1
            IF(PATNM(LL:LL)/='A') NSMOK0=NSMOK0+1
          ENDIF
        ENDDO
        IF(NFIR0>0) THEN
          CALL GXMAKE0(L0FIR,NFIR0,'FIRE')
          CALL GXMAKE0(L0QSOR,NFIR0,'QSOR')
          QTABLE='heat_sources.csv'
          LSTART=FSTEP==1
          IF(MYID==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)=='FIR') THEN
              LL=LENGZZ(PATNM)
              IF(PATNM(LL:LL)=='A') THEN
                NFIR=NFIR+1
                F(L0FIR+NFIR)=IP
                IF(NFIR==NFIR0) EXIT
              ENDIF
            ENDIF
          ENDDO
          IF(LSTART.AND.MYID==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>0) THEN
          CALL GXMAKE0(L0SMK,NSMOK0,'SMK')
          CALL GXMAKE0(L0SSOR,NSMOK0,'SSOR')
          STABLE='smoke_sources.csv'
          LSTART=FSTEP==1
          IF(MYID==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)=='FIR') THEN
              LL=LENGZZ(PATNM)
              IF(PATNM(LL:LL)/='A') THEN
                NSMOK=NSMOK+1
                F(L0SMK+NSMOK)=IP
                IF(NSMOK==NSMOK0) EXIT
              ENDIF
            ENDIF
          ENDDO
          IF(LSTART.AND.MYID==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>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)=='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
C ... Aerosol dispersion & deposition model - initialisation II
      IF(DFMODL) CALL GXDFM
      RETURN
C*****************************************************************
C
C--- GROUP 8. Terms (in differential equations) & devices
C
    8 GO TO (81,81,81,81,81,81,81,81,89,810,81,81,81,81,81,81)
     1,ISC
   81 RETURN
   89 CONTINUE
C   * ------------------- SECTION 9 ---- Diffusion coefficients
      IF(INDVAR==JAGE2) THEN
C.... Set diffusion coefficient to zero at entry and exit from room
        L0RM=L0F(JRM); NXMX=NX; NYMX=NY; NZMX=NZ
        IF(NDIREC==1) THEN ! N-S
          L0DF=L0F(LAN); IPLUS=1;    NYMX=NY-1
        ELSEIF(NDIREC==3) THEN ! E-W
          L0DF=L0F(LAE); IPLUS=NY;   NXMX=NX-1
        ELSEIF(NDIREC==5) THEN ! H-L
          L0DF=L0F(LD11); IPLUS=NFM; NZMX=NZ-1
        ENDIF
!        call writ3i('isweep ',isweep,' ndirec ',ndirec,' iz ',iz)
!        call prn('dif1',-L0df)
        DO IRM=1,NROOM ! loop over rooms
          IPTG=IRMPAT(IRM) ! get patch number for this room
          IF(LPAR) THEN
            IPT=GD_INDPAT(IPTG,1)   !  Local Index of Patch
          ELSE
            IPT=IPTG
          ENDIF
          IF(IPT<0) CYCLE ! not on this processor, skip
          CALL GETPAT(IPT,IR,TYP,IX1,IX2,IY1,IY2,IZ1,IZ2,IT1,IT2) ! get limits
          IZ1=MAX(1,IZ1-1); IZ2=MIN(NZMX,IZ2+1)
          IF(IZIZ2) CYCLE ! slab is not in range of room
          IX1=MAX(1,IX1-1); IX2=MIN(NXMX,IX2+1)
          IY1=MAX(1,IY1-1); IY2=MIN(NYMX,IY2+1)
          DO IX=IX1,IX2
            DO IY=IY1,IY2
              I=(IX-1)*NY+IY
              IF(F(L0RM+I)/=IRM.AND.F(L0RM+I+IPLUS)==IRM) THEN ! at low face
                F(L0DF+I)=0.0
              ELSEIF(F(L0RM+I)==IRM.AND.F(L0RM+I+IPLUS)/=IRM) THEN ! at high face
                F(L0DF+I)=0.0
              ENDIF
            ENDDO ! end IY loop
          ENDDO   ! end IX loop
        ENDDO     ! end room loop
!        call prn('dif2',-L0df)
      ENDIF
      RETURN
  810 CONTINUE
C   * ------------------- SECTION 10 --- Convection neighbours
      IF(INDVAR==JAGE2) THEN
C... Set convection neighbour to zero at inflow to room
        L0RM=L0F(JRM); L0NEI1=L0F(LD7); L0NEI2=L0F(LD8)
        NXMAX=NX; NYMAX=NY
        IF(NDIREC==1) THEN ! N-S
          L0VEL=L0F(V1); IPLUS=1; NYMAX=NY-1; IPLS=IPLUS
        ELSEIF(NDIREC==3) THEN ! E-W
          L0VEL=L0F(U1); IPLUS=NY; NXMAX=NX-1; IPLS=IPLUS
        ELSEIF(NDIREC==5) THEN ! H
          L0VEL=L0F(W1); IPLUS=NFM; L0NEI2=L0F(LD7); IPLS=0
        ELSEIF(NDIREC==6) THEN ! L
          L0VEL=L0F(W1)-NFM; IPLUS=-NFM; L0NEI2=L0F(LD7); IPLS=0
        ENDIF
!        call writ3i('isweep ',isweep,' ndirec ',ndirec,' iz ',iz)
!        call prn('age2',jage2)
!        call prn('ne11',ld7)
!        if(ndirec<5) call prn('ne12',ld8)
        DO IRM=1,NROOM ! loop over rooms
          IPTG=IRMPAT(IRM) ! get patch number for this room
          IF(LPAR) THEN
            IPT=GD_INDPAT(IPTG,1)   !  Local Index of Patch
          ELSE
            IPT=IPTG
          ENDIF
          IF(IPT<0) CYCLE ! not on this processor, skip
          CALL GETPAT(IPT,IR,TYP,IX1,IX2,IY1,IY2,IZ1,IZ2,IT1,IT2) ! get limits
          IZ1=MAX(1,IZ1-1); IZ2=MIN(NZ,IZ2+1)
          IF(IZIZ2) CYCLE ! slab is not in range of room
          IX1=MAX(1,IX1-1); IX2=MIN(NXMAX,IX2+1)
          IY1=MAX(1,IY1-1); IY2=MIN(NYMAX,IY2+1)
          DO IX=IX1,IX2
            DO IY=IY1,IY2
              I=(IX-1)*NY+IY
              IF(F(L0RM+I)/=IRM.AND.F(L0RM+I+IPLUS)==IRM) THEN ! at low face
                IF(F(L0VEL+I)>=0.0) THEN ! Inflow - set to zero
                  F(L0NEI1+I+IPLS)=0.0
                ELSE                     ! Outflow - do nothing
                ENDIF
              ELSEIF(F(L0RM+I)==IRM.AND.F(L0RM+I+IPLUS)/=IRM) THEN ! at high face
                IF(F(L0VEL+I)<=0.0) THEN ! Inflow - set to zero
                  F(L0NEI2+I)=0.0
                ELSE                     ! Outflow - do nothing
                ENDIF
              ENDIF
            ENDDO ! end IY loop
          ENDDO   ! end IX loop
        ENDDO     ! end room loop
!        call prn('ne21',ld7)
!        if(ndirec<5) call prn('ne22',ld8)
      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==12) THEN
C------------------- SECTION 12 ------------------- value = GRND
        IF(NPATCH(1:2)=='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==1.AND.ISINZZ(NPATCH(4:4))) ILEN=2
          IF(ILEN==2.AND.ISINZZ(NPATCH(5:5))) ILEN=3
          IF(ILEN==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==13) THEN
        IF(INDVAR==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
C ... Aerosol dispersion & deposition model - Slip & deposition sources
      IF(DFMODL) CALL GXDFM
      RETURN
C***************************************************************
C
C--- GROUP 19. Special calls to GROUND from EARTH
C
   19 IF(ISC==1) THEN
C   * ------------------- SECTION 1 ---- Start of time step.
C
C.. Initialise link temperature on 1st time step
        IF(.NOT.STEADY.AND.NSPRAY>0) THEN
          IF(ISTEP==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)
                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)>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)>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==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)==0.AND.CSG1(1:2)=='SW') .OR.
     1           MOD(ISWEEP+1,NPRINT)==0) .OR.
     1           ENUFSW.OR.ISWEEP==LSWEEP-1.OR.ISWEEP==LSWEEP
        ELSE
          LOG1= (MOD(ISTEP,IFREQ)==0 .OR.
     1           MOD(ISTEP,NTPRIN)==0) .AND.
     1          (ENUFSW.OR.ISWEEP==LSWEEP-1.OR.ISWEEP==LSWEEP)
        ENDIF
        IF(((LTRES.AND.LOG1).OR.(LPMV.AND.INDVAR==JTEM1).OR.
     1     NSPRAY>0.OR.LPPDR.OR.LTAPP.OR.LTAPR.OR.LTINS.OR.NROOM>0
     1                                   .OR.LUTCI).AND..NOT.LABSV) 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
          L0AVL=L0ABSV+(IZ-1)*NXNY
          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, skipping blocked faces
              IF(NX>1) THEN
                AU1=0.0
                IF(IX==1) THEN
                  IF(.NOT.SLD(I).AND..NOT.SLD(I+NY)) AU1=F(L0U1+I)**2
                ELSEIF(IX==NX) THEN
                  IF(.NOT.SLD(I).AND..NOT.SLD(I-NY)) AU1=F(L0U1+I-NY)**2
                ELSE
                  IDIV=0
                  IF(.NOT.SLD(I)) THEN
                    AU1=F(L0U1+I); IDIV=IDIV+1
                  ENDIF
                  IF(.NOT.SLD(I-NY)) THEN
                    AU1=AU1+F(L0U1+I-NY); IDIV=IDIV+1
                  ENDIF
                  AU1=(AU1/(IDIV+TINY) )**2
                ENDIF
              ENDIF
              IF(NY>1) THEN
                BV1=0.0
                IF(IY==1) THEN
                  IF(.NOT.SLD(I).AND..NOT.SLD(I+1)) BV1=F(L0V1+I)**2
                ELSEIF(IY==NY) THEN
                  IF(.NOT.SLD(I).AND..NOT.SLD(I-1)) BV1=F(L0V1+I-1)**2
                ELSE
                  IDIV=0
                  IF(.NOT.SLD(I)) THEN
                    BV1=F(L0V1+I); IDIV=IDIV+1
                  ENDIF
                  IF(.NOT.SLD(I-1)) THEN
                    BV1=BV1+F(L0V1+I-1); IDIV=IDIV+1
                  ENDIF
                  BV1=(BV1/(IDIV+TINY) )**2
                ENDIF
              ENDIF
              IF(NZ>1) THEN
                CW1=0.0
                IF(IZ==1) THEN
                  IF(.NOT.SLD(I).AND..NOT.SLD(I+NXNY)) CW1=F(L0W1+I)**2
                ELSEIF(IZ==NZ) THEN
                  IF(.NOT.SLD(I).AND..NOT.SLD(I-NXNY)) CW1=F(L0W1L+I)**2
                ELSE
                  IDIV=0
                  IF(.NOT.SLD(I)) THEN
                    CW1=F(L0W1+I); IDIV=IDIV+1
                  ENDIF
                  IF(.NOT.SLD(I-NXNY)) THEN
                    CW1=CW1+F(L0W1L+I); IDIV=IDIV+1
                  ENDIF
                  CW1=(CW1/(IDIV+TINY) )**2
                ENDIF
              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==1) THEN
            VOLSUM=0.
            IF(JAGE>0) AGESUM=0.0
          ENDIF
          IF(JAGE>0) L0AGE=L0F(JAGE)
          DO IX=1,NX
            DO IY=1,NY
              I=(IX-1)*NY+IY
              IF(.NOT.SLD(I)) THEN
                IF(ISACTIVECELLS(IX,IY,IZ)) THEN
                  VOLSUM=VOLSUM+F(L0VL+I)
                  IF(JAGE>0) AGESUM=AGESUM+F(L0AGE+I)*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,NXNY
            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,NXNY
            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(NZ==1.OR.(INDVAR==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) ! in Pa
                    IF(LPSAT) F(L0PSAT+I)=PVSAT
C .. Water vapour Pressure (Pa)
                    PVAP = HRAT*(PRESS0+F(L0P1+I))/(GMWRAT+HRAT)
                    IF(LPVAP) F(L0PVAP+I)=PVAP ! in Pa
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(NZ==1.OR.INDVAR==JTEM1) THEN
          IF(LPMV) THEN
            IF(LABSV) THEN
              L0AVL=L0F(JABSV)
            ELSE
              L0AVL=L0ABSV+(IZ-1)*NXNY
            ENDIF
            L0TEM=L0F(JTEM1)
            IF(LTRAD) L0T3=L0F(JT3)
            IF(JTCL>0) THEN
              L0TCL=L0F(JTCL)
            ELSE
              L0TCL=L0TCLC+(IZ-1)*NXNY
            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
              VMETAB=0.0
            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==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(HCEPS.AND.N<=150) GO TO 19611
                  IF(ABS(XN-XF)>EPS.AND.N<=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>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<=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
          IF(LTAPP.OR.LTAPR) THEN
            L0TEM=L0F(JTEM1)
            IF(LPVAP) THEN
              L0PVAP=L0F(JPVAP)
            ELSE
              IF(JRELH>0) THEN
                L0RELH=L0F(JRELH)
              ELSE
                RELH=RH
              ENDIF
            ENDIF
            IF(LTAPP) THEN
              L0TAPP=L0F(JTAPP)
            ENDIF
            IF(LTAPR) THEN
              L0TAPR=L0F(JTAPR)
              IF(ALLCELLS.AND.JLIT>0) THEN
                L0LIT=L0F(JLIT)
              ELSE
                L0LIT=0; LIT=1
              ENDIF
            ENDIF
            IF(LABSV) THEN
              L0AVL=L0F(JABSV)
            ELSE
              L0AVL=L0ABSV+(IZ-1)*NXNY
            ENDIF
            DO I=1,NXNY
              IF(SLD(I)) CYCLE
              IF(LPVAP) THEN
                VAPR=0.01*F(L0PVAP+I) ! in hPa
              ELSE
                IF(JRELH>0) THEN
                  VAPR=0.01*F(L0RELH+I)*6.105*EXP(17.27*F(L0TEM+I)/
     1                                              (237.7+F(L0TEM+I)))
                ELSE
                  VAPR=0.01*RELH*6.105*EXP(17.27*F(L0TEM+I)/(237.7+
     1                                                      F(L0TEM+I)))
                ENDIF
              ENDIF
              IF(LTAPP) THEN
C...  AT = Ta + 0.33×e - 0.70×ws - 4.0
                F(L0TAPP+I)=F(L0TEM+I)+0.33*VAPR-0.7*F(L0AVL+I)-4.0
              ENDIF
              IF(LTAPR) THEN
C...   ATr = Ta + 0.348×e - 0.70×ws + 0.70×Q/(ws + 10) - 4.25
                IF(L0LIT>0) THEN
                  LIT=ITWO(1,0,F(L0LIT+I)>0.0)
                ENDIF
                IF(LIT>0) THEN
                  F(L0TAPR+I)=F(L0TEM+I)+0.348*VAPR-0.7*F(L0AVL+I)+
     1                                    0.7*QRAD/(F(L0AVL+I)+10.)-4.25
                ELSE
                  F(L0TAPR+I)=F(L0TEM+I)+0.33*VAPR-0.7*F(L0AVL+I)-4.0
                ENDIF
              ENDIF
            ENDDO
          ENDIF
C.... Universal Thermal Climate Index and Physiological Equivalent Temperature
          IF(LUTCI.OR.LPET) THEN
            IF(LABSV) THEN
              L0AVL=L0F(JABSV)
            ELSE
              L0AVL=L0ABSV+(IZ-1)*NXNY
            ENDIF
            L0TEM=L0F(JTEM1)
            IF(LTRAD) THEN
              L0T3=L0F(JT3)
            ELSE
              TRA=TRAD
            ENDIF
            IF(LPVAP) THEN
              L0PVAP=L0F(JPVAP)
            ELSE
              IF(QEQ(RH,-999.0).AND.LRELH) THEN
                L0RELH=L0F(JRELH)
              ELSE
                L0RELH=0; RELH=RH
              ENDIF
            ENDIF
            IF(LUTCI) L0UTCI=L0F(JUTCI)
            IF(LPET) THEN
              L0PET=L0F(JPET); L0P1=L0F(P1)
            ENDIF
            DO I=1,NXNY
              IF(SLD(I)) CYCLE
              VEL=F(L0AVL+I)
              TAIR=F(L0TEM+I)+TEMP0-273.0   ! in deg C
              IF(LPVAP) THEN
                VAPR=0.01*F(L0PVAP+I) ! in hPa, PVAP in Pa
              ELSE
                IF(L0RELH>0) RELH=F(L0RELH+I)
                VAPR=0.01*RELH*PVH2O(TAIR)/100. ! in hPa, PVH2O in Pa
              ENDIF
              IF(LTRAD) TRA=F(L0T3+I)+TEMP0-273.
              IF(LUTCI) F(L0UTCI+I)=UTCI_APPROX(TAIR,VAPR,TRA,VEL)
              IF(LPET) THEN
                P=0.01*(F(L0P1+I)+PRESS0)+TINY ! in hPa
                CALL CALCULATE_PET_STATIC(TAIR,VAPR,VEL,TRA,P,TX)
                F(L0PET+I)=TX
              ENDIF
            ENDDO
          ENDIF
        ENDIF
C...     Thermal Sensation Index (BEAM)
        IF(LTSIB) THEN
          IF(LABSV) THEN
            L0AVL=L0F(JABSV)
          ELSE
            L0AVL=L0ABSV+(IZ-1)*NXNY
          ENDIF
          L0TSIB=L0F(JTSIB)
          if(lbname('TSIF')>0) l0tsif=l0f(lbname('TSIF'))
          T=31.0; RH=72.0; ST=34.0; DHI=254.0; GHI=525.0
          IF(ALLCELLS.AND.JLIT>0) THEN
            L0LIT=L0F(JLIT)
          ELSE
            L0LIT=0; LIT=0
          ENDIF
          DO I=1,NXNY
            IF(SLD(I)) CYCLE
            VEL=F(L0AVL+I)
            IF(L0LIT>0) LIT=ITWO(1,0,F(L0LIT+I)>0.0)
            SR=DHI+(GHI-DHI)*LIT
            GTSIB=1.7+0.118*T+0.0019*SR-0.322*VEL-0.0073*RH
     1                                                 +0.0054*ST
            if(lbname('TSIF')>0) f(l0tsif+i)=gtsib
            GTSIB=MAX(1.0,MIN(GTSIB,7.0))
            F(L0TSIB+I)=INT(GTSIB)
          ENDDO
        ENDIF
C... Draught rating (as per ISO 7730)
        IF(LPPDR.OR.LTINS) THEN
          IF(LPPDR) THEN
            L0PPDR=L0F(JPPDR); L0TEM=L0F(JTEM1)
            VMIN=VARMIN(JPPDR); VMAX=VARMAX(JPPDR)
          ENDIF
          IF(LTINS) L0TINS=L0F(JTINS)
          IF(LABSV) THEN
            L0AVL=L0F(JABSV)
          ELSE
            L0AVL=L0ABSV+(IZ-1)*NXNY
          ENDIF
          IF(ITURB==1) THEN
            L0KE=L0F(KE)
          ELSEIF(ITURB==2) THEN
            L0ENUT=L0F(VIST); L0EL1=L0F(LEN1)
          ENDIF
          DO I=1,NXNY
            IF(.NOT.SLD(I)) THEN
              VABS=AMAX1(0.05,F(L0AVL+I))
              IF(ITURB==0) THEN   ! laminar
                TINS=0.0
              ELSE
                IF(ITURB==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
              TINSP=AMAX1(10.,AMIN1(TINS,60.))
              IF(LPPDR) THEN
                PPDR=(34.0-F(L0TEM+I))*((VABS-0.05)**0.62)*
     1                                        (0.37*VABS*TINSP+3.14)
                F(L0PPDR+I)=AMAX1(VMIN,AMIN1(VMAX,PPDR))
              ENDIF
            ELSE
              IF(LPPDR) F(L0PPDR+I)=0.0; TINS=0.0
            ENDIF
            IF(LTINS) F(L0TINS+I)=MAX(VARMIN(JTINS),
     1                                       MIN(TINS,VARMAX(JTINS)))
          ENDDO
        ENDIF
      ELSEIF(ISC==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)==0.AND.CSG1=='SW') .OR.
          LOG1=((MOD(ISWEEP,IFREQ)==0.AND.CSG1(1:2)=='SW') .OR.
     1           MOD(ISWEEP,NPRINT)==0) .OR.
     1           ENUFSW.OR.ISWEEP==LSWEEP
          RPRINT=((MOD(ISWEEP,IFREQ)==0.AND.CSG1(1:2)=='SW') .OR.
     1           MOD(ISWEEP,NPRINT)==0) .OR.
     1           (MOD(ISWEEP,NPLT)==0) .OR.
     1           ENUFSW.OR.ISWEEP==LSWEEP
        ELSE
          LOG1= (MOD(ISTEP,IFREQ)==0 .OR.
     1           MOD(ISTEP,NTPRIN)==0) .AND.
     1          (ENUFSW.OR.ISWEEP==LSWEEP)
          RPRINT=(MOD(ISTEP,IFREQ)==0 .OR.
     1            MOD(ISTEP,NTPRIN)==0) .AND.
     1           (ENUFSW.OR.ISWEEP==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)*NXNY
            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; TBAR=0.0; CBAR=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(GCO,-999.)) CYCLE
            IF(IR<0) THEN
              SORCE=0.0
            ELSE
              CALL GETSO(IR,R1,SORCE)
            ENDIF
            IF(LPAR) CALL GLSUM(SORCE)
            IF(SORCE>0.0) THEN ! mass inflow
              FMASIN=FMASIN+SORCE  ! sum mass
              IF(JTEM1>0) THEN
                IF(IR>0) THEN
                  CALL GETSO(IR,JTEM1,SOR)
                ELSE
                  SOR=0.0
                ENDIF
                IF(LPAR) CALL GLSUM(SOR)
                TBAR=TBAR+SOR
              ENDIF
              IF(SOLVE(C1)) THEN
                IF(IR>0) THEN
                  CALL GETSO(IR,C1,SOR)
                ELSE
                  SOR=0.0
                ENDIF
                IF(LPAR) CALL GLSUM(SOR)
                CBAR=CBAR+SOR
              ENDIF
            ENDIF
          ENDDO
C
          RHOIN=1.189 ! default inlet density
          IF(CP1>0.AND.JTEM1>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
              IF(NEZ(RHO1A)) THEN
                CBAR=CBAR/FMASIN
                RHOIN=PRESS0/(((1.-CBAR)*RHO1A+CBAR*RHO1B)*TBAR+TINY)
              ELSE
                RHOIN=PRESS0*RHO1B/(TBAR+TINY) ! get average inlet density
              ENDIF
            ENDIF
          ENDIF
          IF(LPAR) CALL GLSUM2(VOLSUM,AGESUM)
          VFLRT=AMAX1(TINY,FMASIN/RHOIN)
          RESTIME=VOLSUM/(VFLRT+TINY)
          IF(LOG1.AND.MYID==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 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
            IF(JAGE>0) THEN
              RMAA=AGESUM/VOLSUM
              ACE=RESTIME/RMAA
              CALL WRIT40('Mean Age of Air in Domain (s)')
              CALL WRIT1R('MAA     ',RMAA)
              CALL WRITBL
              CALL WRIT40('Air Change Effectiveness for domain'
     1                                       //' (RES.TIME/MAA)')
              CALL WRIT1R('ACE     ',ACE)
              CALL WRITBL
            ENDIF
            CALL WRITST
          ENDIF
C
C.... Compute room ACH values
          L0FACZ=L0FACE; IRM=0; IZ=1; I0=1
          jmkx=lbname('MKX'); jmky=lbname('MKY'); jmkz=lbname('MKZ')
          DO IRM=1,NROOM ! loop over rooms
            IPTG=IRMPAT(IRM) ! get patch number for this room
            IF(LPAR) THEN
              IPT=GD_INDPAT(IPTG,1)   !  Local Index of Patch
            ELSE
              IPT=IPTG
            ENDIF
            FLOX=0.0; FLOY=0.0; FLOZ=0.0; SOR=0.0; VROOM=0.0
            CROOM=0.0; FLOXC=0.0; FLOYC=0.0; FLOZC=0.0; CROOM0=0.0
            TOTAE=TINY; TOTAN=TINY; TOTAH=TINY
            TROOM=0.; TMAX=-9999.; TMIN=9999.; COUT=0.0; MOUT=0.0
            VABTOT=0.0; AHTOT=0.0; ZVABS=HVABS(IRM)
            IF(IPT<0) THEN ! not on this processor, skip
!              do nothing
            ELSE  ! process this patch
              CALL GETPAT(IPT,IR,TYP,IX1,IX2,IY1,IY2,IZ1,IZ2,IT1,IT2) ! get limits
              IX1=MAX(1,IX1-1); IX2=MIN(NX,IX2+1)
              IY1=MAX(1,IY1-1); IY2=MIN(NY,IY2+1)
              IZ1=MAX(1,IZ1-1); IZ2=MIN(NZ,IZ2+1)
              ZFLOOR=-999.; ZL=0.0
              DO IZZ=IZ1,IZ2
                IF(NX>1) L0U=L0F(ANYZ(U1,IZZ))
                IF(NY>1) L0V=L0F(ANYZ(V1,IZZ));
                IF(NZ>1) L0W=L0F(ANYZ(W1,IZZ))
                L0RM=L0F(ANYZ(JRM,IZZ)); L0D=L0F(ANYZ(DEN1,IZZ))
                IF(JTEM1>0) L0T=L0F(ANYZ(JTEM1,IZZ))
                IF(JAGE2>0) L0AGE2=L0F(ANYZ(JAGE2,IZZ))
                IF(JAGE>0)  L0AGE0=L0F(ANYZ(JAGE,IZZ))
                IF(LABSV) THEN
                  L0AVL=L0F(ANYZ(JABSV,IZZ))
                ELSE
                  L0AVL=L0ABSV+(IZZ-1)*NXNY
                ENDIF
                if(jmkx>0) l0mkx=l0f(anyz(jmkx,izz))
                if(jmky>0) l0mky=l0f(anyz(jmky,izz))
                if(jmkz>0) l0mkz=l0f(anyz(jmkz,izz))
                IZAD=(IZZ-1)*NXNY
                L0AE=I3DAEX+IZAD; L0AN=I3DANY+IZAD; L0AH=I3DAHZ+IZAD
                L0VOL=I3DVOL+IZAD
                ZH=F(KZZW+IZZ) ! upper face of slab
                DO IX=IX1,IX2
                  IXAD=(IX-1)*NY
                  DO IY=IY1,IY2
                    IF(.NOT.ISACTIVECELLS(IX,IY,IZZ)) CYCLE
                    I=IXAD+IY; IJK=I+IZAD
                    if(jmkx>0) f(l0mkx+i)=0.0
                    if(jmky>0) f(l0mky+i)=0.0
                    if(jmkz>0) f(l0mkz+i)=0.0
                    IF(.NOT.SLD(IJK)) THEN    ! sum open room volume
                      IF(F(L0RM+I)==IRM) THEN ! cell is in room IRM
                        IF(ZFLOOR<0.0) ZFLOOR=ZL ! save height of lower face of room
                        VROOM=VROOM+F(L0VOL+I) ! sum room volume
                        IF(JTEM1>0) THEN       ! get volume-average room temperature
                          TROOM=TROOM+F(L0T+I)*F(L0VOL+I)
                          TMAX=MAX(TMAX,F(L0T+I)) ! get min/max room temperature
                          TMIN=MIN(TMIN,F(L0T+I))
                        ENDIF
                        IF(JAGE2>0) CROOM =CROOM +F(L0AGE2+I)*F(L0VOL+I) ! volume-average age
                        IF(JAGE>0)  CROOM0=CROOM0+F(L0AGE0+I)*F(L0VOL+I)
                        IF(ZL-ZFLOOR<=ZVABS.AND.ZH-ZFLOOR>ZVABS) THEN ! slab is ~1.2m above floor
                          VABTOT=VABTOT+F(L0AVL+I)*F(L0AH+I) ! area-average VABS
                          AHTOT=AHTOT+F(L0AH+I)
                        ENDIF
                      ENDIF ! end cell-in-room block
                    ENDIF
                    IF(NX>1) THEN ! X-directed fluxes
                      IF(.NOT.LSLDIR(IJK,3)) THEN ! skip blocked faces
                        IAD=ITWO(NY,0,IX0) THEN
                            IF(F(L0U+I)<0.0) THEN ! outflow
                              FLOUT=ABS(F(L0AE+I)*F(L0D+I+IAD)*F(L0U+I))
                              FLOXC=FLOXC+FLOUT*F(L0AGE2+I+IAD)
                              TOTAE=TOTAE+FLOUT
                            ENDIF
                          ENDIF
                          if(jmkx>0) f(l0mkx+i)=IRM
                        ELSEIF(F(L0RM+I)==IRM.AND.F(L0RM+I+IAD)/=IRM)
     1                                                              THEN ! high face
                          FLOX=FLOX+F(L0AE+I)*MAX(-F(L0U+I),0.0) ! inflow
                          IF(JAGE2>0) THEN
                            IF(F(L0U+I)>0.0) THEN                ! outflow
                              FLOUT=F(L0AE+I)*F(L0D+I)*F(L0U+I)
                              FLOXC=FLOXC+FLOUT*F(L0AGE2+I)
                              TOTAE=TOTAE+FLOUT
                            ENDIF
                          ENDIF
                          if(jmkx>0) f(l0mkx+i)=-IRM
                        ENDIF ! end faces
                      ENDIF ! end blocked face
                    ENDIF ! end NX>1
                    IF(NY>1) THEN ! Y-directed fluxes
                      IF(.NOT.LSLDIR(IJK,1)) THEN ! skip blocked faces
                        IAD=ITWO(1,0,IY0) THEN
                            IF(F(L0V+I)<0.0) THEN ! outflow
                              FLOUT=ABS(F(L0AN+I)*F(L0D+I+IAD)*F(L0V+I))
                              FLOYC=FLOYC+FLOUT*F(L0AGE2+I+IAD)
                              TOTAN=TOTAN+FLOUT
                            ENDIF
                          ENDIF
                          if(jmky>0) f(l0mky+i)=IRM
                        ELSEIF(F(L0RM+I)==IRM.AND.F(L0RM+I+IAD)/=IRM)
     1                                                              THEN ! high face
                          FLOY=FLOY+F(L0AN+I)*MAX(-F(L0V+I),0.0) ! inflow
                          IF(JAGE2>0) THEN
                            IF(F(L0V+I)>0.0) THEN                ! outflow
                              FLOUT=F(L0AN+I)*F(L0D+I)*F(L0V+I)
                              FLOYC=FLOYC+FLOUT*F(L0AGE2+I)
                              TOTAN=TOTAN+FLOUT
                            ENDIF
                          ENDIF
                          if(jmky>0) f(l0mky+i)=-IRM
                        ENDIF
                      ENDIF
                    ENDIF ! end NY>1
                    IF(NZ>1) THEN ! Z-directed fluxes
                      IF(.NOT.LSLDIR(IJK,5)) THEN ! skip blocked faces
                        IAD=ITWO(NFM,0,IZZ0) THEN
                            IF(F(L0W+I)<0.0) THEN               ! outflow
                              FLOUT=ABS(F(L0AH+I)+F(L0D+I+IAD)*F(L0W+I))
                              FLOZC=FLOZC+FLOUT*F(L0AGE2+I+IAD)
                              TOTAH=TOTAH+FLOUT
                            ENDIF
                          ENDIF
                          if(jmkz>0) f(l0mkz+i)=IRM
                        ELSEIF(F(L0RM+I)==IRM.AND.F(L0RM+I+IAD)/=IRM)
     1                                                              THEN ! high face
                          FLOZ=FLOZ+F(L0AH+I)*MAX(-F(L0W+I),0.0) ! inflow
                          IF(JAGE2>0) THEN
                            IF(F(L0W+I)>0.0) THEN
                              FLOUT=F(L0AH+I)*F(L0D+I)*F(L0W+I)  ! outflow
                              FLOZC=FLOZC+FLOUT*F(L0AGE2+I)
                              TOTAH=TOTAH+FLOUT
                            ENDIF
                          ENDIF
                          if(jmkz>0) f(l0mkz+i)=-IRM
                        ENDIF
                      ENDIF
                    ENDIF ! end NZ>1
                  ENDDO ! end IY loop
                ENDDO   ! end IX loop
                ZL=ZH   ! move upper-face to lower-face for next slab
              ENDDO     ! end IZ loop
              NSOR=IRMSOR(I0) ! get number of internal sources
              IF(NSOR>0) THEN ! there are sources
                DO IS=1,NSOR  ! loop over sources
                  IPS=IRMSOR(IS+I0) ! get patch number for this source
                  CALL GETSO(IPS,P1,SORP) ! get mass source
                  IF(SORP>0) THEN ! source is +ve, so inflow.
                    IF(RHO1>0.0) THEN ! constant density
                      RHOIN=RHO1
                    ELSEIF(QEQ(RHO1,GRND5).AND.JTEM1>0.AND.CP1>0.) THEN ! Ideal Gas
                      CALL GETSO(IPS,JTEM1,GQSOR) ! T1 source=mCpT
                      T1AVE=GQSOR/SORP/CP1        ! ave T=T1sor/mCp
                      IF(EQZ(RHO1A)) THEN         ! ideal gas
                        SPEVOL=1./(RHO1B+TINY)
                      ELSE       ! ideal gas, mixture Gas constant
                        CALL GETSO(IPS,C1,C1SOR) ! inlet C1
                        C1AVE=C1SOR/SORP  ! average C1
                        SPEVOL=RHO1A+(RHO1B-RHO1A)*C1AVE
                      ENDIF
                      RHOIN=PRESS0/(SPEVOL*T1AVE+TINY)
                    ELSE ! unknown density formula - assume air
                      RHOIN=1.189
                    ENDIF
                    SOR=SOR+SORP/RHOIN ! sum volume inflow
                  ELSEIF(SORP<0.0) THEN ! outflow
                    IF(JAGE2>0) THEN
                      CALL GETSO(IPS,JAGE2,GCOUT) ! AGE source is Mout.Cout
                      COUT=COUT+ABS(GCOUT)
                      MOUT=MOUT+ABS(SORP) ! sum mass outflow. Ave C is Mout.Cout/Mout
                    ENDIF
                  ENDIF
                ENDDO
              ENDIF
              I0=I0+NSOR+1 ! increment offset
            ENDIF ! end of patch block
            IF(LPAR) THEN ! in parallel sum over all processors
              VALS(1)=FLOX; VALS(2)=FLOY; VALS(3)=FLOZ; VALS(4)=SOR
              VALS(5)=VROOM; VALS(6)=TROOM; VALS(7)=CROOM;VALS(8)=VABTOT
              VALS(9)=AHTOT; NSUM=9
              IF(JAGE2>0.OR.JAGE>0) THEN ! AEE parameters
                VALS(NSUM+1)=FLOXC; VALS(NSUM+2)=FLOYC
                VALS(NSUM+3)=FLOZC; VALS(NSUM+4)=TOTAE
                VALS(NSUM+5)=TOTAN; VALS(NSUM+6)=TOTAH
                VALS(NSUM+7)=COUT;  VALS(NSUM+8)=MOUT; NSUM=NSUM+8
                IF(JAGE>0) THEN
                  VALS(NSUM+1)=CROOM0; NUM=NSUM+1
                ENDIF
              ENDIF
              CALL DGSUMN(VALS,NSUM)
              FLOX=VALS(1); FLOY=VALS(2); FLOZ=VALS(3); SOR=VALS(4)
              VROOM=VALS(5); TROOM=VALS(6); CROOM=VALS(7);VABTOT=VALS(8)
              AHTOT=VALS(9)
              IF(NSUM>9) THEN
                IF(JAGE2>0) THEN
                  FLOXC=VALS(10); FLOYC=VALS(11); FLOZC=VALS(12)
                  TOTAE=VALS(13); TOTAN=VALS(14); TOTAH=VALS(15)
                  COUT=VALS(16); MOUT=VALS(17)
                ENDIF
                IF(JAGE>0) CROOM0=VALS(NSUM)
              ENDIF
              IF(JTEM1>0) THEN ! min/max TEM1 over all processors
                CALL PAR_MAXR(TMAX); CALL PAR_MINR(TMIN)
              ENDIF
            ENDIF ! end parallel block
            VABTOT=VABTOT/(AHTOT+TINY) ! Area-averaged absolute velocity
            IF(JAGE>0) THEN ! get Air Change Effectiveness based on AGE
              CROOM0=CROOM0/VROOM ! average AGE in room aka MAA
              ACE=VROOM/(SOR+FLOX+FLOY+FLOZ)/CROOM0 ! AGE-based Air Change Effectiveness
            ENDIF
            IF(JAGE2>0) THEN ! get and store Air Exchange Effectiveness
              CROOM=CROOM/VROOM
              AEE=(COUT+FLOXC+FLOYC+FLOZC)/(MOUT+TOTAE+TOTAN+TOTAH)/
     1                                                             CROOM
              ACE2=VROOM/(SOR+FLOX+FLOY+FLOZ)/CROOM ! AGE2-based Air Change Efficiency
              IF((JAEE>0..OR.JACE>0).AND.IPT>0) THEN ! AEE or ACE are STOREd so fill it with AEE for
                DO IZZ=IZ1,IZ2
                  L0RM=L0F(ANYZ(JRM,IZZ))
                  IF(JAEE>0) L0AEE=L0F(ANYZ(JAEE,IZZ))
                  IF(JACE>0) L0ACE=L0F(ANYZ(JACE,IZZ))
                  DO IX=IX1,IX2
                    IXAD=(IX-1)*NY
                    DO IY=IY1,IY2
                      IF(.NOT.ISACTIVECELLS(IX,IY,IZZ)) CYCLE
                      I=IXAD+IY; IF(F(L0RM+I)/=IRM) CYCLE ! cell is not in room IRM
                      IF(JAEE>0) F(L0AEE+I)=AEE
                      IF(JACE>0) F(L0ACE+I)=ACE
                    ENDDO
                  ENDDO
                ENDDO
              ENDIF ! end of AEE block
            ENDIF ! end of JAGE2>0 block
            IF(LOG1.AND.MYID==0.AND..NOT.NULLPR) THEN ! print on master
              FLUXT=FLOX+FLOY+FLOZ+SOR
              IF(LPAR) THEN
                WRITE(14,'(A)') ' For object '//GD_OBJNAM(IPTG)
              ELSE
                WRITE(14,'(A)') ' For object '//OBJNAM(IPT)
              ENDIF
              WRITE(14,'(A)') ' ---------------------'
              CALL WRIT40('Overall residence time calculated as    ')
              CALL WRIT40('free volume/volum.flow-rate (in seconds)')
              CALL WRIT1R('RES.TIME',VROOM/FLUXT)
              CALL WRIT40('Ventilation rate in air changes per hour')
              CALL WRIT1R('ACH     ',3600.*FLUXT/VROOM)
              IF(IRM<10) THEN
                WRITE(BUFF,'(''ACH_'',I1)') IRM
              ELSEIF(IRM<100) THEN
                WRITE(BUFF,'(''ACH_'',I2)') IRM
              ELSEIF(IRM<1000) THEN
                WRITE(BUFF,'(''ACH_'',I3)') IRM
              ENDIF
              LL=LENGZZ(BUFF)
              CALL SET_INF(BUFF(1:LL),FLUXT/VROOM,IERR)
              CALL WRIT40(
     1                     'The total free volume in the room is (m^3)')
              CALL WRIT1R('VOLUME  ',VROOM)
              CALL WRIT40
     1             ('The total volumetric flow rate is (m^3/s, m^3/hr)')
              CALL WRIT2R('VFLRT s ',FLUXT,';VFLRT h',FLUXT*3600.)
              IF(IRM<10) THEN
                WRITE(BUFF,'(''VFLRT_'',I1)') IRM
              ELSEIF(IRM<100) THEN
                WRITE(BUFF,'(''VFLRT_'',I2)') IRM
              ELSEIF(IRM<1000) THEN
                WRITE(BUFF,'(''VFLRT_'',I3)') IRM
              ENDIF
              LL=LENGZZ(BUFF)
              CALL SET_INF(BUFF(1:LENGZZ(BUFF)),FLUXT,IERR)
              IF(NEZ(SOR)) THEN
                CALL WRIT40
     1       ('The total internal volumetric source is (m^3/s, m^3/hr)')
                CALL WRIT2R('VFLRT s ',SOR,';VFLRT h',SOR*3600.)
              ENDIF
              IF(JTEM1>0) THEN
                CALL WRIT40('The volume-weighted average temperature')
                CALL WRIT1R('Ave Temp',TROOM/VROOM)
                CALL WRIT40('The minimum and maximum temperatures')
                CALL WRIT2R('Min Temp',TMIN,';Max temp',TMAX)
              ENDIF
              IF(JAGE2>0) THEN ! get and print Air Exchange Effectiveness
                CALL WRIT40('The Air Exchange Effectiveness')
                CALL WRIT40(' i.e. (ave outflow local age)/'
     1                                  //'(ave local age in room)')
                CALL WRIT1R('AEE',AEE)
                CALL WRITBL
                CALL WRIT40('The local Air Change Effectiveness')
                CALL WRIT40
     1                  (' i.e. RES.TIME/(average local age for room)')
                CALL WRIT40('Local Mean Age of Air for room')
                CALL WRIT1R('MAA local',CROOM)
                CALL WRIT1R('ACE local',ACE2)
              ENDIF ! end of AEE block
              IF(JAGE>0) THEN ! get and print Air Change Effectiveness
                CALL WRITBL
                CALL WRIT40('The global Air Change Effectiveness')
                CALL WRIT40
     1                  (' i.e. RES.TIME/(average domain age for room)')
                CALL WRIT40('Global Mean Age of Air for room')
                CALL WRIT1R('MAA global',CROOM0)
                CALL WRIT1R('ACE global',ACE)
              ENDIF
              CH13=CHGR_13_R(ZVABS); LL=LENGZZ(CH13)
              CALL WRITBL
              CALL WRIT40
     1     ('The area-averaged absolute velocity at '//CH13(1:LL)//'m')
              IF(EQZ(AHTOT)) THEN
                CALL WRIT40('WARNING - this height is above the room!')
              ELSE
                CALL WRIT1R('VAB_AVE',VABTOT)
                CALL WRIT40('The total free area of the room is (m^2)')
                CALL WRIT1R('AREA   ',AHTOT)
              ENDIF
              CALL WRITBL
            ENDIF ! end of print block
          ENDDO  ! end loop over rooms DO IRM
        ENDIF ! end IF(LOG1 .OR. RPRINT
C
C... Print normalised residuals to csv residual table
        IF(RPRINT) THEN
          ISTO=0
          DO IVAR=1,NPHI
            IF(SOLVE(IVAR).AND.NAME(IVAR)/='LTLS') THEN
              ISTO=ISTO+1; TOTR=TOTRES(IVAR)
              RESID(ISTO)=100.*TOTR
            ENDIF
          ENDDO
C
          IF(MYID==0.AND..NOT.NULLPR) THEN
!            CALL WRITST
            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)
          ENDIF
        ENDIF
C
C... Sprinkler-head activation
C
        IF(.NOT.STEADY.AND.NSPRAY>0) 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)*NXNY)
            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))
            IF(LTLINK) THEN
              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
            ENDIF
          ENDDO
        ENDIF
      ELSEIF(ISC==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)==0 .OR.
     1           MOD(ISTEP,NTPRIN)==0) .AND.
     1          (ENUFSW.OR.ISWEEP==LSWEEP-1.OR.ISWEEP==LSWEEP)
        ENDIF
        IF(LOG1.AND.MYID==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>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)>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)
          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==FSTEP) THEN
            LSTART=.TRUE.
            IS1=1;IS2=MIN(NSPRAY,NTB);ITAB=0
            DO WHILE (IS1<=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<=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<=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<0) THEN
                FSOR=0.0
              ELSE
                CALL GETSO(IP,JTEM1,FSOR)
              ENDIF
              CALL GLSUM(FSOR)
              F(L0QSOR+I)=FSOR
            ENDDO
            IF(MYID==0.AND.NFIR>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<0) THEN
                SSOR=0.0
              ELSE
                CALL GETSO(IP,JSMOK,SSOR)
              ENDIF
              CALL GLSUM(SSOR)
              F(L0SSOR+I)=SSOR
            ENDDO
            IF(MYID==0.AND.NSMOK>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>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>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
C... Clean up storage for room ach
        IF(ISTEP==LSTEP) THEN
          IF(ALLOCATED(IRMPAT)) DEALLOCATE(IRMPAT,IRMSOR,STAT=IERR)
          IF(ALLOCATED(HVABS)) DEALLOCATE(HVABS,STAT=IERR)
        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==1) THEN
        IF(LSTART) THEN ! create new file at start of run
          OPEN(LUT,FILE=TABLE,STATUS='NEW',IOSTAT=IOS)
          IF(IOS/=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/=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/=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/=0) GO TO 700
            OPEN(LUT,FILE=TABLE,STATUS='NEW',IOSTAT=IOS)
            LSTART=.TRUE.
          ENDIF
        ENDIF
      ELSE
        LL=LENGZZ(BUFF)
        IF(BUFF(LL:LL)==',') LL=LL-1 ! strip off any trailing ,
        OPEN(LUT,FILE=TABLE,STATUS='OLD',ACCESS='APPEND',IOSTAT=IOS)
        IF(IOS==0) THEN ! file exists
          WRITE(LUT,'(1X,A)',IOSTAT=IOS) BUFF(1:LL)
          IF(IOS/=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/=0) GO TO 700
          OPEN(LUT,FILE=TABLE,STATUS='NEW',IOSTAT=IOS)
          IF(IOS>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/=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==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 'farray'
      INCLUDE 'grdear'
      INCLUDE 'tracmn'
      INCLUDE '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/=0) GO TO 1000 ! problem opening it
      OPEN(LUOUT,FILE='spin_tmp',STATUS='UNKNOWN',IOSTAT=IOS) ! open scratch
      IF(IOS/=0) GO TO 1000 ! problem opening it
      CLOSE(LUOUT,STATUS='DELETE',IOSTAT=IOS) ! delete it
      IF(IOS/=0) GO TO 1000 ! problem opening it
      OPEN(LUOUT,FILE='spin_tmp',STATUS='NEW',IOSTAT=IOS) ! open new copy
      IF(IOS/=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:)==' '.OR.INDEX(LINE,'*')/=0) THEN ! found header or blank line
          LINLEN=LENGZZ(LINE)
          IF(INDEX(LINE,'ports')/=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 line
            IF(ERROR) GO TO 1001
          ELSEIF(INDEX(LINE,'size')/=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 line
            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:)==' '.OR.INDEX(LINE,'*')/=0) GOTO 435 ! skip comment or blank line
          LINLN1=LENGZZ(LINE)
          IF(NINT(F(L0ACTIV+ISP))==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 injection
            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 current time
              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/=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,IERR)         ! delete scratch
      IF(IOS/=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