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 (or k-f) equations so as to
C     define a boundary-layer profile at domain boundaries. At present an
C     option is provided only for a neutral atmospheric boundary layer
C     (wind profile) using either a logarithmic or power-law velocity
C     profile, as follows:
C
C          Q = {Q*/K}*ln((z-d)/zo)  or   Q = Qref*(z/zref)^a
C
C          k = Q*^2/sqrt(cmucd)     or  f =  Q*/(K*sqrt(cmucd)*(z-d))
C
C          e = 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) or COVAL(BLIN1,OMEG,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)  or COVAL(BLIN1,OMEG,  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),
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/TSKEMI/ LBKP,LBKT,LBET,LBVOSQ,LBOMEG
COMMON/LRNTM4/RSCMCD,CDDAK2,TAUDKE,RTTDKE,AKC,EWC
COMMON/RHILO/HI3D,RLO3D
COMMON/IHILO/IXHI3D,IYHI3D,IZHI3D,IXLO3D,IYLO3D,IZLO3D,IMAXC(150)
1        LWP,SWP,WWP,LTZ,SHOWCOMF,WEICOEF,NEN,LPAR,LAWS
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)
REAL, ALLOCATABLE :: LTHRESH(:), LUTHR(:),LPRO(:,:)
LOGICAL*4 EXISTS
COMMON/IPRB/LBPRB1,LBPRB2,LBPRB3,LBPRB4,LBPRB5,LBPRB6,
1            L0PRB1,L0PRB2,L0PRB3,L0PRB4,L0PRB5,L0PRB6
INTEGER LBPRB(6),L0PRB(6)
EQUIVALENCE(LBPRB(1),LBPRB1),(L0PRB(1),L0PRB1)
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,LBLAWS,NCRIT
SAVE AXDIR, RHOIN,WAMPVR, WMAST,WNDA,UTHR,DTHR,SECTP,VSECT
SAVE LTHRESH,LUTHR,LPRO
SAVE SHOWCOMF,WEICOEF,NEN,LAWS
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:                          060217')
CALL WRIT40('GXBLIN of:                          060217')
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
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') THEN
NBLIN(IDMN)=NBLIN(IDMN)+1
F(L0IPAT(IDMN)+IP)=NBLIN(IDMN)
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
ALPHA=1./7.; 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
I1=JXF; I2=JXL; J1=JYF; J2=JYL; 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.)
J1=JYF; J2=JYL; K1=JZF; K2=JZL
NDIM=NY ! need NY edge locations
ELSEIF(TYP<6.) THEN  ! South/North (TYP= 5 or 4)
I1=JXF; I2=JXL
J1=ITWO(JYF,JYL,TYP>4.) ! JYF if TYP=5, JYL if TYP=4
J2=ITWO(JYF,JYL,TYP>4.)
K1=JZF; K2=JZL
NDIM=NX ! need NX edge locations
ELSE ! High/Low (TYP = 6 or 7) for SKY
I1=JXF; I2=JXL; J1=JYF; J2=JYL; 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
C... if domain is split in Z and UP direction is Z, must propagate
C    terrain height to all processors at higher (Z) level
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
I1=JXF; I2=JXL; K1=JZF; K2=JZL; 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.)
J1=JYF; J2=JYL; K1=JZF; K2=JZL
NDIM=NZ
ELSEIF(TYP>5.) THEN          ! Low/High (TYP= 7 or 6)
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=NX
ELSE  ! North/south for SKY
I1=JXF; I2=JXL; K1=JZF; K2=JZL; 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
C... if domain is split in Y and UP direction is Y, must propagate
C    terrain height to all processors at higher (Y) level
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
K1=JZF; K2=JZL; J1=JYF; J2=JYL; I1=1; I2=NX
NDIM=NZ*NY
F(L0SKY(IDMN)+NBLIN(IDMN))=-1
ELSEIF(TYP>5.) THEN ! Low/High
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
I1=JXF; I2=JXL
J1=ITWO(JYF,JYL,TYP>4) ! JYF if TYP=5, JYL if TYP=4
J2=ITWO(JYF,JYL,TYP>4)
K1=JZF; K2=JZL
NDIM=NZ
ELSE          ! East/West for SKY
K1=JZF; K2=JZL; J1=JYF; J2=JYL; 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
C... if domain is split in X and UP direction is X, must propagate
C    terrain height to all processors at higher (X) level
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 WAF (2)
LBWAF=LBNAME('WAF') !  WAF = Vabs/(profile speed at local height)
IF(LBWAF>0.AND.VDIR=='Z') THEN
CALL GXMAKE0(L0HIWAF(IDMN),NXNY,'L0HIWAF') ! store for height of surface
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 ! scan up column at IX,IY
L0FACZ=L0FACE+(IZZ-1)*NX*NY ! increment index for SLD()
IF(IBLK>0) THEN ! have named terrain object
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 ! end IZ loop
1004           CONTINUE
F(L0HIWAF(IDMN)+I)=DIST ! store surface height
ENDDO ! end IY loop
ENDDO   ! end IX loop
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 ! end of split-in-Z block
ENDIF ! end of parallel block
ENDIF ! end of WAF block
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)
ELSEIF(TYPECOMF=='LAWSON') THEN
NCRIT=5; CALL GETSDI('FLAIR','NCRIT',NCRIT)
ALLOCATE(LTHRESH(NCRIT),LUTHR(NCRIT),LPRO(NCRIT,NX*NY),
1                                                        STAT=IERR)
LTHRESH(1)=0.06*100.; LUTHR(1)=10.95 ! B5
LTHRESH(2)=0.02*100.; LUTHR(2)=10.95 ! B5
LTHRESH(3)=0.04*100.; LUTHR(3)=8.25 ! B4
LTHRESH(4)=0.06*100.; LUTHR(4)=5.6 ! B3
LTHRESH(5)=0.01*100.; LUTHR(5)=5.6 ! B3
DO IC=1,NCRIT
WRITE(LINE,'(''THRSH'',I1)') IC
CALL GETSDR('FLAIR',LINE(1:6),LTHRESH(IC))
LTHRESH(IC)=LTHRESH(IC)/100. ! convert back from %
WRITE(LINE,'(''UTHR'',I1)') IC
CALL GETSDR('FLAIR',LINE(1:5),LUTHR(IC))
WRITE(LINE,'(''PRO'',I1)') IC
LBPRB(IC)=LBNAME(LINE(1:4))
ENDDO
LBLAWS=LBNAME('LAWS')
IF(LBLAWS<=0) CALL SET_ERR(584,
1             'Error. STORE LAWS missing for Lawson',1)
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
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
CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL)
DO ISEC=1,NSECT
WNDA(IB,ISEC)=RRDZZZ(ISEC)
ENDDO
ENDDO
LBVAV=LBNAME('VAV')
NBINS=0; NSECT=12; REWIND LU
104         CONTINUE      ! read file to count number of bins
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
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
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 WRITST
CALL SET_ERR(583,
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
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
DO IL=1,NLINEV
CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL)  ! split into words
IF(NWDS/=2) THEN
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
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
DO IL=1,NLINEK
CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LL)  ! split into words
IF(NWDS/=2) THEN
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'
ELSEIF(TYPECOMF=='LAWSON') THEN
LINE='Lawson Comfort Criteria '
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)
IF(TYPECOMF=='LAWSON') THEN
WRITE(LUPR1,
1      '(''   Level Probability (%) Max Speed (m/s)'')')
DO IC=1,NCRIT
CH1=CHGR_13(LTHRESH(IC)*100.)
CH2=CHGR_13(LUTHR(IC))
WRITE(LUPR1,'(5X,I1,7X,A,1X,A)') IC,CH1(1:13),CH2(1:13)
ENDDO
ENDIF
ENDIF ! end of SHOWCOMF section
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 (used for WAMP)'')')
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
DO IY=IYF,IYL
!!!            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
DO IY=IYF,IYL
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
DO IY=IYF,IYL
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
DO IY=IYF,IYL
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
DO IY=IYF,IYL
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
DO IY=IYF,IYL
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
DO IY=IYF,IYL
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.OR.INDVAR.EQ.LBOMEG) 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  ! Nb: Taudke=sqrt(CmuCd)
GEPCON = QTAU2*QTAU/AK
IF(INDVAR==KE) CALL FN1(VAL,GKEIN)
IF(INDVAR==EP.OR.INDVAR.EQ.LBOMEG) THEN
IF(INDVAR==LBOMEG) THEN
L0OMEG=L0F(LBOMEG)
GOMCON=QTAU/(AK*TAUDKE)
ELSE
L0EP=L0F(EP)
ENDIF
DO IX=IXF,IXL
DO IY=IYF,IYL
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
IF(INDVAR==EP) THEN
F(L0VAL+I)=GEPCON/GZMDH     ! e = ustar^3/(ak*(z-d))
ELSE IF(INDVAR==LBOMEG) THEN
F(L0VAL+I)=GOMCON/GZMDH     ! f = ustar/(sqrt(CmuCD)*ak*(z-d))
ENDIF
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)
lbomin=lbname('OMIN')
if(lbomin>0) l0omin=l0f(lbomin)
DO IX=IXF,IXL
DO IY=IYF,IYL
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
ELSEIF (INDVAR==EP) THEN             ! Nb: Cd = (CmuCd)^0.75
F(L0VAL+I)=CD*GKEIN**1.5/(AK*GH)   ! e = Cd*k^1.5/(ak*(z-d))
ELSEIF (INDVAR==LBOMEG) THEN            ! Nb: RTTDKE=SQRT(TAUDKE) & TAUDKE=SQRT(CM
F(L0VAL+I)=SQRT(GKEIN)/(AK*RTTDKE*GH) ! f = sqrt(k)/(ak*(z-d)*(CmuCd)^0.25)
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)
if(indvar==lbomeg.and.lbomin>0) f(l0omin+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(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
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
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)
IF(LBPRO>0) L0PRO=L0F(LBPRO)
NEN=TYPECOMF=='NEN8100'; LAWS=TYPECOMF=='LAWSON'
IF(NEN) THEN
L0DAN=L0F(LBDAN); L0NEN=L0F(LBNEN)
UTHR=5.0; DTHR=15.0
ELSEIF(LAWS) THEN
L0LAWS=L0F(LBLAWS)
DO IC=1,6
IF(LBPRB(IC)>0) L0PRB(IC)=L0F(LBPRB(IC))
ENDDO
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
IF(LAWS) THEN ! Lawson criterion
DO IC=1,NCRIT
VINC=LUTHR(IC)/WRF
LPRO(IC,I)=1.-(1.-EXP(-(VINC/AW)**AKW))
ENDDO
ELSE ! probability of exceeding or NEN
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
ENDIF ! end of .NOT.SLD
IF(.NOT.LAWS) F(L0PRO+I)=PRO; IF(NEN) F(L0DAN+I)=PDAN
ENDDO
ELSE ! Wind data in probability table form
IF(LAWS) LPRO=0.0
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 ! loop over wind speed bins
VI=0.5*(WNDA(IB,1)-WNDA(IB-1,1))+WNDA(IB-1,1) ! mid-bin speed
PROB=WNDA(IB,ISECT+1) ! bin probability
IF(LAWS) THEN ! Lawson
DO IC=1,NCRIT ! loop over criteria
IF(VI*WRF>=LUTHR(IC)) LPRO(IC,I)=LPRO(IC,I)+PROB ! sum probability
if(lbprb(ic)>0) f(l0prb(ic)+i)=LPRO(IC,I)
ENDDO
ELSE ! Probability or NEN
IF(VI*WRF>=UTHR) PRO=PRO+PROB ! sum probability
IF(NEN) THEN; IF(VI*WRF>=DTHR) PDAN=PDAN+PROB; ENDIF
ENDIF
ENDDO ! end loop ovr windspeed bins
IF(.NOT.LAWS) F(L0PRO+I)=PRO; IF(NEN) F(L0DAN+I)=PDAN
IF(LBVAV>0) THEN ! Sector velocities are stored. Prepare for summing over sectors
F(L0VAV+I)=WRF*VSECT*SECTP ! SECTP is the sector probability
IF(LAWS) THEN ! Lawson
DO IC=1,NCRIT
LPRO(IC,I)=LPRO(IC,I)*SECTP
if(lbprb(ic)>0) f(l0prb(ic)+i)=LPRO(IC,I)
ENDDO
ELSE ! Probability of exceeding or NEN
F(L0PRO+I)=F(L0PRO+I)*SECTP
IF(NEN) F(L0DAN+I)=F(L0DAN+I)*SECTP
ENDIF
ENDIF ! end of SECTP block
ENDIF   ! end of .NOT.SLD block
ENDDO     ! end of loop over slab
ENDIF       ! end of probability table block
IF(NEN) THEN ! assign NEN values
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
ELSEIF(LAWS) THEN ! assign Lawson Criterion valus
DO I=1,NXNY
F(L0LAWS+I)=NCRIT
IF(LPRO(IC,I)>=LTHRESH(IC)) F(L0LAWS+I)=IC
ENDDO
ENDDO
ENDIF ! end of NEN or Lawson
ENDIF   ! end of SHOWCOMF block
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(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
```