c

C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C-------------------------------------------------------------------------
      SUBROUTINE GXBLIN
C-------------------------------------------------------------------------
C
C     PHOENICS V2011
C     Author:  M.R.Malin / J C Ludwig
C     Date:    24/09/04 - 30/09/15
C
C     Subroutine GXBLIN is called from Group 1 Section 1 and Group 13 Section
C     19 of GREX3 {using the logical condition IF(NPATCH(1:4=='BLIN')}.
C
C     The purpose of this subroutine is to set mass-flow boundary
C     conditions for the momentum and k-e equations so as to define a
C     boundary-layer profile at domain  boundaries. At present only
C     option is provided for an atmospheric boundary layer (wind profile)
C     using either a logarithmic or power-law velocity profile, as
C     follows:
C
C          Q = {Q*/K}*ln((z-d)/zo)  or   Q = Qref*(z/zref)^a
C
C          k = Q*^2/sqrt(cmucd) and  eps = Q*^3/(K*[z-d])
C
C     where the friction velocity Q* is given by:
C
C          Q*= Qref*K/log([zref-d]/zo)
C
C     and
C
C     K = 0.41, zo is the effective roughness height of the ground terrain,
C     z is the displacement height, Qref is the reference velocity at the
C     reference height zref, a is the power-law exponent, z is the vertical
C     coordinate, and Q is the velocity at the height z from the ground.
C
C     The displacement height d is a positive quantity, although the user
C     can set d=-z0 to create the inlet wind profiles proposed by Richards
C     & Hoxey (1993), which are often used in wind-engineering simulations.
C     ( see Richards,P.J. & Hoxey,R.P., "Appropriate boundary conditions for
C     computational wind engineering models using the k-e turbulence model."
C     J.Wind Engng & Industrial Aerodynamics, 46 & 47, 145-153, (1993) ).
C     	
C
C     The vertical coordinate is taken to start at the lower edge of the first
C     cut or fully-open cell in each column of an inlet plane.
C
C     The facility requires that the user defines a PATCH whose name starts
C     with the 4 characters BLIN at the inlet boundary, such as for example:
C
C   PATCH(BLIN1,LOW,1,NX,1,NY,1,1,1,1)
C   COVAL(BLIN1,P1,FIXFLU, GRND7)
C   COVAL(BLIN1,U1,0.0   , GRND7)
C   COVAL(BLIN1,V1,0.0   , GRND7)
C   COVAL(BLIN1,W1,0.0   , GRND7)
C   COVAL(BLIN1,KE,0.0   , GRND7)
C   COVAL(BLIN1,EP,0.0   , GRND7)
C   COVAL(BLIN1,TEM1,0.0 , GRND7)
C
C     In addition, the following SPEDAT statements must be set in the Q1
C     input file:
C
C   SPEDAT(SET,BLIN1,VELX,R,Wx)   : x-component of Qref
C   SPEDAT(SET,BLIN1,VELY,R,Wy)   : y-component of Qref
C   SPEDAT(SET,BLIN1,VELZ,R,Wz)   : z-component of Qref
C   SPEDAT(SET,BLIN1,TAIR,R,Tair) : Air temperature
C
C   SPEDAT(SET,BLIN1,VDIR,C,Y)     : vertical coordinate direction
C   SPEDAT(SET,BLIN1,RHOIN,R,1.189): inlet density presumed uniform
C   SPEDAT(SET,BLIN1,BLTY,C,LOGL)  : profile type, i.e. LOGL, POWL or TABLE
C   SPEDAT(SET,BLIN1,ZO,R,2.0E-02) : effective roughness height
C   SPEDAT(SET,BLIN1,REFH,R,10.0)  : ref. height for ref. velocity
C   SPEDAT(SET,BLIN1,ALPHA,R,0.21) : power-law index for power-law profile
C   SPEDAT(SET,BLIN,U-TABLE,C,file_name) : name of velocity profile table for
C                                          tabular profile
C   SPEDAT(SET,BLIN,K-TABLE,C,file_name) : name of turbulent kinetic energy
C                                          profile table for tabular profile
C
C   NOTE: for tabular profiles, it is the user's resposibility to ensure that
C         the tables cover the entire height of the domain. If any cells are
C         found below the lowest point in the table or above the highest point
C         in the table, the first or last values will be used. In between,
C         linear interpolation will be used to obtain values at the local height
C         above the terrain. It is the user's responsibility to ensure that
C         the roughness height used in the fully-rough wall function is consistent
C         with the velocity and turbulent kinetic energy profiles.
C
C         The table should contain a header line (giving the titles of the columns)
C         followed by pairs of values in free format with blank or comma as
C         separators. The first value is the height, the second is the velocity or
C         turbulent kinetic energy. The two tables do not have to contain the same
C         number of data points.
C
C   The PIL variable WALLB is used for the displacement height.
C
C   To calculate external values at fixed pressure boundaries, the COVAL for P1
C   can be replaced with a normal pressure boundary condition:
C   COVAL(BLIN1,P1,COef, 0.0)
C   Note that the external pressure is always kept at zero relative to pRESS), and
C   PRESS0 is updated to the current atmospheric pressure.
C
C   In order for the BLIN patch to automatically detect inflow/outflow based
C   on wind direction and switch between fixed mass flow and fixed pressure
C   the COVAL for P1 should be:
C   COVAL(BLIN1,P1,GRND8, GRND7)
C
C   At the SKY boundary, in addition to the fixed pressure condition, the diffusive
C   vertical link can be introduced by using GRND7 as the COefficient for velocities
C   and turbulence variables. The SKY patch is recognised because it is normal to
C   the vertical coordinate direction.
C   COVAL(BLIN1,P1,   COef,  Pext)
C   COVAL(BLIN1,U1,  GRND7, GRND7)
C   COVAL(BLIN1,V1,  GRND7, GRND7)
C   COVAL(BLIN1,W1,    0.0,   0.0) [the vertical component does not need a COVAL]
C   COVAL(BLIN1,KE,  GRND7, GRND7)
C   COVAL(BLIN1,EP,  GRND7, GRND7)
C   COVAL(BLIN1,TEM1,GRND7, GRND7)
C
C   The transient variation of wind speed, direction and air temperature may be taken
C   from  a weather data file in EPW format. In this case, the name of the data file
C   is passed as:
C     SPEDAT(SET,'BLIN','WEATHERFILE',C,filename)
C
C   The data required for the current run will have been extracted from the weather file
C   by the pre-processor. The number of lines of data sent is set as:
C     SPEDAT(SET,'BLIN','NLINES',I,nlines)
C
C   The number of entries per hour (defaulted to 1) is set by:
C     SPEDAT(SET,'BLIN','NPHOUR',I,nphour)
C
C   Each line of data contains the wind speed, wind direction and temperature
C     SPEDAT(SET,'BLIN','LINEn',C, wspeed,wdir,wtemp)
C
C   The direction of the domain relative to North must also be specified as
C     SPEDAT(SET,'BLIN','AXDIR',R,axdir)
C   where axdir is the angle between Y and North (for up Z)
C
C   As an alternative to using a weather file, multiple WIND objects can be specified
C   in a transient. The start and end times must be set to prevent overlap.
C
C   If STORE(WAMP) appears in the Q1, the wind amplification factor, defined
C   as absolute velocity / wind speed at 1.5m, will be calculated and placed in WAMP.
C   If VABS has been stored, that will be used for the absolute velocity. If it
C   is not stored, the absolute velocity will be calculated locally. The reference
C   height for WAMP can be changed using the following spedat:
C     SPEDAT(SET,'BLIN','WAMPH',R,wamph)
C   If the line is absent, a height of 1.5m will be assumed.
C
C     Limitations:
C     No provision has been made as yet for temperature and scalar
C     profiles, and the logarithmic profile is restricted to the fully-
C     rough wall law, as encountered in the atmospheric boundary layer.
C     At present, this facility cannot be used for BFC=T, GCV=T or
C     CCM=T.
C
C
C------------------------------------------------------------------------
      USE weather_data
      INCLUDE 'farray'
      INCLUDE 'patnos'
      INCLUDE 'parallel'
      INCLUDE 'objnam'
      INCLUDE 'd_earth/parvar'
      INCLUDE 'patcmn'
      INCLUDE 'satear'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      INCLUDE 'grdbfc'
      INCLUDE 'd_earth/pbcstr'
      common/topinf/topid,idx,idy,idz
      common/comlnk/llnk,hlnk,slnk,nlnk,elnk,wlnk
      integer topid,llnk,hlnk,slnk,nlnk,elnk,wlnk
      INCLUDE 'parear'
      common/wdomin/nxsd,nysd,nzsd
      COMMON /INTDMN/IDMN,NUMDMN,NXMAX,NYMAX,NZMAX,NXYMAX,NXYZMX,NCLTOT,
     1               NDOMMX
      COMMON/GENI/NXNY,IGFIL1(8),NFM,IGF(21),IPRL,IBTAU,ILTLS,IGFIL(15),
     1 ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,IRADZ,IVFOL
      COMMON/DRHODP/ITEMP,IDEN/DVMOD/IDVCGR
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      COMMON /FACES/L0FACE,L0FACZ
      COMMON /WORDI1/ NWDS,NCHARS(20),NSEMI,H,NLINES
      INTEGER   H
      COMMON /WORDL1/ ERROR,MORWDS
      COMMON /WORDC1/ WD(20),INLINE
      LOGICAL         ERROR,MORWDS
      CHARACTER       WD*20, INLINE*120,CH1*13,CH2*13,CHGR_13*13,DELIM*1
      COMMON/LRNTM4/RSCMCD,CDDAK2,TAUDKE,RTTDKE,AKC,EWC
      COMMON/RHILO/HI3D,RLO3D
      COMMON/IHILO/IXHI3D,IYHI3D,IZHI3D,IXLO3D,IYLO3D,IZLO3D,IMAXC(150)
      LOGICAL MASKPT,SLD,QEQ,QGT,QLE,QNE, LPTDMN,LSOLID,EQZ,
     1        LWP,SWP,WWP,LTZ,SHOWCOMF,WEICOEF,NEN,LPAR
      CHARACTER*14 CTEMP*8,BLTY*5,VDIR*1,LINE*256,CVAL*68,SITENAME*68
      CHARACTER*8 TYPECOMF,WINDFILE*256,TERRNAM*12,COBNM2*12,
     1            VELFILE*256,KEFILE*256
      PARAMETER (MAXDMN = 10, MAXBLIN=20)
      INTEGER L0IPAT(MAXDMN),L0VELX(MAXDMN),L0VELY(MAXDMN),
     1        L0VELZ(MAXDMN),L0QREF(MAXDMN),L0VDIR(MAXDMN),
     1        L0ZREF(MAXDMN),L0HREF(MAXDMN),L0BLTY(MAXDMN),
     1        L0DENS(MAXDMN),L0POWR(MAXDMN),L0HI(MAXDMN,MAXBLIN),L0H,
     1        NBLIN(MAXDMN), L0SKY(MAXDMN),L0TAIR(MAXDMN),
     1        L0PCOE(MAXDMN),L0PVAL(MAXDMN),L0RELH(MAXDMN),
     1        L0WMP(MAXDMN),L0HIWAF(MAXDMN),L0VAB(MAXDMN)
      REAL, ALLOCATABLE :: BUFF(:)
      REAL, ALLOCATABLE :: WNDA(:,:),VEL_TAB(:,:), KE_TAB(:,:)
      REAL NORML(3)
      LOGICAL*4 EXISTS
      SAVE L0IPAT,L0VELX,L0VELY,L0VELZ,L0QREF,L0VDIR,L0ZREF,L0HREF,
     1     L0BLTY,L0DENS,L0POWR,L0HI,L0H,NBLIN,L0SKY,NL,L0TAIR,
     1     L0PCOE,L0PVAL,L0RELH,IFLO,IHUM,IHUNIT,LBWAMP,LBVABS,
     1     IBLIN,INIBUOY,LBCP,LBPRO,LBDAN,LBNEN,ISECT,NBINS,L0WMP,
     1     L0WAMP,LBWAF,L0HIWAF,IERR1,IERR2,LBVAV,L0VAB,
     1     NLINEV,NLINEK
      SAVE AXDIR, RHOIN,WAMPVR, WMAST,WNDA,UTHR,DTHR,SECTP,VSECT
      SAVE SHOWCOMF,WEICOEF,NEN
      SAVE TYPECOMF
      SAVE VEL_TAB, KE_TAB
c***********************************************************************
c
      IXL=IABS(IXL)
      NAMSUB='GXBLIN'; NAMFUN=' '
C*****************************************************************
C--- GROUP 1. Preliminaries
      IF(IGR==1)  THEN
        IF(ISC==1) THEN
          IF(.NOT.NULLPR) THEN
            CALL WRYT40('GXBLIN of:                          180116')
            CALL WRIT40('GXBLIN of:                          180116')
          ENDIF
          CALL GXMAKE(L0H,NXNY,'HDIS') ! storage for grid node heights
        ELSEIF(ISC==2) THEN
  100     CONTINUE
          IF(IDMN>1) CALL INDXDM ! get indices for current FGV domain
          IF(NUMDMN>MAXDMN) THEN
            CALL WRITBL
            CALL WRITST
            CALL WRIT40('Please increase MAXDMN in GXBLIN')
            CALL WRIT2I('Current ',MAXDMN,'; Needed',NUMDMN)
            CALL WRITST
            CALL SET_ERR(557,'MAXDMN too small in GXBLIN',1)
            RETURN
          ENDIF
C... Count how many BLINs
          NBLIN(IDMN)=0
          DO IP=1,NUMPAT
            IF((NAMPAT(IP)(1:4)=='BLIN'.OR.NAMPAT(IP)(1:4)=='GXBL')
     1                       .AND.LPTDMN(IP)) NBLIN(IDMN)=NBLIN(IDMN)+1
          ENDDO
          IF(NBLIN(IDMN)==0) RETURN ! no BLINs so nothing to do
          IF(NBLIN(IDMN)>MAXBLIN) THEN
            CALL WRITBL
            CALL WRITST
            CALL WRIT40('Please increase MAXBLIN in GXBLIN')
            CALL WRIT2I('Current ',MAXBLIN,'; Needed',NBLIN(IDMN))
            CALL WRITST
            CALL SET_ERR(557,'MAXBLIN too small in GXBLIN',1)
            RETURN
          ENDIF
