c<080223
C file-name GXBLIN.HTM 200623
c-->
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C-------------------------------------------------------------------------
SUBROUTINE GXBLIN
C-------------------------------------------------------------------------
C
C PHOENICS V2020
C Author: M.R.Malin / J C Ludwig
C Date: 24/09/04 - 16/04/20
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 = {QTAU/Kappa}*ln((z-d)/zo) or Q = Qref*(z/zref)^a
C
C k = QTAU^2/sqrt(cmucd)
C
C e = QTAU^3/(K*[z-d]) or f = QTAU/(Kappa*sqrt(cmucd)*(z-d))
C
C where the friction velocity QTAU is given by:
C
C Q*= Qref*K/log([zref-d]/zo)
C
C and
C
C Kappa = 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 Pasquill Stability-Class Inlet Profiles
C
C Class-F profiles
C Q = (QTAU/Kappa)*(ln((z-d)/zo) - PSIU)
C T = To - (g/Cp)*(Z-ZTo) + (T*/kappa)*(ln((z-d)/zo) - PSIT)
C k = {QTAU^2/sqrt(cmucd)} * SQRT(1. - ZETA/PHI)
C e = {QTAU^3/(K*[z-d])} *PHI*(1.-ZETA/PHI) or
C f = {QTAU/(Kappa*sqrt(cmucd)*(z-d))} * PHI
C
C where PSIU = -5*z/LMO Monin-Obukhov (MO) Similarity parameter
C PSIT = -5*z/LMO MO Similarity parameter for heat
C PHI = 1. + 5*z/LMO MO Similarity parameter for turbulence
C LMO = MO length scale
C T* = -Qw/(rho*Cp*QTAU) Friction temperature
C
C This facility requires the following SPEDAT statements to be set in the Q1
C input file:
C SPEDAT(SET,BLIN,MOLEN,I,MOLEN) : MOLEN = 0 User-specified Monin-Obukhov length
C = 1 TNO formulae
C = 2 PHAST formulae
C SPEDAT(SET,BLIN,ITPRO,I,ITPRO) : ITPRO = 0 No Pasquill profile, so neutral
C = 1 Class A: Very unstable (LMO<0)
C = 2 Class B: Moderately unstable (LMO<0)
C = 3 Class C: Slightly unstable (LMO<0)
C = 4 Class D: Neutral (LMO=0)
C = 5 Class E: Slightly stable (LMO>0)
C = 6 Class F: Moderately stable (LMO>0)
C
C Nb: ITPRO=4 is neutral with a uniform temperature profile, although this may be used
C in the future to code a neutral logarithmic temperature profile.
C
C SPEDAT(SET,BLIN1,GT0 ,R,To) : T0
C SPEDAT(SET,BLIN1,GTZ0 ,R,ZTo) : Z height for T0
C SPEDAT(SET,BLIN1,QWALL,R,Qw) : heat flux (computed in pre-processor for terrain heat flux)
C
C Limitations:
C No provision has been made as yet for scalar profiles, and the logarithmic
C profile is restricted to the fully-rough wall law, as encountered in the
C atmospheric boundary layer. At present, the GXBLIN facility cannot be used
C for BFC=T, GCV=T or CCM=T.
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/TSKEMI/ LBKP,LBKT,LBET,LBVOSQ,LBOMEG
COMMON/LRNTM4/RSCMCD,CDDAK2,TAUDKE,RTTDKE,AKC,EWC
COMMON/KWMOD3/CWWALL,SIGK1,SIGK2,SIGW1,SIGW2,CWA1,CWA2,CWB1,CWB2,
1 BETAST,CDSIG
COMMON/KWMOD1/LBKWF1(7),LBBF1,LBBF2,LBBF3,LBSIGK,LBSIGW,LBKWF2(5)
COMMON/RHILO/HI3D,RLO3D
COMMON/IHILO/IXHI3D,IYHI3D,IZHI3D,IXLO3D,IYLO3D,IZLO3D,IMAXC(150)
COMMON/PASQLF/ITPRO,PASQBUOY,BUOSSG,MOLEN
LOGICAL PASQBUOY,BUOSSG
LOGICAL MASKPT,SLD,QEQ,QGT,QLE,QNE, LPTDMN,LSOLID,EQZ,
1 LWP,SWP,WWP,LTZ,SHOWCOMF,WEICOEF,NEN,LPAR,LAWS,NEZ,
1 LGWC
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,DLIM*4
CHARACTER PASQUILL*16
PARAMETER (MAXDMN = 10)
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),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(:)
INTEGER, ALLOCATABLE :: L0HI(:,:)
REAL, ALLOCATABLE :: WNDA(:,:),VEL_TAB(:,:), KE_TAB(:,:)
REAL NORML(3)
REAL, ALLOCATABLE :: LTHRESH(:), LUTHR(:),LPRO(:,:)
REAL THRESHD(7),UTHRD(7)
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,LBWAT,LBRHIN,LBTREF,L0TREF0,
1 L0PIN0,L0RHIN0,LBPIN
SAVE AXDIR, RHOIN,WAMPVR, WMAST,WNDA,UTHR,DTHR,SECTP,VSECT,
1 AW,AKW, WDIRN, ANG, GT0,GZT0,GALR,GQWALL,GLMO
SAVE LTHRESH,LUTHR,LPRO
SAVE SHOWCOMF,WEICOEF,NEN,LAWS
SAVE TYPECOMF
SAVE VEL_TAB, KE_TAB, DLIM, VDIR
DATA THRESHD /0.06, 0.02, 0.04, 0.06, 0.01,0.0, 0.0 /
DATA UTHRD /10.95, 10.95, 8.25, 5.6, 5.6, 0.0, 0.0 /
c***********************************************************************
c
IF(USP) RETURN
IXL=IABS(IXL)
NAMSUB='GXBLIN'; NAMFUN=' '
C*****************************************************************
C--- GROUP 1. Preliminaries
IF(IGR==1) THEN
DLIM(1:1)=' '; DLIM(2:2)=','; DLIM(3:3)=';'; DLIM(4:4)=CHAR(9)
C
C... Group 1 Section 1
C =================
IF(ISC==1) THEN
IF(.NOT.NULLPR) THEN
CALL WRYT40('GXBLIN of: 200623')
CALL WRIT40('GXBLIN of: 200623')
ENDIF
CALL GXMAKE(L0H,NXNY,'HDIS') ! storage for grid node heights
C
C... Group 1 Section 2
C =================
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
C
IF(NBLIN(IDMN)==0) RETURN ! no BLINs so nothing to do
IF(.NOT.ALLOCATED(L0HI)) THEN
ALLOCATE(L0HI(MAXDMN,NBLIN(IDMN)),STAT=IERR)
IF(IERR==0) THEN
L0HI=0
ELSE
CALL WRITBL
CALL WRITST
CALL WRIT40('Memory allocation error in GXBLIN')
CALL WRITST
CALL SET_ERR(557,'Memory allocation error in GXBLIN',1)
RETURN
ENDIF
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 SPLTZZZ(CVAL,WD,NWDS,NCHARS,LL,DLIM,20)
WDIR(I)=RRDZZZ(1);WSPD(I)=RRDZZZ(2);AIRTEMP(I)=RRDZZZ(3)
CALL GETSDC('BLINA',LINE(1:LL),CVAL)
LL=LENGZZ(CVAL)
CALL SPLTZZZ(CVAL,WD,NWDS,NCHARS,LL,DLIM,20)
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)
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
IF(ITEM1>0) THEN
ITPRO=0 ! Uniform temperature profile
CALL GETSDI('BLIN','ITPRO',ITPRO)
C
IF(ITPRO.GT.0) THEN
GZT0 = 0.0 ! ZT0 - surface-temperature reference height
CALL GETSDR('BLIN','GZT0',GZT0)
GT0 = TAIR ! Temperature at ZT0
CALL GETSDR('BLIN','GT0',GT0)
GQWALL = 0.0 ! Ground heat flux
CALL GETSDR('BLIN','GQWALL',GQWALL)
GALR = 9.81/CP1 ! Dry adiabatic lapse rate
MOLEN=2 ! User-specified Monin-Obukhov length
CALL GETSDI('BLIN','MOLEN',MOLEN)
C
IF(ITPRO==1) THEN ! Class A: Very unstable
PASQUILL='Pasquill Class A'
IF(MOLEN==1) THEN
GLS=33.162 ; GZS=1117.0 ! TNO
ELSE IF(MOLEN==2) THEN
ALMO=-11.4 ; BLMO=0.1 ! PHAST
ENDIF
ELSEIF(ITPRO==2) THEN ! Class B: Moderately unstable
PASQUILL='Pasquill Class B'
IF(MOLEN==1) THEN
GLS=32.258 ; GZS=11.46 ! TNO
ELSE IF(MOLEN==2) THEN
ALMO=-26.0 ; BLMO=0.17 ! PHAST
ENDIF
ELSEIF(ITPRO==3) THEN ! Class C: Slightly unstable
PASQUILL='Pasquill Class C'
IF(MOLEN==1) THEN
GLS=51.787 ; GZS=1.324 ! TNO
ELSE IF(MOLEN==2) THEN
ALMO=-123.0 ; BLMO=0.3 ! PHAST
ENDIF
ELSEIF(ITPRO==4) THEN ! Class D: Neutral
PASQUILL='Pasquill Class D'
ELSEIF(ITPRO==5) THEN ! Class E: Slightly stable
PASQUILL='Pasquill Class E'
IF(MOLEN==1) THEN
GLS=-48.33 ; GZS=1.262 ! TNO
ELSE IF(MOLEN==2) THEN
ALMO=123.0 ; BLMO=0.3 ! PHAST
ENDIF
ELSEIF(ITPRO==6) THEN ! Class F: Moderately stable
PASQUILL='Pasquill Class F'
C.. Monin Obukhov length
IF(MOLEN==1) THEN
GLS=-31.325 ; GZS=19.36 ! TNO
ELSE IF(MOLEN==2) THEN
ALMO=26.0 ; BLMO=0.17 ! PHAST
ENDIF
ENDIF
IF(MOLEN==0) THEN
GLM0 = 1.0
CALL GETSDR('GLMO','GLMO',GLMO) ! User value
ELSE IF(MOLEN==1) THEN
ZODZS= AMAX1(ZO,0.5)/GZS
GLMO = GLS/LOG10(ZODZS) ! TNO
ELSE IF(MOLEN==2) THEN
GLMO=ALMO*(ZO)**BLMO ! PHAST
ENDIF
PASQBUOY=.FALSE.; BUOSSG=.FALSE.
DO IP=1,NUMPAT
IF(NAMPAT(IP)(1:8)=='BUOYANCY') THEN
PASQBUOY=.TRUE.
DO JIND=U1,W1,2
CALL GETCV(IP,JIND,GCO,GVAL)
IF(QEQ(GVAL,GRND3)) BUOSSG=.TRUE.
ENDDO
ENDIF
ENDDO
LBTREF=LBNAME('TREF')
IF(LBTREF<=0) CALL GXMAKE0(L0TREF0,NX*NY*NZ,'TREF')
IF(.NOT.BUOSSG) THEN
LBPIN =LBNAME('PIN')
IF(LBPIN<=0) CALL GXMAKE0(L0PIN0,NX*NY*NZ,'PIN')
LBRHIN =LBNAME('RHIN')
IF(LBRHIN<=0) CALL GXMAKE0(L0RHIN0,NX*NY*NZ,'RHIN')
ENDIF
CALL WRITST
CALL WRIT40(' Pasquill-Class data ')
CALL WRIT40(' --------------------')
WRITE(14,*) ' ITPRO = ',ITPRO
WRITE(14,*) ' ',PASQUILL
WRITE(14,*) ' Ground surface temperature = ',GT0, ' C'
WRITE(14,*) ' Monin-Obukhov length = ',GLMO,' m'
WRITE(14,*) ' Ground heat flux = ',GQWALL,
1 ' W/m^2'
WRITE(14,*) ' Dry adiabatic lapse rate = ',GALR*1.E3,
1 ' C/km'
IF(BUOSSG) THEN
WRITE(14,*) ' Boussinesq buoyancy'
ELSE
WRITE(14,*) ' Density-difference buoyancy'
ENDIF
CALL WRITBL
ENDIF
ENDIF
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)
C... Wind Attenuation Coefficient WAT
LBWAT=LBNAME('WAT') ! WAT = (Vabs/(profile speed at local height))-1
IF((LBWAF>0.OR.LBWAT>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)
NCRIT=MIN(NCRIT,7)
ALLOCATE(LTHRESH(NCRIT),LUTHR(NCRIT),LPRO(NCRIT,NX*NY),
1 STAT=IERR)
LTHRESH=0.0; LUTHR=0.0; LPRO=0.0 ! initialise
DO IC=1,NCRIT
LTHRESH(IC)=THRESHD(IC)*100.; LUTHR(IC)=UTHRD(IC) ! take default from data statement
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
LBVAV=LBNAME('VAV')
IF(WEICOEF) THEN ! read Weibull parameters
LGWC=.FALSE.;CALL GETSDL('FLAIR','WEIB-GWC',LGWC)
NSECT=12; REWIND LU
READ(LU,'(A)') LINE; SITENAME=LINE ! site name
CALL GETSDR('FLAIR','WEIB-H',DHEIGHT)
CALL GETSDR('FLAIR','WEIB-AW',AW)
CALL GETSDR('FLAIR','WEIB-AKW',AKW)
CALL GETSDR('FLAIR','WEIB-SECTP',SECTP)
ELSE ! Read probability bins
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 SPLTZZZ(LINE,WD,NWDS,NCHARS,LL,DLIM,20)
DHEIGHT=RRDZZZ(3)
READ(LU,'(A)') LINE; LL=LENGZZ(LINE) ! no of sectors
CALL SPLTZZZ(LINE,WD,NWDS,NCHARS,LL,DLIM,20)
NSECT=IRDZZZ(1)
ALLOCATE(WNDA(NBINS+1,NSECT+1),STAT=IERR)
WNDA=0.0 ! initialise bins
READ(LU,'(A)') LINE; LL=LENGZZ(LINE)
CALL SPLTZZZ(LINE,WD,NWDS,NCHARS,LL,DLIM,20)
DO ISEC=1,NSECT
WNDA(1,ISEC+1)=RRDZZZ(ISEC)
ENDDO
PSUM=TINY
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 ! read the probabilities for each bin
READ(LU,'(A)') LINE
LL=LENGZZ(LINE)
CALL SPLTZZZ(LINE,WD,NWDS,NCHARS,LL,DLIM,20)
DO ISEC=1,NSECT+1
WNDA(IB,ISEC)=RRDZZZ(ISEC)
ENDDO
ENDDO
DO ISEC=2,NSECT+1 ! normalise probabilities for each direction
PSUM=TINY
DO IB=2,NBINS ! sum probabilities over all bins
PSUM=PSUM+WNDA(IB,ISEC)
ENDDO
DO IB=2,NBINS ! 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 average velocity stored
IF(WEICOEF) THEN ! Weibull coefficients
VSECT=AW*GAMMA(1.0+1.0/AKW) ! Mean value of Weibull distribution
!!! VSECT=AW*LOG(2.0)**(1.0/AKW) ! Median value of Weibull distribution
ELSE ! Probability table
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
ENDIF
QREF=VSECT; REFH=DHEIGHT
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
DO II=1,NBLIN(IDMN)
F(L0QREF(IDMN)+II)=QREF
F(L0HREF(IDMN)+II)=REFH
C... Update values for BLIN patches for use as inlet values in Group 13
IUP=F(L0VDIR(IDMN)+II)
IF(IUP==3) THEN ! Up Z
VELX=-VSECT*SIN(ANGR)
VELY=-VSECT*COS(ANGR)
VELZ=0.0
ELSEIF(IUP==2) THEN ! Up Y
VELX=-VSECT*COS(ANGR)
VELY=0.0
VELZ=-VSECT*SIN(ANGR)
ELSE ! Up X
VELX=0.0
VELY=-VSECT*SIN(ANGR)
VELZ=-VSECT*COS(ANGR)
ENDIF
C... Set inlet velocity components
F(L0VELX(IDMN)+II)=VELX
F(L0VELY(IDMN)+II)=VELY
F(L0VELZ(IDMN)+II)=VELZ
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 SPLTZZZ(LINE,WD,NWDS,NCHARS,LL,DLIM,20) ! 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 SPLTZZZ(LINE,WD,NWDS,NCHARS,LL,DLIM,20) ! 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'
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)
C
C... Group 1 Section 3
C =================
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
C
C The followed yet to be coded.
C IF(ITPRO > 0 .AND. ITPRO /= 4) THEN
C GPHIM = 1.0-PSIFT(GZMDH,GLMO)
C ENUT=ENUT/GPHIM
C ENDIF
C
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
IF((INDVAR.EQ.KE.OR.INDVAR.EQ.LBOMEG).AND.
1 (IENUTA.GE.17.AND.IENUTA.LE.20)) THEN
L0VIST=L0F(VIST); L0BF1 =L0F(LBBF1)
IF(INDVAR.EQ.KE) THEN
GSIG1=SIGK1 ; GSIG2=SIGK2
ELSE ! w
GSIG1=SIGW1 ; GSIG2=SIGW2
ENDIF
DO IX=IXF,IXL
IADD=NY*(IX-1)
DO IY=IYF,IYL
I=IY+IADD
GBF1 = F(L0BF1+I)
GPRTRB = (GBF1*GSIG1 + (1.-GBF1)*GSIG2)
F(L0CO+I)=F(L0DEN+I)*F(L0VIST+I)/(GPRTRB*GDH)
ENDDO
ENDDO
ELSE
DO IX=IXF,IXL
IADD=NY*(IX-1)
DO IY=IYF,IYL
I=IY+IADD
F(L0CO+I)=F(L0DEN+I)*F(L0VIST+I)/(PRT(INDVAR)*GDH)
ENDDO
ENDDO
ENDIF
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
C ====================================
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
C ---------------
C.. VELCON = RHOIN*(QTAU/KAPPA)
IF(ITPRO > 0 .AND. ITPRO /= 4) THEN
VELCON=RHOIN*VELIN/(LOG(REFH/ZO)-PSIF(REFH,GLMO))
ELSE ! Neutral
VELCON=RHOIN*VELIN/LOG(REFH/ZO)
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) 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
C.. VELCON = RHOIN*(QTAU/KAPPA)
IF(ITPRO > 0 .AND. ITPRO /= 4) THEN
F(L0VAL+I)=VELCON*(LOG(GHDZO)-PSIF(GH,GLMO))*AMULT
ELSE ! Neutral
F(L0VAL+I)=VELCON*LOG(GHDZO)*AMULT
ENDIF
ELSE
F(L0VAL+I)=0.
ENDIF
ENDDO
ENDDO
ENDIF
C
C... Inlet velocity components
C =========================
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
C ---------------
ELSEIF(IBLTY==1) THEN
C.. VELCON = QTAU/KAPPA
IF(ITPRO > 0 .AND. ITPRO /= 4) THEN
VELCON=VELIN/(LOG((REFH-WALLB)/ZO)-PSIF(REFH,GLMO))
ELSE ! Neutral
VELCON=VELIN/LOG((REFH-WALLB)/ZO)
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)
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.. VELCON = QTAU/KAPPA
IF(ITPRO > 0 .AND. ITPRO /= 4) THEN
F(L0VAL+I)=VELCON*(LOG(GHDZO)-PSIF(GH,GLMO))
ELSE ! Neutral
F(L0VAL+I)=VELCON*LOG(GHDZO)
ENDIF
ELSE
F(L0VAL+I)=0.0
ENDIF
if(l0velin>0) f(l0velin+i)=f(l0val+i)
ENDDO
ENDDO
ENDIF
C
C... inlet k and ep/omeg values
C ==========================
C
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
REFHMD = REFH-WALLB
IF(ITPRO > 0 .AND. ITPRO /= 4) THEN
QTAU = ABS(AK*QREF/(LOG((REFHMD)/ZO)-PSIF(REFHMD,GLMO)))
ELSE
QTAU = ABS(AK*QREF/(LOG((REFHMD)/ZO)))
ENDIF
QTAU2 = QTAU*QTAU
C... Taudke=cmucd**0.5
GKEIN = QTAU2/TAUDKE ! Nb: Taudke=sqrt(CmuCd)
GEPCON = QTAU2*QTAU/AK
IF(INDVAR==KE) THEN
CALL FN1(VAL,GKEIN)
IF(ITPRO > 0 .AND. ITPRO /= 4) THEN
lbkein=lbname('KEIN') ; if(lbkein>0) l0kein=l0f(lbkein)
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
GPHIE = PHIE(GZMDH,GLMO)
GPHIM = 1.0-PSIFT(GZMDH,GLMO)
F(L0VAL+I)=GKEIN*(GPHIE/GPHIM)**0.5
ELSE
F(L0VAL+I)=0.0
ENDIF
if(lbkein>0) f(l0kein+i)=f(l0val+i)
ENDDO
ENDDO
ENDIF
ENDIF
IF(INDVAR==EP.OR.INDVAR.EQ.LBOMEG) THEN
IF(INDVAR==LBOMEG) THEN
lbomin=lbname('OMIN') ; if(lbomin>0) l0omin=l0f(lbomin)
L0OMEG=L0F(LBOMEG)
GOMCON=QTAU/(AK*TAUDKE)
ELSE
lbepin=lbname('EPIN') ; if(lbepin>0) l0epin=l0f(lbepin)
L0EP=L0F(EP)
ENDIF
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
IF(INDVAR==EP) THEN
GEPIN=GEPCON/GZMDH ! e = ustar^3/(ak*(z-d))
IF(ITPRO > 0 .AND. ITPRO /= 4) THEN
GEPIN=GEPIN*PHIE(GZMDH,GLMO)
ENDIF
F(L0VAL+I)=GEPIN
ELSE IF(INDVAR==LBOMEG) THEN
GOMEGI=GOMCON/GZMDH ! f = ustar/(sqrt(CmuCD)*ak*(z-d))
IF(ITPRO > 0 .AND. ITPRO /= 4) THEN
GPHIE = PHIE(GZMDH,GLMO)
GPHIM = 1.0-PSIFT(GZMDH,GLMO)
F(L0VAL+I)=GOMEGI*(GPHIE*GPHIM)**0.5
ELSE
F(L0VAL+I)=GOMEGI
ENDIF
ENDIF
ELSE
F(L0VAL+I)=0.0
ENDIF
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
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
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
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
C... External air temperature
C ========================
ELSEIF(INDVAR==ITEM1) THEN
IF(ITPRO==0.OR.ITPRO==4) THEN ! Uniform profile or Neutral
CALL FN1(VAL,TAIR)
ELSE IF(ITPRO > 0 ) THEN
C... Pasquill temperature profiles
lbtin=lbname('TIN')
if(lbtin>0) l0tin=l0f(lbtin)
IF(ITPRO > 0 .AND. ITPRO /= 4) THEN
C.. friction velocity
REFHMD = REFH-WALLB ! (z-d)
QTAU = AK*QREF/(LOG((REFHMD)/ZO)-PSIF(REFHMD,GLMO))
C... Friction temperature
GTSTAR= -GQWALL/(RHOIN*CP1*QTAU) ; GTSDAK=GTSTAR/AK
ENDIF
L0VAL=L0F(VAL)
IPBC=0
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
C IF(ITPRO==) THEN
C.. linear profile:- T = T0 +gam*(z-z0)
C GTLIN = GT0 - GALR*(GH-GZT0)
C F(L0VAL+I)= GTLIN
IF(ITPRO > 0 .AND. ITPRO /= 4) THEN ! Pasquill - Stable (LMO > 0)
C.. log profile:- T = T0-gam*(z-z0)+(T*/k)*log(z/z0)
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
GTLIN = GT0 - GALR*(GH-GZT0)
GPSIT = PSIFT(GH-WALLB,GLMO)
F(L0VAL+I)=GTLIN+GTSDAK*(LOG(GHDZO)-GPSIT)
ENDIF
ELSE
F(L0VAL+I)=0.0
ENDIF
if(lbtin>0) f(l0tin+i)=f(l0val+i)
ENDDO
ENDDO
ENDIF
C
C... External air humidity
C =====================
ELSEIF(INDVAR==IHUM) THEN
CALL FN1(VAL,HUMIN)
ENDIF
C
C... Group 13, Section 14 Density-difference buoyancy forces
C
C... Buoyancy Forces
C ===============
ELSEIF(IGR==13.AND.ISC==14) THEN ! Density difference
C.. Entry here is from GXBUOY with VAL=-g ; Rhoref = f(height)
C for all velocities, so need to identify active gravity direction
IF(INDVAR==3) THEN
L0VAL= L0F(VAL) ; L0D9 = L0F(LD9)
C DO I=1,NXNY
C F(L0VAL+I) = F(L0VAL+I)*(F(L0D9+I) - BUOYD)/F(L0D9+I)
C ENDDO
L0DEN=L0F(DEN1)
IF(LBRHIN>0) THEN
L0RHIN=L0F(LBRHIN)
ELSE
L0RHIN=(IZ-1)*NXNY+L0RHIN0
ENDIF
GGRAVY=BUOYA
DO JX=1,NX-1
JXAD=(JX-1)*NY
DO JY=1,NY
J=JXAD+JY
IF(SLD(J)) THEN
F(L0VAL+J)=0.0
ELSE
GRHO=F(L0DEN+J)
GRHON=F(L0DEN+J+NY)
GRHOAV=0.5*(GRHO+GRHON)
GRHRAV=0.5*(F(L0RHIN+J)+F(L0RHIN+J+NY))
C GRHOREF = BUOYD
F(L0VAL+J)=(1.-GRHRAV/GRHOAV)*GGRAVY
ENDIF
ENDDO
ENDDO
ELSEIF(INDVAR==5) THEN
L0VAL=L0F(VAL)
L0DEN=L0F(DEN1)
IF(LBRHIN>0) THEN
L0RHIN=L0F(LBRHIN)
ELSE
L0RHIN=(IZ-1)*NXNY+L0RHIN0
ENDIF
GGRAVY=BUOYB
DO JX=1,NX
JXAD=(JX-1)*NY
DO JY=1,NY-1
J=JXAD+JY
IF(SLD(J)) THEN
F(L0VAL+J)=0.0
ELSE
GRHO=F(L0DEN+J)
GRHON=F(L0DEN+J+1)
GRHOAV=0.5*(GRHO+GRHON)
GRHRAV=0.5*(F(L0RHIN+J)+F(L0RHIN+J+1))
C GRHOREF = BUOYD
F(L0VAL+J)=(1.-GRHRAV/GRHOAV)*GGRAVY
ENDIF
ENDDO
ENDDO
ELSEIF(INDVAR==7) THEN
L0VAL=L0F(VAL) ; L0D9 = L0F(LD9)
L0DEN=L0F(DEN1); L0DENH=L0F(HIGH(DEN1))
IF(LBRHIN>0) THEN
L0RHIN=L0F(LBRHIN); L0RHINH=L0F(HIGH(LBRHIN))
ELSE
L0RHIN=(IZ-1)*NXNY+L0RHIN0
L0RHINH=L0RHIN+NXNY
ENDIF
GGRAVY=BUOYC
DO JX=1,NX
JXAD=(JX-1)*NY
DO JY=1,NY
J=JXAD+JY
IF(SLD(J)) THEN
F(L0VAL+J)=0.0
ELSE
GRHO=F(L0DEN+J) ; GRHON=F(L0DENH+J)
GRHOAV=0.5*(GRHO+GRHON)
C.. NB: 1st visit rhin,high=0.0
GRHRAV=0.5*(F(L0RHIN+J)+F(L0RHINH+J))
F(L0VAL+J)=(1.-GRHRAV/GRHOAV)*GGRAVY
ENDIF
ENDDO
ENDDO
ENDIF
C-------------------------------------------------------------------------
C... Group 13, Section 15 Boussinesq buoyancy forces
C
ELSEIF(IGR==13.AND.ISC==15) THEN ! Boussinesq buoyancy
L0TEM1=L0F(ITEM1)
IF(LBTREF>0) THEN
L0TREF=L0F(LBTREF)
ELSE
L0TREF=L0TREF0+(IZ-1)*NXNY
ENDIF
L0VAL= L0F(VAL) ;L0D9 = L0F(LD9)
IF(INDVAR==3) THEN ! U1
C DO I=1,NXNY
C F(L0VAL+I) = F(L0VAL+I)*(F(L0D9+I) - BUOYD)/F(L0D9+I)
C ENDDO
GGRAVY=BUOYA
DO JX=1,NX-1
JXAD=(JX-1)*NY
DO JY=1,NY
J=JXAD+JY
IF(SLD(J)) THEN
F(L0VAL+J)=0.0
ELSE
GTEM=F(L0TEM1+J)
GTEMN=F(L0TEM1+J+NY)
GTEMAV=0.5*(GTEM+GTEMN)
GTRFAV=0.5*(F(L0TREF+J)+F(L0TREF+J+NY))
F(L0VAL+J)=(GTEMAV-GTRFAV)*GGRAVY*DVO1DT
ENDIF
ENDDO
ENDDO
ELSEIF(INDVAR==5) THEN ! V1
GGRAVY=BUOYB
DO JX=1,NX
JXAD=(JX-1)*NY
DO JY=1,NY-1
J=JXAD+JY
IF(SLD(J)) THEN
F(L0VAL+J)=0.0
ELSE
GTEM=F(L0TEM1+J)
GTEMN=F(L0TEM1+J+1)
GTEMAV=0.5*(GTEM+GTEMN)
GTRFAV=0.5*(F(L0TREF+J)+F(L0TREF+J+1))
C GRHOREF = BUOYD
F(L0VAL+J)=(GTEMAV-GTRFAV)*GGRAVY*DVO1DT
ENDIF
ENDDO
ENDDO
ELSEIF(INDVAR==7) THEN ! W1
L0TEMH=L0F(HIGH(ITEM1))
IF(LBTREF>0) THEN
L0TREFH=L0F(HIGH(LBTREF))
ELSE
L0TREFH=L0TREF+NXNY
ENDIF
GGRAVY=BUOYC
DO JX=1,NX
JXAD=(JX-1)*NY
DO JY=1,NY
J=JXAD+JY
IF(SLD(J)) THEN
F(L0VAL+J)=0.0
ELSE
GTEM=F(L0TEM1+J) ; GTEMN=F(L0TEMH+J)
GTEMAV=0.5*(GTEM+GTEMN)
C.. NB: 1st visit rhin,high=0.0
GTRFAV=0.5*(F(L0TREF+J)+F(L0TREFH+J))
F(L0VAL+J)=(GTEMAV-GTRFAV)*GGRAVY*DVO1DT
ENDIF
ENDDO
ENDDO
ENDIF
C-----------------------------------------------------------------------------
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==4) THEN
C============================================================
C... Group 19, Section 4. End of IZ slab
C Compute Pasquill-F Buoyancy-reference profiles for use
C in Group 13 buoyancy source terms
C Entry here from Grex3 only if PASQBUOY=T
C============================================================
IBLIN=1
C.. inlet velocity magnitude at reference height
QREF=F(L0QREF(IDMN)+IBLIN)
C.. inlet temperature
TAIR=F(L0TAIR(IDMN)+IBLIN)
C.. vertical coordinate direction
IVDIR=NINT(F(L0VDIR(IDMN)+IBLIN))
C.. inlet density
RHOIN=F(L0DENS(IDMN)+IBLIN)
C.. effective roughness length
ZO=F(L0ZREF(IDMN)+IBLIN)
C.. reference height for wind reference velocity
REFH=F(L0HREF(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.. Compute inlet (reference) temperature profile
C =============================================
IF(LBTREF>0) THEN
L0TREF=L0F(LBTREF)
ELSE
L0TREF=L0TREF0+(IZ-1)*NXNY
ENDIF
C... Reference Temperature
c... Datum is at y=0 ; p=p0 rho=rhoo T=To
GY0 = 0.0 ; GP0 = 0.0
GRGAS= 286.7
GRHO0= (PRESS0+GP0)/(GRGAS*(GT0+TEMP0))
C GZG = GH
REFHMD = REFH-WALLB ! (z-d)
QTAU = AK*QREF/(LOG((REFHMD)/ZO)-PSIF(REFHMD,GLMO))
GTSTAR= -GQWALL/(RHOIN*CP1*QTAU) ; GTSDAK=GTSTAR/AK
IPBC=0
DO IX=1,NX
IADD=NY*(IX-1)
DO IY=1,NY
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
GTLIN = GT0 - GALR*(GH-GZT0)
GPSIT = PSIFT(GH-WALLB,GLMO)
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
ENDIF
F(L0TREF+I)=GTLIN+GTSDAK*(LOG(GHDZO)-GPSIT)
ENDDO
ENDDO
C
IF(BUOSSG) RETURN
C
C Density-difference buoyancy profiles
C
C.. Compute inlet (reference) pressure & density profiles
C =====================================================
IF(LBPIN>0) THEN
L0PIN = L0F(LBPIN)
L0PINL = L0F(LOW(LBPIN))
ELSE
L0PIN = (IZ-1)*NXNY+L0PIN0
L0PINL = L0PIN-NXNY
ENDIF
IF(LBRHIN>0) THEN
L0RHIN = L0F(LBRHIN)
L0RHINL = L0F(LOW(LBRHIN))
ELSE
L0RHIN = (IZ-1)*NXNY+L0RHIN0
L0RHINL = L0RHIN-NXNY
ENDIF
IF(IVDIR==1) THEN ! Up X
GGRAVY=BUOYA
L0XG = L0F(XG2D) ; L0DXG=L0F(DXG2D)
C... Pressure & density at first vertical cell
JX=1
DO JY=1,NY
J=(JX-1)*NY+JY
GDYP=F(L0XG+1)
GTERM1 = 0.5*GRHO0*GGRAVY*GDYP
GTERM2 = 0.5*GGRAVY*GDYP/(GRGAS*(F(L0TREF+J)+TEMP0))
F(L0PIN+J)=(GP0+GTERM1)/(1.-GTERM2)
F(L0RHIN+J)=(PRESS0+F(L0PIN+J))/(GRGAS*(F(L0TREF+J)+TEMP0))
ENDDO
C... Remainder of reference pressure & density profiles
DO JX=2,NX
JXAD=(JX-1)*NY
DO JY=1,NY
J=JXAD+JY
JS=J-NY
GDYGS=F(L0DXG+JS)
GTERM1 = 0.5*F(L0RHIN+JS)*GGRAVY*GDYGS
GTERM2 = 0.5*GGRAVY*GDYGS/(GRGAS*(F(L0TREF+J)+TEMP0))
F(L0PIN+J)=(F(L0PIN+JS)+GTERM1)/(1.-GTERM2)
F(L0RHIN+J)=(PRESS0+F(L0PIN+J))/(GRGAS*(F(L0TREF+J)+TEMP0))
ENDDO
ENDDO
C
ELSEIF(IVDIR==2) THEN ! Up y
C
GGRAVY=BUOYB
L0YG = L0F(YG2D) ; L0DYG=L0F(DYG2D)
C... Pressure & density at first vertical cell
DO JX=1,NX
J=(JX-1)*NY+1
GDYP=F(L0YG+1)
GTERM1 = 0.5*GRHO0*GGRAVY*GDYP
GTERM2 = 0.5*GGRAVY*GDYP/(GRGAS*(F(L0TREF+J)+TEMP0))
F(L0PIN+J)=(GP0+GTERM1)/(1.-GTERM2)
F(L0RHIN+J)=(PRESS0+F(L0PIN+J))/(GRGAS*(F(L0TREF+J)+TEMP0))
ENDDO
C... Remainder of reference pressure & density profiles
DO JX=1,NX
JXAD=(JX-1)*NY
DO JY=2,NY
J=JXAD+JY
JS=J-1
GDYGS=F(L0DYG+JS)
GTERM1 = 0.5*F(L0RHIN+JS)*GGRAVY*GDYGS
GTERM2 = 0.5*GGRAVY*GDYGS/(GRGAS*(F(L0TREF+J)+TEMP0))
F(L0PIN+J)=(F(L0PIN+JS)+GTERM1)/(1.-GTERM2)
F(L0RHIN+J)=(PRESS0+F(L0PIN+J))/(GRGAS*(F(L0TREF+J)+TEMP0))
ENDDO
ENDDO
C
ELSEIF(IVDIR==3) THEN ! Up z
C
GGRAVY=BUOYC
L0ZG = L0F(ZGNZ) ; L0DZG=L0F(DZGNZ)
C... Pressure & density at first vertical cell
IF(IZ.EQ.1) THEN
DO J=1,NXNY
GDZP=F(L0ZG+IZ)
GTREF=F(L0TREF+J)
GTERM1 = 0.5*GRHO0*GGRAVY*GDZP
GTERM2 = 0.5*GGRAVY*GDZP/(GRGAS*(F(L0TREF+J)+TEMP0))
F(L0PIN+J)=(GP0+GTERM1)/(1.-GTERM2)
F(L0RHIN+J)=(PRESS0+F(L0PIN+J))/(GRGAS*(GTREF+TEMP0))
ENDDO
C... Remainder of reference pressure & density profiles
ELSE
DO J=1,NXNY
GDZGS=F(L0DZG+IZ-1)
GTERM1 = 0.5*F(L0RHINL+J)*GGRAVY*GDZGS
GTERM2 = 0.5*GGRAVY*GDZGS/(GRGAS*(F(L0TREF+J)+TEMP0))
F(L0PIN+J)=(F(L0PINL+J)+GTERM1)/(1.-GTERM2)
F(L0RHIN+J)=(PRESS0+F(L0PIN+J))/(GRGAS*(F(L0TREF+J)+TEMP0))
ENDDO
ENDIF
C
ENDIF
C-------------------------------------------------------------------------
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.LBWAT>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)
ENDIF
IF(LBWAT>0) THEN
L0WAT=L0F(LBWAT); L0ZG=L0F(ZGNZ)
ENDIF
IF(LBWAF>0.OR.LBWAT>0) THEN
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
lbghh=lbname('GHH')
if(lbghh>0) l0ghh=l0f(lbghh)
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.OR.LBWAT>0) THEN
GH=F(L0ZG+IZ)-F(L0HIWAF(IDMN)+I)
if(lbghh>0) f(l0ghh+i)=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
IF(LBWAF>0) F(L0WAF+I)=VABS/(WAFVR+TINY)
IF(LBWAT>0) F(L0WAT+I)=(VABS/(WAFVR+TINY))-1.0
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
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))
if(lbprb(ic)>0) f(l0prb(ic)+i)=LPRO(IC,I)
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
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
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
INEN=1
ELSEIF(PRO>=0.025.AND.PRO<0.05) THEN
INEN=2
ELSEIF(PRO>=0.05.AND.PRO<0.1) THEN
INEN=3
ELSEIF(PRO>=0.1.AND.PRO<0.2) THEN
INEN=4
ELSEIF(PRO>=0.2) THEN
INEN=5
ENDIF
IF(PDAN>=0.05.AND.PDAN<=0.3) THEN
INEN=6
ELSEIF(PDAN>0.3) THEN
INEN=7
ENDIF
F(L0NEN+I)=INEN
ENDDO
ELSEIF(LAWS) THEN ! assign Lawson Criterion valus
DO I=1,NXNY
F(L0LAWS+I)=1
DO IC=NCRIT,1,-1 ! start with most comfortable and work up
IF(LPRO(IC,I)>=LTHRESH(IC)) F(L0LAWS+I)=NCRIT-IC+1
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(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
C--------------------------------------
C Monin-Obukhov similarity-parameter:- psif (momentum)
FUNCTION PSIF(GH,GLMO)
ZETA=GH/GLMO
IF(GLMO.GT.0.0) THEN ! Stable
PSIF = -5.*ZETA
ELSE ! Unstable
PI = 3.141592654
X = (1.0-16.*ZETA)**0.25
T1 = 2.*LOG(0.5*(1.0+X))
T2 = LOG(0.5*(1.0+X*X))
PSIF = T1 + T2 - 2.0*ATAN(X)+0.5*PI
ENDIF
END
C Monin-Obukhov similarity-parameter:- psi fif= 1+5*z/L
FUNCTION FIF (GH,GLMO)
ZETA=GH/GLMO
FIF=1+5.*GH/GLMO ! Stable only
END
C Monin-Obukhov similarity-parameter:- psift (energy)
FUNCTION PSIFT(GH,GLMO)
ZETA=GH/GLMO
IF(GLMO.GT.0.0) THEN ! Stable
PSIFT = -5*ZETA
ELSE ! Unstable
X=(1.0-16.*ZETA)**0.25
PSIFT=2.*LOG(0.5*(1.0+X*X))
ENDIF
END
C Monin-Obukhov similarity-parameter:- phie (turbulence)
FUNCTION PHIE(GH,GLMO)
ZETA=GH/GLMO
IF(GLMO.GT.0.0) THEN ! Stable
PHIE = 1.0 + 4.0*ZETA
ELSE ! Unstable
PHIE = 1.0 - ZETA
ENDIF
END