C... check for humidity and VABS
          IHUM=0; IHUNIT=2; LBVABS=0
          DO IPH=NPHI,16,-1
            IF(NAME(IPH)=='MH2O') THEN
              IHUM=IPH
              CALL GETSDI('BLIN','HUNIT',IHUNIT)
            ELSEIF(NAME(IPH)=='VABS') THEN
              LBVABS=IPH
            ENDIF
          ENDDO
C... Allocate storage for BLIN parameters
          CALL GXMAKE0(L0IPAT(IDMN),NUMPAT,'L0IPAT')
          CALL GXMAKE(L0VELX(IDMN),NBLIN(IDMN),'L0VELX')
          CALL GXMAKE(L0VELY(IDMN),NBLIN(IDMN),'L0VELY')
          CALL GXMAKE(L0VELZ(IDMN),NBLIN(IDMN),'L0VELZ')
          CALL GXMAKE(L0QREF(IDMN),NBLIN(IDMN),'L0QREF')
          CALL GXMAKE(L0TAIR(IDMN),NBLIN(IDMN),'L0TAIR')
          CALL GXMAKE(L0PCOE(IDMN),NBLIN(IDMN),'L0PCOE')
          CALL GXMAKE(L0PVAL(IDMN),NBLIN(IDMN),'L0PVAL')
          CALL GXMAKE(L0RELH(IDMN),NBLIN(IDMN),'L0RELH')
          CALL GXMAKE(L0VDIR(IDMN),NBLIN(IDMN),'L0VDIR')
          CALL GXMAKE(L0ZREF(IDMN),NBLIN(IDMN),'L0ZREF')
          CALL GXMAKE(L0HREF(IDMN),NBLIN(IDMN),'L0HREF')
          CALL GXMAKE(L0BLTY(IDMN),NBLIN(IDMN),'L0BLTY')
          CALL GXMAKE(L0DENS(IDMN),NBLIN(IDMN),'L0DENS')
          CALL GXMAKE(L0POWR(IDMN),NBLIN(IDMN),'L0POWR')
          CALL GXMAKE0(L0SKY(IDMN), NBLIN(IDMN),'L0SKY')
          IF(LBVABS<=0) CALL GXMAKE(L0VAB(IDMN),NX*NY,'VABS')
C... Get Weather data
          EPWNAM=' ';CALL GETSDC('BLIN','WEATHERFILE',EPWNAM)
          NL=0
          IF(EPWNAM/=' '.AND..NOT.STEADY) THEN
            LOCATION=' '; CALL GETSDC('BLIN','LOCATION',LOCATION) ! Location
            CALL GETSDI('BLIN','NLINES',NL)
            NPHOUR=1; CALL GETSDI('BLIN','NPHOUR',NPHOUR)
            IF(NL>0) THEN
              ALLOCATE(WSPD(NL),WDIR(NL),AIRTEMP(NL),ATMPRES(NL),
     1                                             RELHUM(NL),STAT=IERR)
              WSDP=0.; WDIR=0.; AIRTEMP=0.; ATMPRES=0.; RELHUM=0.
              DO I=1,NL
                WRITE(LINE,'(''LINE'',I4)') I
                LL=8; CALL REMSPC(LINE,LL)
                CALL GETSDC('BLIN',LINE(1:LL),CVAL)
                LL=LENGZZ(CVAL)
                CALL SPLTZZ(CVAL,WD,NWDS,NCHARS,LL)
                WDIR(I)=RRDZZZ(1);WSPD(I)=RRDZZZ(2);AIRTEMP(I)=RRDZZZ(3)
                CALL GETSDC('BLINA',LINE(1:LL),CVAL)
                LL=LENGZZ(CVAL)
                CALL SPLTZZ(CVAL,WD,NWDS,NCHARS,LL)
                ATMPRES(I)=RRDZZZ(1); RELHUM(I)=RRDZZZ(2)
              ENDDO
            ENDIF
          ENDIF
          WDIRN=-999.0; CALL GETSDR('BLIN','WDIR',WDIRN)
          AXDIR=0.0; CALL GETSDR('BLIN','AXDIR',AXDIR)
          INIBUOY=1; CALL GETSDI('BLIN','INIBUOY',INIBUOY)
C... Now loop over BLINs, extract parameters and save
          NBLIN(IDMN)=0
          DO IP=1,NUMPAT
            IF(.NOT.LPTDMN(IP)) GO TO 103
C... allow GENTRA particles to pass through WIND/WIND_PROFILE
            IF(NAMPAT(IP)(1:4)=='BLIN'.OR.NAMPAT(IP)(1:4)=='GXBL')
     1                                                              THEN
              NBLIN(IDMN)=NBLIN(IDMN)+1
              F(L0IPAT(IDMN)+IP)=NBLIN(IDMN)
              IF(MASKPT(IP)) THEN
                CTEMP='^'//NAMPAT(IP)
              ELSE
                CTEMP=NAMPAT(IP)
              ENDIF
              CALL GETCV(IP,1,GCO,GVAL)
              IF(QEQ(GCO,-999.0)) CYCLE ! no settings for ground plane
C.. inlet velocity vector at reference height (will be overwritten for tabular input)
              VELX=0.0; VELY=0.0; VELZ=0.0
              CALL GETSDR(CTEMP,'VELX',VELX)
              F(L0VELX(IDMN)+NBLIN(IDMN))=VELX
              CALL GETSDR(CTEMP,'VELY',VELY)
              F(L0VELY(IDMN)+NBLIN(IDMN))=VELY
              CALL GETSDR(CTEMP,'VELZ',VELZ)
              F(L0VELZ(IDMN)+NBLIN(IDMN))=VELZ
C.. compute inlet velocity magnitude at reference height
              QREF = SQRT ( VELX*VELX + VELY*VELY + VELZ*VELZ )
              F(L0QREF(IDMN)+NBLIN(IDMN))=QREF
C... external temperature
              TAIR=20.0; CALL GETSDR(CTEMP,'TAIR',TAIR)
              F(L0TAIR(IDMN)+NBLIN(IDMN))=TAIR
C... pressure coefficient & value
              PCOEF=1000.0; CALL GETSDR(CTEMP,'PCOEF',PCOEF)
              F(L0PCOE(IDMN)+NBLIN(IDMN))=PCOEF
              PEXT=0.0;  CALL GETSDR(CTEMP,'PEXT',PEXT)
C... External pressure always 0 relative to PRESS0, which is set to atmospheric
              F(L0PVAL(IDMN)+NBLIN(IDMN))=PEXT
C... humidity
              IF(IHUM>0) THEN
                HUMIN=0.0; CALL GETSDR(CTEMP,'HUMIN',HUMIN)
                IF(IHUNIT>0) THEN
                  IF(IHUNIT==1) THEN ! convert from humidity ratio
                    HUMIN=1.E-3*HUMIN
                  ELSEIF(IHUNIT==2) THEN ! convert from relative humidity
                    TIN=TAIR+TEMP0-273.    ! convert to deg C
                    CALL GETSDR(CTEMP,'PEXT',PEXT)
                    PIN=PEXT
                    PVAP = 1.E-2*HUMIN*PVH2O(TIN)
C... save Relative humidity for printout
                    GWRAT = 18.015/28.96; RELH=HUMIN
                    HUMIN = GWRAT*PVAP/(PIN-PVAP)
                  ENDIF
                  HUMIN=HUMIN/(1+HUMIN) ! convert to mass fraction
                ENDIF
                F(L0RELH(IDMN)+NBLIN(IDMN))=HUMIN
              ENDIF
C.. vertical coordinate direction
              VDIR='Z'
              CALL GETSDC(CTEMP,'VDIR',VDIR)
              IF(VDIR=='X') THEN
                 F(L0VDIR(IDMN)+NBLIN(IDMN))=1.
              ELSEIF(VDIR=='Y') THEN
                 F(L0VDIR(IDMN)+NBLIN(IDMN))=2.
              ELSEIF(VDIR=='Z') THEN
                 F(L0VDIR(IDMN)+NBLIN(IDMN))=3.
              ENDIF
C.. inlet density
              RHOIN=1.189
              CALL GETSDR(CTEMP,'RHOIN',RHOIN)
              F(L0DENS(IDMN)+NBLIN(IDMN))=RHOIN
C.. profile type = TABLE (3) for tabular,
C                  POWL  (2) for power law &
C                  LOGL  (1) for log law
              BLTY='LOGL'
              CALL GETSDC(CTEMP,'BLTY',BLTY)
              F(L0BLTY(IDMN)+NBLIN(IDMN))=1
              IF(BLTY=='POWL') F(L0BLTY(IDMN)+NBLIN(IDMN))=2
              IF(BLTY=='TABLE') F(L0BLTY(IDMN)+NBLIN(IDMN))=3
C.. effective roughness length
              ZO=0.1
              CALL GETSDR(CTEMP,'ZO',ZO)
              F(L0ZREF(IDMN)+NBLIN(IDMN))=ZO
C.. reference height for wind reference velocity
              REFH=10.0
              CALL GETSDR(CTEMP,'REFH',REFH)
              F(L0HREF(IDMN)+NBLIN(IDMN))=REFH
C
              IF(BLTY=='POWL') THEN
                CALL GETSDR(CTEMP,'ALPHA',ALPHA)
                F(L0POWR(IDMN)+NBLIN(IDMN))=ALPHA
              ENDIF
              CALL GETPAT(IP,IREG,TYP,JXF,JXL,JYF,JYL,JZF,JZL,JTF,JTL)
              IPBC=0
C... Get lower limit of open area for each BLIN
              IF(VDIR=='Z') THEN ! UP direction is Z. Deal with X-Z and Y-Z
                IF(QEQ(TYP,1.0)) THEN  ! CELL, initialisation
                  CALL SUB4(I1,JXF,I2,JXL,J1,JYF,J2,JYL)
                  CALL SUB2(K1,1,K2,NZ)
                  NDIM=NX*NY ! need NXNY edge locations
                  F(L0SKY(IDMN)+NBLIN(IDMN))=-1
                ELSEIF(TYP<4.) THEN ! West/East (TYP= 3 or 2)
                  I1=ITWO(JXF,JXL,TYP>2.) ! JXF if TYP=3, JXL if TYP=2
                  I2=ITWO(JXF,JXL,TYP>2.)
                  CALL SUB4(J1,JYF,J2,JYL,K1,JZF,K2,JZL)
                  NDIM=NY ! need NY edge locations
                ELSEIF(TYP<6.) THEN  ! South/North (TYP= 5 or 4)
                  CALL SUB2(I1,JXF,I2,JXL)
                  J1=ITWO(JYF,JYL,TYP>4.) ! JYF if TYP=5, JYL if TYP=4
                  J2=ITWO(JYF,JYL,TYP>4.)
                  CALL SUB2(K1,JZF,K2,JZL)
                  NDIM=NX ! need NX edge locations
                ELSE ! High/Low (TYP = 6 or 7) for SKY
                  CALL SUB4(I1,JXF,I2,JXL,J1,JYF,J2,JYL)
                  CALL SUB2(K1,1,K2,NZ)
                  NDIM=NX*NY ! need NXNY edge locations
                  F(L0SKY(IDMN)+NBLIN(IDMN))=1
                ENDIF
                CALL GXMAKE0(L0HI(IDMN,NBLIN(IDMN)),NDIM,'L0HI')
                L0DIST=L0F(ZWNZ)
                ISKY=F(L0SKY(IDMN)+NBLIN(IDMN))
C... loop over area of BLIN to find first open cell
                DO IX=I1,I2
                  DO IY=J1,J2
                    DIST=0.0
                    DO IZZ=K1,K2
                      IZZM1=(IZZ-1)*NXNY
                      L0FACZ= L0FACE+IZZM1
                      I=(IX-1)*NY+IY
                      IJKDM=I+IZZM1
                      IF(ISKY==0) THEN
                        IB=ITWO(IX,IY,TYP>3.)
                      ELSE
                        IB=I
                      ENDIF
                      IF(PARSOL) IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
                      IF(IPBC>0.AND.ISKY==0) THEN          ! dealing with cut cel
                        IT1= ISHPB(IDMN,IPBC)
                        DIST = F(IT1+PBC_CTCEN(3)) ! cut-face centre in Z
                        F(L0HI(IDMN,NBLIN(IDMN))+IB)=DIST
                        GO TO 1001 ! next column
                      ELSEIF(.NOT.SLD(I)) THEN ! open cell
                        IF(IZZ>1.AND.(LWP(I).OR.SLD(I-NXNY)))
     1                                              DIST=F(L0DIST+IZZ-1)
                        F(L0HI(IDMN,NBLIN(IDMN))+IB)=DIST
                        GO TO 1001 ! next column
                      ENDIF
                    ENDDO
 1001               CONTINUE
                  ENDDO
                ENDDO
                IF(MIMD.AND.NPROC>1.AND.TYP<6.) THEN ! parallel & E,W,N or S
C... Must send terrain height from lowest (in Z) processors to upper ones.
                  IF(NZSD>1) THEN ! domain is split in Z
                    ITAG=4880; IXSD=IDX+1; IYSD=IDY+1; IZSD=IDZ+1
                    NCELLS=ITWO(I2-I1+1, J2-J1+1,TYP>3.)
                    IF(IZSD0) THEN ! send to above
                        CALL SENDI_FOR(NCELLS,1,ITAG,IDD)
                        CALL SENDR_FOR(BUFF,NCELLS,ITAG,IDD)
                      ENDIF
                    ENDIF
                    DEALLOCATE(BUFF,STAT=IERR)
                  ENDIF
                ENDIF
              ELSEIF(VDIR=='Y') THEN ! UP direction is Y. Deal with Z-Y and X-Y
                IF(QEQ(TYP,1.0)) THEN  ! CELL, initialisation
                  CALL SUB4(I1,JXF,I2,JXL,K1,JZF,K2,JZL)
                  CALL SUB2(J1,1,J2,NY)
                  NDIM=NX*NZ
                  F(L0SKY(IDMN)+NBLIN(IDMN))=-1
                ELSEIF(TYP<4.) THEN ! West/East (TYP= 3 or 2)
                  I1=ITWO(JXF,JXL,TYP>2.)
                  I2=ITWO(JXF,JXL,TYP>2.)
                  CALL SUB4(J1,JYF,J2,JYL,K1,JZF,K2,JZL)
                  NDIM=NZ
                ELSEIF(TYP>5.) THEN          ! Low/High (TYP= 7 or 6)
                  CALL SUB2(I1,JXF,I2,JXL)
                  CALL SUB2(J1,JYF,J2,JYL)
                  K1=ITWO(JZF,JZL,TYP>6.) ! JZF if TYP=7, JZL if TYP=6
                  K2=ITWO(JZF,JZL,TYP>6.)
                  NDIM=NX
                ELSE  ! North/south for SKY
                  CALL SUB4(I1,JXF,I2,JXL,K1,JZF,K2,JZL)
                  CALL SUB2(J1,1,J2,NY)
                  NDIM=NX*NZ
                  F(L0SKY(IDMN)+NBLIN(IDMN))=1
                ENDIF
                CALL GXMAKE0(L0HI(IDMN,NBLIN(IDMN)),NDIM,'L0HI')
                L0DIST=L0F(YV2D)
                ISKY=F(L0SKY(IDMN)+NBLIN(IDMN))
C... loop over area of BLIN to find first open cell
                DO IX=I1,I2
                  DO IZZ=K1,K2
                    IZZM1=(IZZ-1)*NXNY
                    L0FACZ= L0FACE+IZZM1
                    DIST=0.0
                    DO IY=J1,J2
                      I=(IX-1)*NY+IY
                      IJKDM=I+IZZM1
                      IF(ISKY==0) THEN
                        IB=ITWO(IX,IZZ,TYP>5.)
                      ELSE
                        IB=(IX-1)*NZ+IZZ
                      ENDIF
                      IF(PARSOL) IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
                      IF(IPBC>0.AND.ISKY==0) THEN          ! dealing with cut cel
                        IT1= ISHPB(IDMN,IPBC)
                        DIST = F(IT1+PBC_CTCEN(2)) ! cut-face centre in Y
                        F(L0HI(IDMN,NBLIN(IDMN))+IB)=DIST
                        GO TO 1002
                      ELSEIF(.NOT.SLD(I)) THEN ! open cell
                        IF(IY>1.AND.(SWP(I).OR.SLD(I-1)))
     1                                                DIST=F(L0DIST+I-1)
                        F(L0HI(IDMN,NBLIN(IDMN))+IB)=DIST
                        GO TO 1002
                      ENDIF
                    ENDDO
 1002               CONTINUE
                  ENDDO
                ENDDO
                IF(MIMD.AND.NPROC>1.AND.(TYP<4..OR.TYP>5.)) THEN ! parallel & E,W,H or L
C... Must send terrain height from lowest (in Z) processors to upper ones.
                  IF(NYSD>1) THEN ! domain is split in Z
                    ITAG=4880; IXSD=IDX+1; IYSD=IDY+1; IZSD=IDZ+1
                    NCELLS=ITWO(I2-I1+1, K2-K1+1,TYP>5.)
                    IF(IYSD0) THEN ! send to above
                        CALL SENDI_FOR(NCELLS,1,ITAG,IDD)
                        CALL SENDR_FOR(BUFF,NCELLS,ITAG,IDD)
                      ENDIF
                    ENDIF
                    DEALLOCATE(BUFF,STAT=IERR)
                  ENDIF
                ENDIF
              ELSEIF(VDIR=='X') THEN ! UP direction is X. Deal with Z-X and Y-X
                IF(QEQ(TYP,1.0)) THEN  ! CELL, initialisation
                  CALL SUB4(K1,JZF,K2,JZL,J1,JYF,J2,JYL)
                  CALL SUB2(I1,1,I2,NX)
                  NDIM=NZ*NY
                  F(L0SKY(IDMN)+NBLIN(IDMN))=-1
                ELSEIF(TYP>5.) THEN ! Low/High
                  CALL SUB4(I1,JXF,I2,JXL,J1,JYF,J2,JYL)
                  K1=ITWO(JZF,JZL,TYP>6.) ! JZF if TYP=7, JZL if TYP=6
                  K2=ITWO(JZF,JZL,TYP>6.)
                  NDIM=NY
                ELSEIF(TYP>3.) THEN               ! South/North
                  CALL SUB2(I1,JXF,I2,JXL)
                  J1=ITWO(JYF,JYL,TYP>4) ! JYF if TYP=5, JYL if TYP=4
                  J2=ITWO(JYF,JYL,TYP>4)
                  CALL SUB2(K1,JZF,K2,JZL)
                  NDIM=NZ
                ELSE          ! East/West for SKY
                  CALL SUB4(K1,JZF,K2,JZL,J1,JYF,J2,JYL)
                  CALL SUB2(I1,1,I2,NX)
                  NDIM=NZ*NY
                  F(L0SKY(IDMN)+NBLIN(IDMN))=1
                ENDIF
                CALL GXMAKE0(L0HI(IDMN,NBLIN(IDMN)),NDIM,'L0HI')
                L0DIST=L0F(XU2D)
                ISKY=F(L0SKY(IDMN)+NBLIN(IDMN))
C... loop over area of BLIN to find first open cell
                DO IZZ=K1,K2
                  IZZM1=(IZZ-1)*NXNY
                  L0FACZ= L0FACE+IZZM1
                  DO IY=J1,J2
                    DIST=0.0
                    DO IX=I1,I2
                      I=(IX-1)*NY+IY
                      IJKDM=I+IZZM1
                      IF(ISKY==0) THEN
                        IB=ITWO(IY,IZZ,TYP>5.)
                      ELSE
                        IB=(IY-1)*NZ+IZZ
                      ENDIF
                      IF(PARSOL) IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
                      IF(IPBC>0.AND.ISKY==0) THEN          ! dealing with cut cell
                        IT1= ISHPB(IDMN,IPBC)
                        DIST = F(IT1+PBC_CTCEN(1)) ! cut-face centre in X
                        F(L0HI(IDMN,NBLIN(IDMN))+IB)=DIST
                        GO TO 1003
                      ELSEIF(.NOT.SLD(I)) THEN ! open cell
                        IF(IX>1.AND.(WWP(I).OR.SLD(I-NY)))
     1                                               DIST=F(L0DIST+I-NY)
                        F(L0HI(IDMN,NBLIN(IDMN))+IB)=DIST
                        GO TO 1003
                      ENDIF
                    ENDDO
 1003               CONTINUE
                  ENDDO
                ENDDO
                IF(MIMD.AND.NPROC>1.AND.TYP>3.) THEN ! parallel & N, S, H or L
C... Must send terrain height from lowest (in X) processors to upper ones.
                  IF(NXSD>1) THEN ! domain is split in X
                    ITAG=4880; IXSD=IDX+1; IYSD=IDY+1; IZSD=IDZ+1
                    NCELLS=ITWO(J2-J1+1, K2-K1+1,TYP>5.)
                    IF(IXSD0) THEN ! send to above
                        CALL SENDI_FOR(NCELLS,1,ITAG,IDD)
                        CALL SENDR_FOR(BUFF,NCELLS,ITAG,IDD)
                      ENDIF
                    ENDIF
                    DEALLOCATE(BUFF,STAT=IERR)
                  ENDIF
                ENDIF
              ENDIF
            ENDIF
  103       CONTINUE
          ENDDO
C... Wind amplification factor (1)
          SHOWCOMF=.FALSE.; CALL GETSDL('FLAIR','SHOWCOMF',SHOWCOMF)
C... Wind amplification factor (2)
          LBWAF=LBNAME('WAF')
          IF(LBWAF>0.AND.VDIR=='Z') THEN
            CALL GXMAKE0(L0HIWAF(IDMN),NXNY,'L0HIWAF')
            L0DIST=L0F(ZWNZ)
            TERRNAM=' '; IBLK=0
            CALL GETSDC('FLAIR','TERRAIN',TERRNAM) ! get name of parent
            IF(TERRNAM.NE.' ') THEN  ! parent name set, so find patch number
              LPAR=MIMD.AND.NPROC.GT.1
              ILIM=ITWO(GD_NUMPAT,NUMPAT,LPAR) ! for parallel loop over global patches
              DO III=1,ILIM
                IF(LPAR) THEN ! get object name for parallel
                  COBNM2=GD_OBJNAM(III)
                ELSE
                  COBNM2= OBJNAM(III)
                ENDIF
                IF(COBNM2.EQ.TERRNAM) THEN
                  IF(LPAR) THEN
                    IBLK=GD_INDPAT(III,1)
                  ELSE
                    IBLK=III
                  ENDIF
                  EXIT
                ENDIF
              ENDDO
            ENDIF
            DO IX=1,NX
              DO IY=1,NY
                I=(IX-1)*NY+IY
                DIST=0.0
                DO IZZ=1,NZ
                  L0FACZ=L0FACE+(IZZ-1)*NX*NY ! increment index for SLD()
                  IF(IBLK>0) THEN
                    L0PAT=L0PATNO(IDMN)+(IZZ-1)*NXNY ! index for patch-number store
                    CALL GET_SURFACE(IX,IY,IZZ,IBLK,L0PAT,CAREA,
     1                                                    NORML,IPLUS)
                    IF(CAREA>0.0) THEN
                      IZZM1=(IZZ-1)*NXNY; L0FACZ= L0FACE+IZZM1
                      IJKDM=I+IZZM1
                      IF(PARSOL) IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
                      IF(IPBC>0) THEN        ! dealing with cut cel
                        IT1= ISHPB(IDMN,IPBC)
                        DIST = F(IT1+PBC_CTCEN(3)) ! cut-face centre in Z
                        GO TO 1004 ! next column
                      ELSEIF(.NOT.SLD(I)) THEN ! open cell
                        IF(IZZ>1.AND.(LWP(I).OR.SLD(I-NXNY))) THEN
                          DIST=F(L0DIST+IZZ-1)
                        ENDIF
                        GO TO 1004 ! next column
                      ENDIF
                    ENDIF
                  ELSEIF(.NOT.SLD(I)) THEN ! open cell
                    IF(IZZ>1.AND.(LWP(I).OR.SLD(I-NXNY))) THEN
                      DIST=F(L0DIST+IZZ-1)
                    ENDIF
                    GO TO 1004 ! next column
                  ENDIF
                ENDDO
 1004           CONTINUE
                F(L0HIWAF(IDMN)+I)=DIST
              ENDDO
            ENDDO
            IF(MIMD.AND.NPROC>1) THEN ! parallel
C... Must send terrain height from lowest (in Z) processors to upper ones.
              IF(NZSD>1) THEN ! domain is split in Z
                ITAG=4880; IXSD=IDX+1; IYSD=IDY+1; IZSD=IDZ+1
                NCELLS=NXNY
                IF(IZSD0) THEN ! send to above
                    CALL SENDI_FOR(NCELLS,1,ITAG,IDD)
                    CALL SENDR_FOR(BUFF,NCELLS,ITAG,IDD)
                  ENDIF
                ENDIF
                DEALLOCATE(BUFF,STAT=IERR)
              ENDIF
            ENDIF
          ENDIF
C... loop back for next FGV domain
          IF(IDMN1) THEN
            IDMN=1
            CALL INDXDM
          ENDIF
C... Wind comfort
          IF(SHOWCOMF) THEN
            LBPRO=LBNAME('PRO')
            TYPECOMF='NEN8100'; CALL GETSDC('FLAIR','TYPECOMF',TYPECOMF)
            IF(TYPECOMF=='NEN8100') THEN
              LBDAN=LBNAME('PDAN'); LBNEN=LBNAME('NEN')
              IF(LBDAN+LBNEN<=0) CALL SET_ERR(584,
     1             'Error. STORE NEN and/or PDAN missing for NEN8100',1)
              UTHR=5; DTHR=15; CALL GETSDR('FLAIR','DTHR',DTHR)
            ELSE
              UTHR=6; CALL GETSDR('FLAIR','UTHR',UTHR)
            ENDIF
            WEICOEF=.FALSE.; CALL GETSDL('FLAIR','WEICOEF',WEICOEF)
            LU=111
            WINDFILE='NOTSET'
            CALL GETSDCLONG('FLAIR','WINDFILE',WINDFILE) ! name of file
            LL=LENGZZ(WINDFILE); CALL LOWCSE(WINDFILE,-LL)
            IF(MIMD.AND.NPROC>1) THEN ! must copy remote file to local
              CALL SEND_FILE(LU,WINDFILE,IERR)
              IF(IERR.NE.0) GOTO 110
            ENDIF
            OPEN(LU,FILE=WINDFILE,STATUS='OLD',ERR=110)
C
            DHEIGHT=-999.0 ! initialise to silly value
            IF(WEICOEF) THEN ! read Weibull parameters
              NSECT=12; REWIND LU
              READ(LU,'(A)') LINE; SITENAME=LINE   ! site name
              READ(LU,'(A)') LINE; LL=LENGZZ(LINE) ! mast location & height
              CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL); DHEIGHT=RRDZZZ(3)
              READ(LU,'(A)') LINE; LL=LENGZZ(LINE) ! no of sectors
              CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL); NSECT=IRDZZZ(1)
              ALLOCATE(WNDA(3,NSECT),STAT=IERR)
              WNDA=0.0
              READ(LU,'(A)') LINE; LL=LENGZZ(LINE)
              CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL)
              PSUM=0.0
              DO ISEC=1,NSECT ! read & sum probability
                WNDA(1,ISEC)=RRDZZZ(ISEC); PSUM=PSUM+WNDA(1,ISEC)
              ENDDO
              DO ISEC=1,NSECT ! normalise
                WNDA(1,ISEC)=WNDA(1,ISEC)/PSUM
              ENDDO
              DO IB=2,3
                READ(LU,'(A)') LINE; LL=LENGZZ(LINE)
                CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL)
                DO ISEC=1,NSECT
                  WNDA(IB,ISEC)=RRDZZZ(ISEC)
                ENDDO
              ENDDO
            ELSE ! Read probability bins
              LBVAV=LBNAME('VAV')
              NBINS=0; NSECT=12; REWIND LU
  104         CONTINUE      ! read file to count number of bins
              READ(LU,'(A)',END=105) LINE; LL=LENGZZ(LINE)
              IF(LL==0) GO TO 105 ! catch empty lines
              NBINS=NBINS+1; GO TO 104
  105         CONTINUE      ! 1st three lines are headers
              NBINS=NBINS-3
              REWIND(LU)    ! now read file properly
              READ(LU,'(A)') LINE; SITENAME=LINE ! site name
              READ(LU,'(A)') LINE; LL=LENGZZ(LINE) ! mast location & height
              CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL); DHEIGHT=RRDZZZ(3)
              READ(LU,'(A)') LINE; LL=LENGZZ(LINE) ! no of sectors
              CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL); NSECT=IRDZZZ(1)
              ALLOCATE(WNDA(NBINS+1,NSECT+1),STAT=IERR)
              WNDA=0.0 ! initialise bins
              READ(LU,'(A)') LINE; LL=LENGZZ(LINE)
              CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL)
              DO ISEC=1,NSECT
                WNDA(1,ISEC+1)=RRDZZZ(ISEC)
              ENDDO
              PSUM=0.0
              DO ISEC=2,NSECT+1
                PSUM=PSUM+WNDA(1,ISEC)
              ENDDO
              DO ISEC=2,NSECT+1
                WNDA(1,ISEC)=WNDA(1,ISEC)/(PSUM+TINY)
              ENDDO
              DO IB=2,NBINS+1 ! read the probabilities for each bin
                READ(LU,'(A)') LINE
                LL=LENGZZ(LINE)
                CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL)
                DO ISEC=1,NSECT+1
                  WNDA(IB,ISEC)=RRDZZZ(ISEC)
                ENDDO
              ENDDO
              DO ISEC=2,NSECT+1 ! normalise probabilities for each direction
                PSUM=0.0
                DO IB=2,NBINS+1 ! sum probabilities over all bins
                  PSUM=PSUM+WNDA(IB,ISEC)
                ENDDO
                DO IB=2,NBINS+1 ! divide through by sum
                  WNDA(IB,ISEC)=WNDA(IB,ISEC)/PSUM
                ENDDO
              ENDDO
            ENDIF
            IF(DHEIGHT<=0.0) THEN
              CALL WRITST
              CALL WRIT40('ERROR! Mast (measurement) height not found')
              CALL WRITST
              CALL SET_ERR(583,
     1                   'Error. Mast (measurement) height not found',1)
            ENDIF
            SECSIZE=360./NSECT ! sector size
            IF(WDIRN>=360.0) THEN
              WDIRN=WDIRN-360.0
            ELSEIF(WDIRN<0.0) THEN
              WDIRN=360.0-WDIRN
            ENDIF
C... Wind speed sectors are centred on direction, so for 12 sectors of 30deg each,
C... sector 1 will be from -15 to +15
C           2               15 to  45 etc.
            ISECT=MIN(NSECT,INT((WDIRN+SECSIZE/2.)/SECSIZE)+1) ! current sector
            IF(LBVAV>0) THEN ! sector aaverage velocity stored
              SECTP=WNDA(1,ISECT+1); VSECT=0.0
              DO IB=2,NBINS+1
                VI=0.5*(WNDA(IB,1)-WNDA(IB-1,1))+WNDA(IB-1,1)
                PROB=WNDA(IB,ISECT+1)
                VSECT=VSECT+VI*PROB
              ENDDO
              QREF=VSECT; REFH=DHEIGHT
              DO II=1,NBLIN(IDMN)
                F(L0QREF(IDMN)+II)=QREF
                F(L0HREF(IDMN)+II)=REFH
              ENDDO
            ENDIF
          ENDIF
C... Inlet profile tables
          IF(BLTY=='TABLE') THEN
            CLOSE(LU,IOSTAT=IOS); LU=111; VELFILE='NOTSET'
            CALL GETSDCLONG('BLIN','U-TABLE',VELFILE) ! name of velocity file
            LL=LENGZZ(VELFILE); CALL LOWCSE(VELFILE,-LL)
            IF(MIMD.AND.NPROC>1) THEN ! must copy remote file to local
              CALL SEND_FILE(LU,VELFILE,IERR)
              IF(IERR.NE.0) THEN
                WINDFILE=VELFILE; GOTO 110
              ENDIF
            ENDIF
            OPEN(LU,FILE=VELFILE,STATUS='OLD',IOSTAT=IOS)
            IF(IOS/=0) THEN
              WINDFILE=VELFILE; GO TO 110
            ENDIF
            NLINES=0; REWIND LU
 1041       CONTINUE      ! read file to count number of lines
            READ(LU,'(A)',END=1051) LINE; LL=LENGZZ(LINE)
            IF(LL==0) GO TO 1051 ! catch empty lines
            NLINES=NLINES+1; GO TO 1041
 1051       CONTINUE      ! 1st lines is header
            NLINEV=NLINES-1
            ALLOCATE(VEL_TAB(NLINEV,4),STAT=IERR)
            REWIND LU
            ANG=WDIRN-AXDIR
            IF(ANG<0.0) THEN
              ANG=360.0+ANG
C... if wind angle > 360, subtract 360
            ELSEIF(ANG>360) THEN
              ANG=ANG-360.0
            ENDIF
            ANGR=ANG*ATAN(1.)/45. ! make radians
            READ(LU,'(A)') LINE   ! skip header line
            DO IL=1,NLINEV
              READ(LU,'(A)') LINE; LL=LENGZZ(LINE) ! read a line
              CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL)  ! split into words
              IF(NWDS/=2) THEN
                WRITE(LUPR1,'('' ERROR reading '',A)') VELFILE
                CALL WRIT40
     1          (' Velocity profile file must have 2 values per line.')
                CALL WRIT40('Item 1 - Height, Item 2 - velocity')
                CALL SET_ERR(583,
     1          'Error. velocity profile must have 2 values per line',1)
              ENDIF
              VEL_TAB(IL,1)=RRDZZZ(1); VEL=RRDZZZ(2) ! get height and velocity
              IF(VDIR=='Z') THEN     ! Up Z, get velocity components
                VEL_TAB(IL,2)=-VEL*SIN(ANGR) ! U
                VEL_TAB(IL,3)=-VEL*COS(ANGR) ! V
                VEL_TAB(IL,4)= 0.0           ! W
              ELSEIF(VDIR=='Y') THEN ! Up Y, get velocity components
                VEL_TAB(IL,2)=-VEL*COS(ANGR) ! U
                VEL_TAB(IL,3)= 0.0           ! V
                VEL_TAB(IL,4)=-VEL*SIN(ANGR) ! W
              ELSE                  ! Up X, get velocity components
                VEL_TAB(IL,2)= 0.0           ! U
                VEL_TAB(IL,3)=-VEL*SIN(ANGR) ! V
                VEL_TAB(IL,4)=-VEL*COS(ANGR) ! W
              ENDIF
            ENDDO
            QREF=SQRT(VEL_TAB(NLINEV,2)**2+VEL_TAB(NLINEV,3)**2+
     1                VEL_TAB(NLINEV,4)**2)
            DO II=1,NBLIN(IDMN)
              F(L0VELX(IDMN)+II)=VEL_TAB(NLINEV,2) ! fill these to get flow
              F(L0VELY(IDMN)+II)=VEL_TAB(NLINEV,3) ! direction at boundaries
              F(L0VELZ(IDMN)+II)=VEL_TAB(NLINEV,4)
              F(L0QREF(IDMN)+II)=QREF
            ENDDO
            CLOSE(LU,IOSTAT=IOS); LU=111; KEFILE='NOTSET'
            CALL GETSDCLONG('BLIN','K-TABLE',KEFILE) ! name of KE file
            LL=LENGZZ(KEFILE); CALL LOWCSE(KEFILE,-LL)
            IF(MIMD.AND.NPROC>1) THEN ! must copy remote file to local
              CALL SEND_FILE(LU,KEFILE,IERR)
              IF(IERR.NE.0) THEN
                WINDFILE=KEFILE; GOTO 110
              ENDIF
            ENDIF
            OPEN(LU,FILE=KEFILE,STATUS='OLD',IOSTAT=IOS)
            IF(IOS/=0) THEN
              WINDFILE=KEFILE; GO TO 110
            ENDIF
            NLINES=0; REWIND LU
 1042       CONTINUE      ! read file to count number of lines
            READ(LU,'(A)',END=1052) LINE; LL=LENGZZ(LINE)
            IF(LL==0) GO TO 1052 ! catch empty lines
            NLINES=NLINES+1; GO TO 1042
 1052       CONTINUE      ! 1st lines is header
            NLINEK=NLINES-1
            ALLOCATE(KE_TAB(NLINEK,2),STAT=IERR)
            REWIND LU
            READ(LU,'(A)') LINE ! skip header
            DO IL=1,NLINEK
              READ(LU,'(A)') LINE; LL=LENGZZ(LINE) ! read line
              CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL)  ! split into words
              IF(NWDS/=2) THEN
                WRITE(LUPR1,'('' ERROR reading '',A)') KEFILE
                CALL WRIT40
     1           (' K profile file must have 2 values per line.')
                CALL WRIT40('Item 1 - Height, Item 2 - KE')
                CALL SET_ERR(583,
     1           'Error. K profile must have 2 values per line',1)
              ENDIF
              KE_TAB(IL,1)=RRDZZZ(1); KE_TAB(IL,2)=RRDZZZ(2) ! save height & KE
            ENDDO
            CLOSE(LU,IOSTAT=IOS)
          ENDIF
C
          CALL WRITST
          CALL WRIT40(' Wind profile data')
          CALL WRIT40(' -----------------')
          WRITE(LUPR1,'(''  Up direction '',A1)') VDIR
          CH1=CHGR_13(REFH); LT=LENGZZ(CH1)
          IF(BLTY/='TABLE')
     1     WRITE(LUPR1,'(''  Reference height:    '',A,'' m'')')
     1                                                        CH1(1:LT)
          CH1=CHGR_13(ZO); LT=LENGZZ(CH1)
          WRITE(LUPR1,'(''  Roughness height:    '',A,'' m'')')
     1                                                        CH1(1:LT)
          IF(NEZ(WALLB)) THEN
            CH1=CHGR_13(WALLB); LT=LENGZZ(CH1)
            WRITE(LUPR1,'(''  Displacement height: '',A,'' m'')')
     1                                                        CH1(1:LT)
          ENDIF
          IF(BLTY=='POWL') THEN
            WRITE(LUPR1,'(''  Profile type - Power Law'')')
            CH1=CHGR_13(ALPHA)
            WRITE(LUPR1,'(''  Power law constant: '',A)') CH1
          ELSEIF(BLTY=='LOGL') THEN
            WRITE(LUPR1,'(''  Profile type - Logarithmic Law'')')
          ELSE
            WRITE(LUPR1,'(''  Profile type - from tables'')')
            LL=LENGZZ(VELFILE)
            WRITE(LUPR1,'(''  Velocity file used: '',A)') VELFILE(1:LL)
            LL=LENGZZ(KEFILE)
            WRITE(LUPR1,'(''  KE file used:       '',A)') KEFILE(1:LL)
          ENDIF
          CH1=CHGR_13(PCOEF)
          WRITE(LUPR1,'(''  Pressure coeff at outflow boundaries: '',
     1              A)') CH1
          IF(EPWNAM/=' ') THEN ! Weather file specified
            LL=LENGZZ(EPWNAM)
            WRITE(LUPR1,'(''  Using weather data file '',A)')
     1                                                      EPWNAM(1:LL)
            IF(NL==0) THEN
              CH1=CHGR_13(WDIRN); LT=LENGZZ(CH1)
              WRITE(LUPR1,
     1         '(''  Wind direction:      '',A,A)') CH1(1:LT),CHAR(176)
              CH1=CHGR_13(QREF); LT=LENGZZ(CH1)
              IF(BLTY/='TABLE') WRITE(LUPR1,
     1         '(''  Wind speed:          '',A,'' m/s'')')   CH1(1:LT)
              IF(ITEM1>0) THEN
                CH1=CHGR_13(TAIR); LT=LENGZZ(CH1)
                WRITE(LUPR1,
     1         '(''  Air temperature:     '',A,A1,''C'')')   CH1(1:LT),
     1                                                         CHAR(176)
              ENDIF
              CH1=CHGR_13(PEXT); LT=LENGZZ(CH1)
              WRITE(LUPR1,
     1         '(''  External pressure:   '',A,'' Pa'')')    CH1(1:LT)
              IF(IHUM>0) THEN
                CH1=CHGR_13(RELH); LT=LENGZZ(CH1)
                WRITE(LUPR1,
     1         '(''  Relative humidity:   '',A,'' %'')')     CH1(1:LT)
              ENDIF
              CH1=CHGR_13(RHOIN); LT=LENGZZ(CH1)
              WRITE(LUPR1,
     1         '(''  Air density:         '',A,'' kg/m^3'')') CH1(1:LT)
            ELSE
              WRITE(LUPR1,
     1   '(''  Time  Direction     Speed        Temperature'',
     1                       ''    Pressure      Humidity'')')
              DO I=1,NL
                WRITE(LUPR1,
     1     '(I4,'' '',1PE13.6,'' '',1PE13.6,'' '',1PE13.6,'' '',1PE13.6,
     1     '' '',1PE13.6)') I-1, WDIR(I),WSPD(I),AIRTEMP(I),ATMPRES(I),
     1                      RELHUM(I)
              ENDDO
            ENDIF
          ELSEIF(STEADY) THEN ! manual wind settings
            IF(QNE(WDIRN,-999.0)) THEN
              CH1=CHGR_13(WDIRN); LT=LENGZZ(CH1)
              WRITE(LUPR1,
     1         '(''  Wind direction:      '',A,A)') CH1(1:LT),CHAR(176)
            ENDIF
            CH1=CHGR_13(QREF); LT=LENGZZ(CH1)
            IF(BLTY/='TABLE') THEN
              WRITE(LUPR1,
     1         '(''  Wind speed:          '',A,'' m/s'')')    CH1(1:LT)
            ELSE
              CH2=CHGR_13(VEL_TAB(NLINEV,1)); L2=LENGZZ(CH2)
              WRITE(LUPR1,
     1     '(''  Wind speed:          '',A,
     1       '' m/s (at top of table '',A,'' m)'')') CH1(1:LT),CH2(1:L2)
            ENDIF
            IF(ITEM1>0) THEN
              CH1=CHGR_13(TAIR); LT=LENGZZ(CH1)
              WRITE(LUPR1,
     1         '(''  Air temperature:     '',A,A1,''C'')')    CH1(1:LT),
     1                                                         CHAR(176)
            ENDIF
            CH1=CHGR_13(PEXT); LT=LENGZZ(CH1)
            WRITE(LUPR1,
     1         '(''  External pressure:   '',A,'' Pa'')')     CH1(1:LT)
            IF(IHUM>0) THEN
              IF(IHUNIT==0) THEN
                CH1=CHGR_13(HUMIN); LT=LENGZZ(CH1)
                WRITE(LUPR1,
     1         '(''  H2O mass fraction:   '',A,'' '')')       CH1(1:LT)
              ELSEIF(IHUNIT==1) THEN
                CH1=CHGR_13(HUMIN*1000.); LT=LENGZZ(CH1)
                WRITE(LUPR1,
     1         '(''  Humidity ratio:      '',A,'' g/kg'')')   CH1(1:LT)
              ELSE
                CH1=CHGR_13(RELH); LT=LENGZZ(CH1)
                WRITE(LUPR1,
     1         '(''  Relative humidity:   '',A,'' %'')')      CH1(1:LT)
              ENDIF
            ENDIF
            CH1=CHGR_13(RHOIN); LT=LENGZZ(CH1)
            WRITE(LUPR1,
     1         '(''  Air density:         '',A,'' kg/m^3'')') CH1(1:LT)
          ENDIF
C... Pressure coefficient
          LBCP=LBNAME('CP')
          IF(SHOWCOMF) THEN
            CALL WRITBL
            WRITE(14,'(''  Wind comfort input data'')')
            WRITE(14,'(''  -----------------------'')')
            LL=LENGZZ(SITENAME)
            WRITE(14,'(''  Site name:  '',A)') SITENAME(1:LL)
            WRITE(14,'(''  Using data for sector: '',I2)') ISECT
            CH1=CHGR_13(DHEIGHT); L1=LENGZZ(CH1)
            WRITE(14,'(''  Mast height:'',A)') CH1(1:L1)
            IF(LBVAV>0) THEN
              WRITE(14,'(''  ## Multi-sector run ##'')')
              CH1=CHGR_13(SECTP); L1=LENGZZ(CH1)
              CH2=CHGR_13(VSECT); L2=LENGZZ(CH2)
              WRITE(LUPR1,'(''  Probability of wind in sector'',I3,
     1         '' is: '',A)') ISECT,CH1(1:L1)
              WRITE(LUPR1,
     1        '(''  Measured average wind speed at mast height is: '',A,
     1                                             '' m/s'')') CH2(1:L2)
              WRITE(14,'(''  ## Wind speed was set to '',
     1     ''average velocity at mast height: '',A,'' m/s'')') CH2(1:L2)
              CH1=CHGR_13(DHEIGHT); L1=LENGZZ(CH1)
              WRITE(14,'(''  ## Reference height was set to ''
     1                           ''mast height: '',A,'' m'')') CH1(1:L1)
            ENDIF
            IF(BLTY=='POWL') THEN      ! power-law profile
              WMAST=QREF*(DHEIGHT/REFH)**ALPHA
            ELSEIF(BLTY=='LOGL') THEN   ! log profile
              IF(EQZ(WALLB)) THEN
                GHDZO=AMAX1(DHEIGHT/ZO,2.0)
              ELSE
                IF(LTZ(WALLB)) THEN
                  GHDZO=(DHEIGHT-WALLB)/ZO
                ELSE
                  GHDZO=AMAX1((DHEIGHT-WALLB)/ZO,2.0)
                ENDIF
              ENDIF
              WMAST=QREF*LOG(GHDZO)/LOG(REFH/ZO)
            ELSE                       ! tabular profile
              FACT=-999. ! initialize
              GH=DHEIGHT ! height from probability table file
              DO IL=1,NLINEV-1 ! loop over table
                IF(GH>=VEL_TAB(IL,1).AND.GH<=VEL_TAB(IL+1,1)) THEN
                  FACT=(GH             -VEL_TAB(IL,1))/
     1                 (VEL_TAB(IL+1,1)-VEL_TAB(IL,1))
                  ILP1=IL+1; EXIT
                ENDIF
              ENDDO
              IF(QEQ(FACT,-999.)) THEN ! no value was found
                IF(GH>=VEL_TAB(NLINEV,1)) THEN ! above top of table
                  IL=NLINEV ! use top value
                ELSE                           ! below bottom
                  IL=1      ! use bottom value
                ENDIF
                FACT=0.0; ILP1=IL
              ENDIF
              UCOMP=VEL_TAB(IL,2)+FACT*(VEL_TAB(ILP1,2)-VEL_TAB(IL,2))
              VCOMP=VEL_TAB(IL,3)+FACT*(VEL_TAB(ILP1,3)-VEL_TAB(IL,3))
              WCOMP=VEL_TAB(IL,4)+FACT*(VEL_TAB(ILP1,4)-VEL_TAB(IL,4))
              WMAST=SQRT(UCOMP*UCOMP+VCOMP*VCOMP+WCOMP*WCOMP)
            ENDIF
            CH2=CHGR_13(WMAST); L2=LENGZZ(CH2)
            WRITE(LUPR1,
     1         '(''  Wind-profile speed at mast height: '',A,'' m/s'')')
     1                                                        CH2(1:L2)
            IF(TYPECOMF=='NEN8100') THEN
              LINE='Dutch standard NEN8100'
            ELSE
              CH1=CHGR_13(UTHR); LL=LENGZZ(CH1)
              LINE='probability of exceeding '//CH1(1:LL)//' m/s'
            ENDIF
            LL=LENGZZ(LINE)
            WRITE(LUPR1,'(''  Comfort type: '',A)') LINE(1:LL)
          ENDIF
C... Wind amplification factor (1)
          LBWAMP=LBNAME('WAMP'); !LBVABS=LBNAME('VABS')
          IF(LBWAMP>0) THEN ! WAMP is stored
            WAMPH=1.5; CALL GETSDR('BLIN','WAMPH',WAMPH) ! get ref height
            IF(BLTY=='POWL') THEN       ! power-law profile
              WAMPVR=QREF*(WAMPH/REFH)**ALPHA
            ELSEIF(BLTY=='LOGL') THEN    ! log profile
              IF(EQZ(WALLB)) THEN
                GHDZO=AMAX1(WAMPH/ZO,2.0)
              ELSE
                IF(LTZ(WALLB)) THEN
                  GHDZO=(WAMPH-WALLB)/ZO
                ELSE
                  GHDZO=AMAX1((WAMPH-WALLB)/ZO,2.0)
                ENDIF
              ENDIF
              WAMPVR=QREF*LOG(GHDZO)/LOG(REFH/ZO)
            ELSE                       ! tabular profile
              FACT=-999. ! initialize
              GH=WAMPH
              DO IL=1,NLINEV-1 ! loop over table
                IF(GH>=VEL_TAB(IL,1).AND.GH<=VEL_TAB(IL+1,1)) THEN
                  FACT=(GH             -VEL_TAB(IL,1))/
     1                 (VEL_TAB(IL+1,1)-VEL_TAB(IL,1))
                  ILP1=IL+1; EXIT
                ENDIF
              ENDDO
              IF(QEQ(FACT,-999.)) THEN ! no value was found
                IF(GH>=VEL_TAB(NLINEV,1)) THEN ! above top of table
                  IL=NLINEV ! use top value
                ELSE                           ! below bottom
                  IL=1      ! use bottom value
                ENDIF
                FACT=0.0; ILP1=IL
              ENDIF
              UCOMP=VEL_TAB(IL,2)+FACT*(VEL_TAB(ILP1,2)-VEL_TAB(IL,2))
              VCOMP=VEL_TAB(IL,3)+FACT*(VEL_TAB(ILP1,3)-VEL_TAB(IL,3))
              WCOMP=VEL_TAB(IL,4)+FACT*(VEL_TAB(ILP1,4)-VEL_TAB(IL,4))
              WAMPVR=SQRT(UCOMP*UCOMP+VCOMP*VCOMP+WCOMP*WCOMP)
            ENDIF
            CH1=CHGR_13(WAMPH); CH2=CHGR_13(WAMPVR)
            WRITE(LUPR1,
     1       '(''  Wind speed at '',A,''m:  '',A,'' m/s'')')
     1                             CH1(1:LENGZZ(CH1)),CH2(1:LENGZZ(CH2))
          ENDIF
          CALL WRITST
          RETURN
  110     CONTINUE
          LL=LENGZZ(WINDFILE)
          CALL WRIT40('Cannot open wind data file '//WINDFILE(1:LL))
          CALL SET_ERR(583,'Error. Cannot open wind data file',1)
        ELSEIF(ISC==3) THEN
        ENDIF
      ENDIF
C*****************************************************************
C--- GROUP 13. Boundary conditions and special sources
C                                       Index for Coefficient - CO
C                                       Index for Value       - VAL
      IF(IGR==13.AND.ISC==8) THEN
C------------------- SECTION 8 ------------------- coef = GRND7
C... Set diffusive coefficient at SKY
C      CO = density*turbulent viscosity/distance
!!!C         = density*(Q*)*K*z/dz
        IBLIN=NINT(F(L0IPAT(IDMN)+IPAT))
!!!C.. inlet velocity magnitude at reference height
!!!        QREF=F(L0QREF(IDMN)+IBLIN)
!!!C.. reference height for wind reference velocity
!!!        REFH=F(L0HREF(IDMN)+IBLIN)
!!!C.. effective roughness length
!!!        ZO=F(L0ZREF(IDMN)+IBLIN)
C.. vertical coordinate direction
        IVDIR=NINT(F(L0VDIR(IDMN)+IBLIN))
C... store vertical height of grid nodes in L0H
        IF(IVDIR==1) THEN     ! Up X
!!!          CALL FN0(-L0H,XU2D)
          GDH=0.5*F(L0F(DXU2D)+NX)
        ELSEIF(IVDIR==2) THEN ! Up Y
!!!          CALL FN0(-L0H,YV2D)
          GDH=0.5*F(L0F(DYV2D)+NY)
        ELSE                    ! Up Z
!!!          GH=F(L0F(ZWNZ)+IZ)
          GDH=0.5*F(L0F(DZWNZ)+NZ)
!!!          CALL FN1(-L0H,GH)
        ENDIF
        L0CO=L0F(CO)
        L0DEN=L0F(DEN1); L0VIST=L0F(VIST)
!!!C#### JCL/MRM 19.08.14 allow for 'd' constant in log profile
!!!        QTAU   = AK*QREF/(LOG((REFH-WALLB)/ZO)) ! friction velocity
        DO IX=IXF,IXL
          IADD=NY*(IX-1)
          DO IY=IYF,IYL
            I=IY+IADD
!!!            IF(IVDIR==1) THEN     ! up X
!!!              IB=(IY-1)*NZ+IZ
!!!            ELSEIF(IVDIR==2) THEN ! up Y
!!!              IB=(IX-1)*NZ+IZ
!!!            ELSE                    ! up Z
!!!              IB=I
!!!            ENDIF
!!!            F(L0CO+I)=F(L0DEN+I)*QTAU*AK*AMAX1(F(L0H+I)-
!!!     1                                  F(L0HI(IDMN,IBLIN)+IB),0.0)/GDH
            F(L0CO+I)=F(L0DEN+I)*F(L0VIST+I)/(PRT(INDVAR)*GDH)
          ENDDO
        ENDDO
      ELSEIF(IGR==13.AND.ISC==9) THEN
C------------------- SECTION 8 ------------------- coef = GRND8
C... Set COefficient for mass-flow boundaries.
C      Set to FIXFLU if inflow, set to PCOEF if outflow.
        IF(INDVAR==P1) THEN
          IBLIN=NINT(F(L0IPAT(IDMN)+IPAT))
          VELX=F(L0VELX(IDMN)+IBLIN)
          VELY=F(L0VELY(IDMN)+IBLIN)
          VELZ=F(L0VELZ(IDMN)+IBLIN)
C         check cell face to check inflow/outflow
C e=2, w=3, n=4, s=5, h=6, l=7
          IF(INTTYP==2.OR.INTTYP==3) THEN     ! E or W
            VELIN=VELX
          ELSEIF(INTTYP==4.OR.INTTYP==5) THEN ! N or S
            VELIN=VELY
          ELSEIF(INTTYP==6.OR.INTTYP==7) THEN ! H or L
            VELIN=VELZ
          ENDIF
          IF(ABS(VELIN)<=1.0E-6) VELIN=0.0
          IF(INTTYP==2.OR.INTTYP==4.OR.INTTYP==6) THEN
            IFLO=ITWO(0,1,VELIN<0.0)
          ELSE
            IFLO=ITWO(0,1,VELIN>0.0)
          ENDIF
          IF(IFLO==0) THEN ! inflow at this face
            CALL FN1(CO,FIXFLU)
          ELSE               ! outflow, treat as fixed pressure
            CALL FN1(CO,F(L0PCOE(IDMN)+IBLIN))
          ENDIF
        ENDIF
      ELSEIF(IGR==13.AND.ISC==19) THEN
C------------------- SECTION 19 ------------------- value = GRND7
C     power-law form:  Uz=Uh*(z/h)**a
C     log-law   form:  Uz=Uh*log(z/zo)/log(h/zo)
c
C.. inlet velocity vector at reference height
        IBLIN=NINT(F(L0IPAT(IDMN)+IPAT))
        VELX=F(L0VELX(IDMN)+IBLIN)
        VELY=F(L0VELY(IDMN)+IBLIN)
        VELZ=F(L0VELZ(IDMN)+IBLIN)
C.. inlet velocity magnitude at reference height
        QREF=F(L0QREF(IDMN)+IBLIN)
C.. inlet temperature
        TAIR=F(L0TAIR(IDMN)+IBLIN)
C.. inlet humidity
        HUMIN=F(L0RELH(IDMN)+IBLIN)
C.. vertical coordinate direction
        IVDIR=NINT(F(L0VDIR(IDMN)+IBLIN))
C.. inlet density
        RHOIN=F(L0DENS(IDMN)+IBLIN)
C.. profile type = TABLE (3) for table, POWL (2) for power law & = LOGL (1) for log law
        IBLTY=NINT(F(L0BLTY(IDMN)+IBLIN))
C.. effective roughness length
        ZO=F(L0ZREF(IDMN)+IBLIN)
C.. reference height for wind reference velocity
        REFH=F(L0HREF(IDMN)+IBLIN)
C
        IF(IBLTY==2) ALPHA=F(L0POWR(IDMN)+IBLIN)
C
        ISKY=NINT(F(L0SKY(IDMN)+IBLIN))
C... store vertical height of grid nodes in L0H
        IF(IVDIR==1) THEN     ! up X
          IF(ISKY==1) THEN
            CALL FN0(-L0H,XU2D)
          ELSE
            CALL FN0(-L0H,XG2D)
          ENDIF
        ELSEIF(IVDIR==2) THEN ! up Y
          IF(ISKY==1) THEN
            CALL FN0(-L0H,YV2D)
          ELSE
            CALL FN0(-L0H,YG2D)
          ENDIF
        ELSEIF(IVDIR==3) THEN ! up Z
          IF(ISKY==1) THEN
            GH=F(L0F(ZWNZ)+IZ)
          ELSE
            GH=F(L0F(ZGNZ)+IZ)
          ENDIF
          CALL FN1(-L0H,GH)
        ENDIF
C
C... Inlet mass flux or pressure boundary
        IF(INDVAR==P1) THEN
C         check cell face to calculate flux
C e=2, w=3, n=4, s=5, h=6, l=7
          IF(INTTYP==2.OR.INTTYP==3) THEN     ! E or W
            VELIN=VELX
          ELSEIF(INTTYP==4.OR.INTTYP==5) THEN ! N or S
            VELIN=VELY
          ELSEIF(INTTYP==6.OR.INTTYP==7) THEN ! H or L
            VELIN=VELZ
          ENDIF
          IF(ABS(VELIN)<=1.0E-6) VELIN=0.0
          IF(INTTYP==2.OR.INTTYP==4.OR.INTTYP==6) THEN
            IFLO=ITWO(0,1,VELIN<0.0)
          ELSE
            IFLO=ITWO(0,1,VELIN>0.0)
          ENDIF
          IF(IFLO==1) THEN ! fixed pressure boundary
            CALL FN1(VAL,0.0)
            RETURN
          ENDIF
          CALL GETCV(IPAT,INDVAR,GCO,GVAL)
          IF(QEQ(GCO,GRND8)) THEN
            AMULT=1./FIXFLU
          ELSE
            AMULT=1.
          ENDIF
C... Set sign of mass-flux
          IF(INTTYP==2.OR.INTTYP==4.OR.INTTYP==6) THEN ! E, N or H
            IF(VELIN>0) THEN ! +ve vel points out, so mass is -ve
              VELIN=-VELIN
            ELSE                ! -ve vel points in, so mass is +ve
              VELIN=ABS(VELIN)
            ENDIF
          ENDIF
C... Inlet mass flux
          L0VAL=L0F(VAL); IPBC=0
          IF(IBLTY==3) THEN ! Table profile
            IF(INTTYP==2.OR.INTTYP==3) THEN     ! E or W
              IVELIN=2 ! U component
            ELSEIF(INTTYP==4.OR.INTTYP==5) THEN ! N or S
              IVELIN=3 ! V component
            ELSEIF(INTTYP==6.OR.INTTYP==7) THEN ! H or L
              IVELIN=4 ! W component
            ENDIF
            DO IX=IXF,IXL
              IADD=NY*(IX-1)
              DO IY=IYF,IYL
                I=IY+IADD
                IF(.NOT.SLD(I)) THEN
                  IJKDM=I+(IZ-1)*NX*NY
                  IF(PARSOL) IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
                  IF(IPBC>0.AND.ISKY==0) THEN    ! dealing with cut cell
                    IT1= ISHPB(IDMN,IPBC)
                    GH=F(IT1+PBC_SNYDIST)   ! dist from sunny sub-cell cen.
                                            ! to cut-plane
                  ELSE
                    IF(IVDIR==1) THEN     ! up X
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IY,ITYPE==4.OR.ITYPE==5)
                      ELSE
                        IB=(IY-1)*NZ+IZ
                      ENDIF
                    ELSEIF(IVDIR==2) THEN ! up Y
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=(IX-1)*NZ+IZ
                      ENDIF
                    ELSE                    ! up Z
                      IF(ISKY==0) THEN
                        IB=ITWO(IY,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=I
                      ENDIF
                    ENDIF
                    GH=AMAX1(F(L0H+I)-F(L0HI(IDMN,IBLIN)+IB),0.0)
                  ENDIF
                  VELIN=-999. ! initialize
                  DO IL=1,NLINEV-1 ! loop over table
                    IF(GH>=VEL_TAB(IL,1).AND.GH<=VEL_TAB(IL+1,1)) THEN ! linear
                      VELIN=VEL_TAB(IL,IVELIN)+(GH-VEL_TAB(IL,1))      ! interpolation
     1                 *(VEL_TAB(IL+1,IVELIN)-VEL_TAB(IL,IVELIN))/
     1                  (VEL_TAB(IL+1,1)-VEL_TAB(IL,1))
                      EXIT
                    ENDIF
                  ENDDO
                  IF(QEQ(VELIN,-999.)) THEN ! no value was found
                    IF(GH>=VEL_TAB(NLINEV,1)) THEN ! above top of table
                      VELIN=VEL_TAB(NLINEV,IVELIN) ! use top value
                    ELSE                           ! below bottom
                      VELIN=VEL_TAB(1,IVELIN)      ! use bottom value
                    ENDIF
                  ENDIF
C... Set sign of mass-flux
                  IF(INTTYP==2.OR.INTTYP==4.OR.INTTYP==6) THEN ! E, N or H
                    IF(VELIN>0) THEN ! +ve vel points out, so mass is -ve
                      VELIN=-VELIN
                    ELSE             ! -ve vel points in, so mass is +ve
                      VELIN=ABS(VELIN)
                    ENDIF
                  ENDIF
                  F(L0VAL+I)=RHOIN*VELIN*AMULT
                ELSE
                  F(L0VAL+I)=0.0
                ENDIF
              ENDDO
            ENDDO
          ELSEIF(IBLTY==2) THEN ! Power-Law profile
            VELCON=RHOIN*VELIN/REFH**ALPHA
            DO IX=IXF,IXL
              IADD=NY*(IX-1)
              DO IY=IYF,IYL
                I=IY+IADD
                IF(.NOT.SLD(I)) THEN
                  IJKDM=I+(IZ-1)*NX*NY
                  IF(PARSOL) IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
                  IF(IPBC>0) THEN        ! dealing with cut cell
                    IT1= ISHPB(IDMN,IPBC)
                    GH=F(IT1+PBC_SNYDIST)   ! dist from sunny sub-cell cen.
                                            ! to cut-plane
                  ELSE
                    IF(IVDIR==1) THEN     ! up X
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IY,ITYPE==4.OR.ITYPE==5)
                      ELSE
                        IB=(IY-1)*NZ+IZ
                      ENDIF
                    ELSEIF(IVDIR==2) THEN ! up Y
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=(IX-1)*NZ+IZ
                      ENDIF
                    ELSE                    ! up Z
                      IF(ISKY==0) THEN
                        IB=ITWO(IY,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=I
                      ENDIF
                    ENDIF
                    GH=AMAX1(F(L0H+I)-F(L0HI(IDMN,IBLIN)+IB),0.0)
                  ENDIF
C... divide by FIXFLU if CO was GRND8
                  F(L0VAL+I)=VELCON*(GH**ALPHA)*AMULT
                ELSE
                  F(L0VAL+I)=0.0
                ENDIF
              ENDDO
            ENDDO
          ELSEIF(IBLTY==1) THEN
C... log-law profile
            VELCON=RHOIN*VELIN/LOG(REFH/ZO)
            DO IX=IXF,IXL
              IADD=NY*(IX-1)
              DO IY=IYF,IYL
                I=IY+IADD
                IF(.NOT.SLD(I)) THEN
                  IJKDM=I+(IZ-1)*NX*NY
                  IF(PARSOL) IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
                  IF(IPBC>0) THEN        ! dealing with cut cell
                    IT1= ISHPB(IDMN,IPBC)
                    GH=F(IT1+PBC_SNYDIST)   ! dist from sunny sub-cell cen.
                                            ! to cut-plane
                  ELSE
                    IF(IVDIR==1) THEN     ! up X
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IY,ITYPE==4.OR.ITYPE==5)
                      ELSE
                        IB=(IY-1)*NZ+IZ
                      ENDIF
                    ELSEIF(IVDIR==2) THEN ! up Y
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=(IX-1)*NZ+IZ
                      ENDIF
                    ELSE                    ! up Z
                      IF(ISKY==0) THEN
                        IB=ITWO(IY,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=I
                      ENDIF
                    ENDIF
                    GH=AMAX1(F(L0H+I)-F(L0HI(IDMN,IBLIN)+IB),0.0)
                  ENDIF
                  IF(EQZ(WALLB)) THEN
                    GHDZO=AMAX1(GH/ZO,2.0)
                  ELSE
                    IF(LTZ(WALLB)) THEN
                      GHDZO=(GH-WALLB)/ZO
                    ELSE
                      GHDZO=AMAX1((GH-WALLB)/ZO,2.0)
                    ENDIF
                  ENDIF
C... divide by FIXFLU if CO was GRND8
                  F(L0VAL+I)=VELCON*LOG(GHDZO)*AMULT
                ELSE
                  F(L0VAL+I)=0.
                ENDIF
              ENDDO
            ENDDO
          ENDIF
C
C... Inlet velocity components
        ELSEIF(INDVAR==U1.OR.INDVAR==V1.OR.INDVAR==W1) THEN
          L0VAL=L0F(VAL); l0velin=0
          IF(INDVAR==U1) THEN
            VELIN=VELX
            lbvelin=lbname('UIN')
          ELSEIF(INDVAR==V1) THEN
            VELIN=VELY
            lbvelin=lbname('VIN')
          ELSEIF(INDVAR==W1) THEN
            VELIN=VELZ
            lbvelin=lbname('WIN')
          ENDIF
          if(lbvelin>0) l0velin=l0f(lbvelin)
C
          IPBC=0
          IF(IBLTY==3) THEN ! Table profile
            IVELIN=(INDVAR-1)/2+1 ! get component in table, column 2,3 or 4
            DO IX=IXF,IXL
              IADD=NY*(IX-1)
              DO IY=IYF,IYL
                I=IY+IADD
                IF(.NOT.SLD(I)) THEN
                  IJKDM=I+(IZ-1)*NX*NY
                  IF(PARSOL) IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
                  IF(IPBC>0.AND.ISKY==0) THEN          ! dealing with cut cell
                    IT1= ISHPB(IDMN,IPBC)
                    GH=F(IT1+PBC_SNYDIST)   ! dist from sunny sub-cell cen.
                                            ! to cut-plane
                  ELSE
                    IF(IVDIR==1) THEN     ! up X
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IY,ITYPE==4.OR.ITYPE==5)
                      ELSE
                        IB=(IY-1)*NZ+IZ
                      ENDIF
                    ELSEIF(IVDIR==2) THEN ! up Y
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=(IX-1)*NZ+IZ
                      ENDIF
                    ELSE                    ! up Z
                      IF(ISKY==0) THEN
                        IB=ITWO(IY,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=I
                      ENDIF
                    ENDIF
                    GH=AMAX1(F(L0H+I)-F(L0HI(IDMN,IBLIN)+IB),0.0)
                  ENDIF
                  VELIN=-999. ! initalise
                  DO IL=1,NLINEV-1 ! loop over table
                    IF(GH>=VEL_TAB(IL,1).AND.GH<=VEL_TAB(IL+1,1)) THEN ! linear
                      VELIN=VEL_TAB(IL,IVELIN)+(GH-VEL_TAB(IL,1))      ! interpolation
     1                 *(VEL_TAB(IL+1,IVELIN)-VEL_TAB(IL,IVELIN))/
     1                  (VEL_TAB(IL+1,1)-VEL_TAB(IL,1))
                      EXIT
                    ENDIF
                  ENDDO
                  IF(QEQ(VELIN,-999.)) THEN ! no value found
                    IF(GH>=VEL_TAB(NLINEV,1)) THEN ! above top of table
                      VELIN=VEL_TAB(NLINEV,IVELIN) ! use top value
                    ELSE                           ! below bottom of table
                      VELIN=VEL_TAB(1,IVELIN)      ! use bottom value
                    ENDIF
                  ENDIF
                  F(L0VAL+I)=VELIN
                ELSE
                  F(L0VAL+I)=0.0
                ENDIF
                if(l0velin>0) f(l0velin+i)=f(l0val+i)
              ENDDO
            ENDDO
          ELSEIF(IBLTY==2) THEN ! Power-Law profile
            VELCON=VELIN/REFH**ALPHA
            DO IX=IXF,IXL
              IADD=NY*(IX-1)
              DO IY=IYF,IYL
                I=IY+IADD
                IF(.NOT.SLD(I)) THEN
                  IJKDM=I+(IZ-1)*NX*NY
                  IF(PARSOL) IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
                  IF(IPBC>0.AND.ISKY==0) THEN        ! dealing with cut cell
                    IT1= ISHPB(IDMN,IPBC)
                    GH=F(IT1+PBC_SNYDIST)
                  ELSE
                    IF(IVDIR==1) THEN     ! up X
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IY,ITYPE==4.OR.ITYPE==5)
                      ELSE
                        IB=(IY-1)*NZ+IZ
                      ENDIF
                    ELSEIF(IVDIR==2) THEN ! up Y
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=(IX-1)*NZ+IZ
                      ENDIF
                    ELSE                    ! up Z
                      IF(ISKY==0) THEN
                        IB=ITWO(IY,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=I
                      ENDIF
                    ENDIF
                    GH=AMAX1(F(L0H+I)-F(L0HI(IDMN,IBLIN)+IB),0.0)
                  ENDIF
                  F(L0VAL+I)=VELCON*GH**ALPHA
                ELSE
                  F(L0VAL+I)=0.0
                ENDIF
              ENDDO
            ENDDO
C... log-law profile
          ELSEIF(IBLTY==1) THEN
            VELCON=VELIN/LOG((REFH-WALLB)/ZO)
            DO IX=IXF,IXL
              IADD=NY*(IX-1)
              DO IY=IYF,IYL
                I=IY+IADD
                IF(.NOT.SLD(I)) THEN
                  IJKDM=I+(IZ-1)*NX*NY
                  IF(PARSOL) IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
                  IF(IPBC>0.AND.ISKY==0) THEN          ! dealing with cut cell
                    IT1= ISHPB(IDMN,IPBC)
                    GH=F(IT1+PBC_SNYDIST)
                  ELSE
                    IF(IVDIR==1) THEN     ! up X
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IY,ITYPE==4.OR.ITYPE==5)
                      ELSE
                        IB=(IY-1)*NZ+IZ
                      ENDIF
                    ELSEIF(IVDIR==2) THEN ! up Y
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=(IX-1)*NZ+IZ
                      ENDIF
                    ELSE                    ! up Z
                      IF(ISKY==0) THEN
                        IB=ITWO(IY,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=I
                      ENDIF
                    ENDIF
                    GH=AMAX1(F(L0H+I)-F(L0HI(IDMN,IBLIN)+IB),0.0)
                  ENDIF
                  IF(EQZ(WALLB)) THEN
                    GHDZO=AMAX1(GH/ZO,2.0)
                  ELSE
                    IF(LTZ(WALLB)) THEN
                      GHDZO=(GH-WALLB)/ZO
                    ELSE
                      GHDZO=AMAX1((GH-WALLB)/ZO,2.0)
                    ENDIF
                  ENDIF
                  F(L0VAL+I)=VELCON*LOG(GHDZO)
                ELSE
                  F(L0VAL+I)=0.0
                ENDIF
              ENDDO
            ENDDO
          ENDIF
C
C... inlet k and ep values
        ELSEIF(INDVAR==KE.OR.INDVAR==EP) THEN
          L0VAL=L0F(VAL)
          IPBC=0
          IF(IBLTY==1.OR.IBLTY==2) THEN ! log and power law
            QTAU   = AK*QREF/(LOG((REFH-WALLB)/ZO))
            QTAU2  = QTAU*QTAU
C... Taudke=cmucd**0.5
            GKEIN  = QTAU2/TAUDKE
            GEPCON = QTAU2*QTAU/AK
            IF(INDVAR==KE) CALL FN1(VAL,GKEIN)
            IF(INDVAR==EP) THEN
              L0EP=L0F(EP)
              DO IX=IXF,IXL
                IADD=NY*(IX-1)
                DO IY=IYF,IYL
                  I=IY+IADD
                  IF(.NOT.SLD(I)) THEN ! open cell
                    IJKDM=I+(IZ-1)*NX*NY
                    IF(PARSOL) IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
                    IF(IPBC>0.AND.ISKY==0) THEN          ! dealing with cut cell
                      IT1= ISHPB(IDMN,IPBC)
                      GH=F(IT1+PBC_SNYDIST)
                    ELSE
                      IF(IVDIR==1) THEN     ! up X
                        IF(ISKY==0) THEN
                          IB=ITWO(IZ,IY,ITYPE==4.OR.ITYPE==5)
                        ELSE
                          IB=(IY-1)*NZ+IZ
                        ENDIF
                      ELSEIF(IVDIR==2) THEN ! up Y
                        IF(ISKY==0) THEN
                          IB=ITWO(IZ,IX,ITYPE==2.OR.ITYPE==3)
                        ELSE
                          IB=(IX-1)*NZ+IZ
                        ENDIF
                      ELSE                    ! up Z
                        IF(ISKY==0) THEN
                          IB=ITWO(IY,IX,ITYPE==2.OR.ITYPE==3)
                        ELSE
                          IB=I
                        ENDIF
                      ENDIF
                      GH=AMAX1(F(L0H+I)-F(L0HI(IDMN,IBLIN)+IB),0.0)
                    ENDIF
                    IF(LTZ(WALLB)) THEN
                      GZMDH=(GH-WALLB)
                    ELSE
                      GZMDH=AMAX1(2.0*ZO,(GH-WALLB))
                    ENDIF
                    F(L0VAL+I)=GEPCON/GZMDH
                  ELSE
                    F(L0VAL+I)=0.0
                  ENDIF
                ENDDO
              ENDDO
            ENDIF
          ELSE                ! TABLE profile
            lbkein=lbname('KEIN')
            if(lbkein>0) l0kein=l0f(lbkein)
            lbepin=lbname('EPIN')
            if(lbepin>0) l0epin=l0f(lbepin)
            DO IX=IXF,IXL
              IADD=NY*(IX-1)
              DO IY=IYF,IYL
                I=IY+IADD
                IF(.NOT.SLD(I)) THEN ! open cell
                  IJKDM=I+(IZ-1)*NX*NY
                  IF(PARSOL) IPBC = IPBSEQ(IDMN,IJKDM) ! get cut-cell number
                  IF(IPBC>0.AND.ISKY==0) THEN          ! dealing with cut cell
                    IT1= ISHPB(IDMN,IPBC)
                    GH=F(IT1+PBC_SNYDIST)
                  ELSE
                    IF(IVDIR==1) THEN     ! up X
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IY,ITYPE==4.OR.ITYPE==5)
                      ELSE
                        IB=(IY-1)*NZ+IZ
                      ENDIF
                    ELSEIF(IVDIR==2) THEN ! up Y
                      IF(ISKY==0) THEN
                        IB=ITWO(IZ,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=(IX-1)*NZ+IZ
                      ENDIF
                    ELSE                    ! up Z
                      IF(ISKY==0) THEN
                        IB=ITWO(IY,IX,ITYPE==2.OR.ITYPE==3)
                      ELSE
                        IB=I
                      ENDIF
                    ENDIF
                    GH=AMAX1(F(L0H+I)-F(L0HI(IDMN,IBLIN)+IB),0.0)
                  ENDIF
                  GKEIN=-999. ! initialise
                  DO IL=1,NLINEK-1
                    IF(GH>=KE_TAB(IL,1).AND.GH<=KE_TAB(IL+1,1)) THEN ! linear
                      GKEIN=KE_TAB(IL,2)+(GH-KE_TAB(IL,1))           ! interpolation
     1                 *(KE_TAB(IL+1,2)-KE_TAB(IL,2))/
     1                  (KE_TAB(IL+1,1)-KE_TAB(IL,1))
                      EXIT
                    ENDIF
                  ENDDO
                  IF(QEQ(GKEIN,-999.)) THEN ! no value found
                    IF(GH>=KE_TAB(NLINEK,1)) THEN ! above top of table
                      GKEIN=KE_TAB(NLINEK,2)      ! use top value
                    ELSE                          ! below bottom of table
                      GKEIN=KE_TAB(1,2)           ! use bottom value
                    ENDIF
                  ENDIF
                  IF(INDVAR==KE) THEN
                    F(L0VAL+I)=GKEIN
                  ELSE
                    F(L0VAL+I)=CD*GKEIN**1.5/(AK*GH)
                  ENDIF
                ELSE
                  F(L0VAL+I)=0.0
                ENDIF
                if(indvar==ke.and.lbkein>0) f(l0kein+i)=f(l0val+i)
                if(indvar==ep.and.lbepin>0) f(l0epin+i)=f(l0val+i)
              ENDDO
            ENDDO
          ENDIF
C... External air temperature
        ELSEIF(INDVAR==ITEM1) THEN
          CALL FN1(VAL,TAIR)
C... External air humidity
        ELSEIF(INDVAR==IHUM) THEN
          CALL FN1(VAL,HUMIN)
        ENDIF
      ELSEIF(IGR==19.AND.ISC==1) THEN
C... Group 19, Section 1. Start of time step
C... Echo transient variation of wind parameters
        IF(STEADY) RETURN
        IF(ASSOCIATED(WDIR)) THEN
C... A weather file is in use. Scan through data and interpolate
          DO I=1,NL
            TIM1=(I-1)*3600/NPHOUR; TIM2=I*3600/NPHOUR; DELT=TIM2-TIM1
            IF(QGT(TIM,TIM1).AND.QLE(TIM,TIM2)) THEN
              ANG=WDIR(I)+(WDIR(I+1)-WDIR(I))*(TIM-TIM1)/DELT
              VEL=WSPD(I)+(WSPD(I+1)-WSPD(I))*(TIM-TIM1)/DELT
              TAIR=AIRTEMP(I)+(AIRTEMP(I+1)-AIRTEMP(I))*(TIM-TIM1)/DELT
              PAIR=ATMPRES(I)+(ATMPRES(I+1)-ATMPRES(I))*(TIM-TIM1)/DELT
              IF(IHUM>0) THEN
                HUMIN=RELHUM(I)+(RELHUM(I+1)-RELHUM(I))*(TIM-TIM1)/DELT
              ENDIF
              EXIT
            ENDIF
          ENDDO
          ANG=ANG-AXDIR
          IF(ANG<0.0) THEN
            ANG=360.0+ANG
C... if wind angle > 360, subtract 360
          ELSEIF(ANG>360) THEN
            ANG=ANG-360.0
          ENDIF
          ANGR=ANG*ATAN(1.)/45. ! make radians
          IF(IHUM>0) THEN
            IF(IHUNIT>0) THEN
              IF(IHUNIT==1) THEN ! convert from humidity ratio
                HUMIN=1.E-3*HUMIN
              ELSEIF(IHUNIT==2) THEN ! convert from relative humidity
                PVAP = 1.E-2*HUMIN*PVH2O(TAIR)
                GWRAT = 18.015/28.96
                RELH=HUMIN ! save for printout
                HUMIN = GWRAT*PVAP/(PAIR-PVAP)
              ENDIF
              HUMIN=HUMIN/(1+HUMIN)    ! convert to mass fraction
            ENDIF
          ENDIF
C... Reset buoyancy parameters
          IF(INIBUOY==1) THEN
            IF(QEQ(RHO1,GRND5)) THEN ! Ideal Gas used
              IF(NEZ(RHO1A).AND.IHUM==C1) THEN
                GASCON=RHO1A*(1-HUMIN)+RHO1B*HUMIN
              ELSE
                GASCON=1./RHO1B
              ENDIF
              RHOIN=PAIR/(GASCON*(TAIR+TEMP0))
C.... reset reference density for buoyancy to inlet density
              BUOYD=RHOIN
            ELSEIF(RHO1>0.0) THEN ! density is constant
C.... reset reference temperature
              BUOYE=TAIR
            ENDIF
C.... reset PRESS0 to external pressure
            PRESS0=PAIR
          ENDIF
C... Update values for BLIN patches
          DO IDOM=1,NUMDMN
            DO IBL=1,NBLIN(IDOM)
              IUP=F(L0VDIR(IDOM)+IBL)
              IF(IUP==3) THEN     ! Up Z
                VELX=-VEL*SIN(ANGR)
                VELY=-VEL*COS(ANGR)
                VELZ=0.0
              ELSEIF(IUP==2) THEN ! Up Y
                VELX=-VEL*COS(ANGR)
                VELY=0.0
                VELZ=-VEL*SIN(ANGR)
              ELSE                  ! Up X
                VELX=0.0
                VELY=-VEL*SIN(ANGR)
                VELZ=-VEL*COS(ANGR)
              ENDIF
C... Set inlet velocity components
              F(L0VELX(IDOM)+IBL)=VELX
              F(L0VELY(IDOM)+IBL)=VELY
              F(L0VELZ(IDOM)+IBL)=VELZ
C.. Set inlet velocity magnitude at reference height
              F(L0QREF(IDMN)+IBL)=VEL
              F(L0TAIR(IDMN)+IBL)=TAIR
              F(L0PVAL(IDMN)+IBL)=PAIR
              IF(QEQ(RHO1,GRND5)) F(L0DENS(IDMN)+IBL)=RHOIN
              IF(IHUM>0) F(L0RELH(IDMN)+IBL)=HUMIN
            ENDDO
          ENDDO
        ELSE
C... No Weather file - search for first BLIN active on this time step
          DO IP=1,NUMPAT
            IBLIN=NINT(F(L0IPAT(IDMN)+IP))
            IF(IBLIN>0) THEN
              CALL GETPAT(IP,IREG,TYPE,JXF,JXL,JYF,JYL,JZF,JZL,JTF,JTL)
              IF(ISTEP>=JTF.AND.ISTEP<=JTL) THEN ! active this time step
                PAIR=F(L0PVAL(IDMN)+IBLIN); RHOIN=F(L0DENS(IDMN)+IBLIN)
                TAIR=F(L0TAIR(IDMN)+IBLIN); HUMIN=F(L0RELH(IDMN)+IBLIN)
                VEL=F(L0QREF(IDMN)+IBLIN); RELH=-999.0
                CALL GETSDR(NAMPAT(IP),'WDIR',ANG);ANG=ANG-AXDIR
                IF(INIBUOY==1) THEN
                  PRESS0=PAIR   ! update reference pressure
                  IF(QEQ(RHO1,GRND5)) THEN ! Ideal Gas
                    BUOYD=RHOIN ! update reference density
                  ELSEIF(RHO1>0.0) THEN    ! Constant density
                    BUOYE=TAIR  ! update reference temperature
                  ENDIF
                ENDIF
                EXIT
              ENDIF
            ENDIF
          ENDDO
C... No wind boundary found for this step, so nothing to print
          IF(IP>NUMPAT) THEN
            WRITE(14,'(''No wind boundary conditions found at step '',
     1        I6)') ISTEP
            RETURN
          ENDIF
        ENDIF
C... Echo inputs to RESULT
        CALL WRITST
        CH1=CHGR_13(TIM); LT=LENGZZ(CH1)
        WRITE(LUPR1,'(''  Wind profile data for time step '',I6,
     1   '' at time '',A,''s'')') ISTEP, CH1(1:LT)
        IHR=INT(TIM/3600.); TIM1=TIM-IHR*3600.; IMN=INT(TIM1/60.)
        WRITE(LUPR1,'(''  Elapsed time: '',I4,'' hrs '',I2,'' mins'')'
     1                                                        ) IHR,IMN
        CH1=CHGR_13(ANG+AXDIR); LT=LENGZZ(CH1)
        WRITE(LUPR1,
     1   '(''  Wind direction:    '',A,A)') CH1(1:LT),CHAR(176)
        CH1=CHGR_13(VEL); LT=LENGZZ(CH1)
        WRITE(LUPR1,
     1   '(''  Wind speed:        '',A,'' m/s'')')   CH1(1:LT)
        IF(ITEM1>0) THEN
          CH1=CHGR_13(TAIR); LT=LENGZZ(CH1)
          WRITE(LUPR1,
     1   '(''  Air temperature:   '',A,A1,''C'')')   CH1(1:LT),CHAR(176)
        ENDIF
        CH1=CHGR_13(PAIR); LT=LENGZZ(CH1)
        WRITE(LUPR1,
     1   '(''  External pressure: '',A,'' Pa'')')    CH1(1:LT)
        IF(IHUM>0) THEN
          IF(IHUNIT==0) THEN     ! mass fraction
            CH1=CHGR_13(HUMIN); LT=LENGZZ(CH1)
            WRITE(LUPR1,
     1   '(''  H2O mass fraction: '',A,''  '')')     CH1(1:LT)
          ELSEIF(IHUNIT==1) THEN ! humidity ratio
            CH1=CHGR_13(HUMIN*1000.); LT=LENGZZ(CH1)
            WRITE(LUPR1,
     1   '(''  Humidity ratio:    '',A,'' g/kg'')')  CH1(1:LT)
          ELSE                   ! relative humidity
            IF(QEQ(RELH,-999.0)) THEN ! recover from mass fraction
              HUMIN=AMIN1(1.0,AMAX1(0.0,HUMIN))
              IF(QEQ(HUMIN,1.0)) THEN
                HRAT=1./TINY
              ELSE
                HRAT = HUMIN/(1.-HUMIN+1.E-10)
              ENDIF
              PVSAT = PVH2O(TAIR) ! Water Vapour Saturation Pressure (Pa)
              GWRAT = 18.015/28.96
              PVAP = HRAT*PAIR/(GWRAT+HRAT) ! Water vapour Pressure (Pa)
              RELH = AMAX1(0.0, AMIN1(100.0, 100.*PVAP/PVSAT)) ! Relative humidity (%)
            ENDIF
            CH1=CHGR_13(RELH); LT=LENGZZ(CH1)
            WRITE(LUPR1,
     1   '(''  Relative humidity: '',A,'' %'')')  CH1(1:LT)
          ENDIF
        ENDIF
        CH1=CHGR_13(RHOIN); LT=LENGZZ(CH1)
        WRITE(LUPR1,
     1   '(''  Air density:       '',A,'' kg/m^3'')') CH1(1:LT)
        CALL WRITST
      ELSEIF(IGR==19.AND.ISC==6) THEN
C... Group 19, Section 6. End of IZ slab
C... get velocity amplification factor
        IF(LBWAMP>0.OR.LBWAF>0.OR.LBCP>0.OR.SHOWCOMF) THEN
C... Calculate and store Wind Velocity Amplification Factor
C... Factor is defined as Absolute velocity / wind speed
          IF(STEADY.OR.ASSOCIATED(WDIR)) THEN
            IBL=1
          ELSE
            IBL=IBLIN
          ENDIF
          ZO=F(L0ZREF(IDMN)+IBL)
          IBLTY=NINT(F(L0BLTY(IDMN)+IBL))
          QREF=F(L0QREF(IDMN)+IBL)
          REFH=F(L0HREF(IDMN)+IBL)
          IF(IBLTY==2) ALPHA=F(L0POWR(IDMN)+IBL)
          IF(LBWAMP>0) THEN
            L0WAMP=L0F(LBWAMP)  ! index for amplification store
          ENDIF
          IF(LBWAF>0) THEN
            L0WAF=L0F(LBWAF); L0ZG=L0F(ZGNZ)
            lbvref=lbname('VREF'); lbhref=lbname('HREF')
            if(lbvref>0) l0vrf=l0f(lbvref)
            if(lbhref>0) l0hrf=l0f(lbhref)
          ENDIF
          IF(LBVABS>0) THEN ! VABS is stored, so just use it
            L0VABS=L0F(LBVABS)
          ELSE
            L0VABS=L0VAB(IDMN)
            L0U=0; L0V=0; L0W=0; L0WL=0
            IF(NX>1) L0U=L0F(U1)
            IF(NY>1) L0V=L0F(V1)
            IF(NZ>1) THEN
              L0W=L0F(W1); IF(IZ>1) L0WL=L0F(LOW(W1))
            ENDIF
          ENDIF
          DO IX=1,NX
            DO IY=1,NY
              I=(IX-1)*NY+IY
              IF(LBVABS>0) THEN
                VABS=F(L0VABS+I)
              ELSE
                DU=TINY; UVEL=0.; DV=TINY; VVEL=0.; DW=TINY; WVEL=0.
                IF(.NOT.SLD(I)) THEN ! current cell fluid
C... get cell centre velocity as average of E/W, N/S, H/L faces where not blocked
                  IF(L0U>0) THEN
                    IF(IX 1) THEN ! not at IX=1
                      IF(.NOT.SLD(I-NY)) THEN ! West cell fluid
                        UVEL=UVEL+F(L0U+I-NY); DU=DU+1.
                      ENDIF
                    ENDIF
                  ENDIF
                  IF(L0V>0) THEN
                    IF(IY 1) THEN ! not at IY=1
                      IF(.NOT.SLD(I-1)) THEN ! South cell fluid
                        VVEL=VVEL+F(L0V+I-1); DV=DV+1.
                      ENDIF
                    ENDIF
                  ENDIF
                  IF(L0W>0) THEN
                    IF(IZ 1) THEN ! not at IZ=1
                      IF(.NOT.SLD(I-NXNY)) THEN ! Low cell fluid
                        WVEL=WVEL+F(L0WL+I); DW=DW+1.
                      ENDIF
                    ENDIF
                    UVEL=UVEL/DU; VVEL=VVEL/DV; WVEL=WVEL/DW ! average front/back
                  ENDIF
                ENDIF
C... now get absolute velocity and then amplification factor
                VABS=(UVEL**2+VVEL**2+WVEL**2)**.5
                F(L0VABS+I)=VABS
              ENDIF
              IF(LBWAMP>0) F(L0WAMP+I)=VABS/(WAMPVR+TINY)
              IF(LBWAF>0) THEN
                GH=F(L0ZG+IZ)-F(L0HIWAF(IDMN)+I)
                IF(IBLTY==2) THEN ! power-law profile
                  WAFVR=QREF*(GH/REFH)**ALPHA
                ELSEIF(IBLTY==1) THEN    ! log profile
                  IF(EQZ(WALLB)) THEN
                    GHDZO=AMAX1(GH/ZO,2.0)
                  ELSE
                    IF(LTZ(WALLB)) THEN
                      GHDZO=(GH-WALLB)/ZO
                    ELSE
                      GHDZO=AMAX1((GH-WALLB)/ZO,2.0)
                    ENDIF
                  ENDIF
                  WAFVR=QREF*LOG(GHDZO)/LOG(REFH/ZO)
                ELSE
                  FACT=-999. ! initialize
                  DO IL=1,NLINEV-1 ! loop over table
                    IF(GH>=VEL_TAB(IL,1).AND.GH<=VEL_TAB(IL+1,1)) THEN
                      FACT=(GH             -VEL_TAB(IL,1))/
     1                     (VEL_TAB(IL+1,1)-VEL_TAB(IL,1))
                      ILP1=IL+1; EXIT
                    ENDIF
                  ENDDO
                  IF(QEQ(FACT,-999.)) THEN ! no value was found
                    IF(GH>=VEL_TAB(NLINEV,1)) THEN ! above top of table
                      IL=NLINEV ! use top value
                    ELSE                           ! below bottom
                      IL=1      ! use bottom value
                    ENDIF
                    FACT=0.0; ILP1=IL
                  ENDIF
                  UCOMP=VEL_TAB(IL,2)+FACT*(VEL_TAB(ILP1,2)-
     1                                                    VEL_TAB(IL,2))
                  VCOMP=VEL_TAB(IL,3)+FACT*(VEL_TAB(ILP1,3)-
     1                                                    VEL_TAB(IL,3))
                  WCOMP=VEL_TAB(IL,4)+FACT*(VEL_TAB(ILP1,4)-
     1                                                    VEL_TAB(IL,4))
                  WAFVR=SQRT(UCOMP*UCOMP+VCOMP*VCOMP+WCOMP*WCOMP)
                ENDIF
                F(L0WAF+I)=VABS/(WAFVR+TINY)
                if(lbvref>0) f(l0vrf+i)=WAFvr
                if(lbhref>0) f(l0hrf+i)=gh
              ENDIF
            ENDDO
          ENDDO
C... Pressure coefficient
          IF(LBCP>0) THEN
            QREFSQ=F(L0QREF(IDMN)+IBL)**2 ! current wind speed
            L0CP=L0F(LBCP); L0P1=L0F(P1); L0D1=L0F(DEN1)
            DO I=1,NXNY
              IF(.NOT.SLD(I)) THEN
                F(L0CP+I)=F(L0P1+I)/(0.5*F(L0D1+I)*QREFSQ)
              ELSE
                F(L0CP+I)=0.0
              ENDIF
            ENDDO
          ENDIF
        ENDIF
c... Wind comfort
        IF(SHOWCOMF) THEN
          IF(LBWAMP>0) L0WAMP=L0F(LBWAMP)
          IF(LBVAV>0) L0VAV=L0F(LBVAV)
          lbwrf=lbname('WRF')
          if(lbwrf>0) l0wrf=l0f(lbwrf)
          L0PRO=L0F(LBPRO)
          NEN=TYPECOMF=='NEN8100'
          IF(NEN) THEN
            L0DAN=L0F(LBDAN); L0NEN=L0F(LBNEN)
            UTHR=5.0; DTHR=15.0
          ENDIF
          IF(WEICOEF) THEN  ! Wind data in Weibul coefficients
            AW=WNDA(2,ISECT); AKW=WNDA(3,ISECT)
            DO I=1,NXNY
              PRO=0.0; IF(NEN) PDAN=0.0
              IF(.NOT.SLD(I)) THEN
                WRF=F(L0VABS+I)/WMAST; VINC=UTHR/WRF
                PRO=1.-(1.-EXP(-(VINC/AW)**AKW))
                IF(NEN) THEN
                  DINC=DTHR/WRF; PDAN=1.-(1.-EXP(-(DINC/AW)**AKW))
                ENDIF
              ENDIF
              F(L0PRO+I)=PRO; IF(NEN) F(L0DAN+I)=PDAN
            ENDDO
          ELSE ! Wind data in probability form
            DO I=1,NXNY
              PRO=0.0; IF(NEN) PDAN=0.0
              IF(.NOT.SLD(I)) THEN
                WRF=F(L0VABS+I)/WMAST
                if(lbwrf>0) f(l0wrf+i)=wrf
                DO IB=2,NBINS+1
                  VI=0.5*(WNDA(IB,1)-WNDA(IB-1,1))+WNDA(IB-1,1)
                  PROB=WNDA(IB,ISECT+1)
                  IF(VI*WRF>=UTHR) PRO=PRO+PROB
                  IF(NEN) THEN; IF(VI*WRF>=DTHR) PDAN=PDAN+PROB; ENDIF
                ENDDO
                F(L0PRO+I)=PRO; IF(NEN) F(L0DAN+I)=PDAN
                IF(LBVAV>0) THEN
                  F(L0VAV+I)=WRF*VSECT*SECTP
                  F(L0PRO+I)=F(L0PRO+I)*SECTP
                  IF(NEN) F(L0DAN+I)=F(L0DAN+I)*SECTP
                ENDIF
              ENDIF
            ENDDO
          ENDIF
          IF(NEN) THEN
            DO I=1,NXNY
              PRO=F(L0PRO+I); PDAN=F(L0DAN+I)
              IF(PRO<0.025) THEN
                NEN=1
              ELSEIF(PRO>=0.025.AND.PRO<0.05) THEN
                NEN=2
              ELSEIF(PRO>=0.05.AND.PRO<0.1) THEN
                NEN=3
              ELSEIF(PRO>=0.1.AND.PRO<0.2) THEN
                NEN=4
              ELSEIF(PRO>=0.2) THEN
                NEN=5
              ENDIF
              IF(PDAN>=0.05.AND.PDAN<=0.3) THEN
                NEN=6
              ELSEIF(PDAN>0.3) THEN
                NEN=7
              ENDIF
              F(L0NEN+I)=NEN
            ENDDO
          ENDIF
        ENDIF
      ELSEIF(IGR==19.AND.ISC==7) THEN
C... Group 19, Section 7. End of sweep
        IF(LBWAF>0.AND.IERR1==0) THEN ! IERR1 is error flag from setting InForm below
          CALL HILO3D(LBWAF) ! get domain HI and LO values
          IF(MIMD.AND.NPROC>1) THEN
            CALL PAR_MAXR(HI3D); CALL PAR_MINR(RLO3D)
          ENDIF
          CALL SET_INF('WAFM',HI3D,IERR1)
          IHI=(IXHI3D-1)*NY+IYHI3D
          CALL SET_INF('XWFM',F(L0F(XG2D)+IHI),IERR)
          CALL SET_INF('YWFM',F(L0F(YG2D)+IHI),IERR)
          CALL SET_INF('ZWFM',F(L0F(ZGNZ)+IZHI3D),IERR)
        ENDIF
        IF(LBCP>0.AND.IERR2==0) THEN
          CALL HILO3D(LBCP) ! get domain HI and LO values
          IF(MIMD.AND.NPROC>1) THEN
            CALL PAR_MAXR(HI3D); CALL PAR_MINR(RLO3D)
          ENDIF
          CALL SET_INF('CPM',HI3D,IERR2)
          IHI=(IXHI3D-1)*NY+IYHI3D
          CALL SET_INF('XCPM',F(L0F(XG2D)+IHI),IERR)
          CALL SET_INF('YCPM',F(L0F(YG2D)+IHI),IERR)
          CALL SET_INF('ZCPM',F(L0F(ZGNZ)+IZHI3D),IERR)
        ENDIF
      ELSEIF(IGR==19.AND.ISC==8) THEN
C... Group 19, Section 8. End of time step
        IF(STEADY.OR.ISTEP==LSTEP) THEN
          IF(ALLOCATED(WNDA)) DEALLOCATE(WNDA,STAT=IERR)
          IF(ALLOCATED(VEL_TAB)) DEALLOCATE(VEL_TAB,STAT=IERR)
          IF(ALLOCATED(KE_TAB)) DEALLOCATE(KE_TAB,STAT=IERR)
        ENDIF
      ENDIF
      NAMSUB='gxblin'
      END
C---------------------------------------------------------------------
      SUBROUTINE SEND_FILE(LU,FILNAM,IERR)
      INCLUDE 'parear'
      CHARACTER FILNAM*(*),CVAL*256, DELIM*1
      LOGICAL*4 EXISTS
      OPEN(LU,FILE=FILNAM,STATUS='OLD',IOSTAT=IERR)
      CALL GLSUMI(IERR)
      IF(IERR==0) THEN
        CLOSE(LU); RETURN ! file exists on all processors, no need to copy
      ENDIF
      CALL GETSYSID(ISYSID)
      DELIM='/'; IF(ISYSID.LE.2) DELIM=CHAR(92)
      II=LAST_INDEX(FILNAM,DELIM) ! find last delimiter
      IF(II>0) THEN ! there was one, remove local copies
        CVAL=FILNAM(II+1:) ! name without path i.e. local name
        OPEN(90,FILE=CVAL,STATUS='OLD',ERR=102) ! does it exist
        CLOSE(90,STATUS='DELETE',IOSTAT=IOS)    ! delete if it did
  102   CONTINUE
        EXISTS=.FALSE. ! flag to copy to master as well as slaves
      ELSE ! no delimiter, file is already local
        CVAL=FILNAM
        IF(MYID>0) THEN ! delete copies on slaves
          OPEN(90,FILE=CVAL,STATUS='OLD',ERR=1022) ! does it exist
          CLOSE(90,STATUS='DELETE',IOSTAT=IOS)    ! delete if it did
 1022     CONTINUE
        ENDIF
        EXISTS=.TRUE. ! flag to only copy to slaves
      ENDIF
      II=1 ! now copy central file to local
      CALL COPYMOFFILE(II,DELIM,CVAL,FILNAM,IERR,IOS,EXISTS)
      END