c
C file name ctgrd.ftn................................. 220306
c
C This file contains two subroutines used for TACT evaporative
C cooling model. First is TACTGR which calculates sources during
C heat and mass exchange; second is modified GXBFC, now called
C PROFIL used to prescribe profiled inlet velocity.
C***************************************************************
SUBROUTINE TACTGR
INCLUDE 'patcmn'
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
PARAMETER (NLG=20, NIG=20, NRG=100, NCG=10)
PARAMETER (NPNAM=200)
C
COMMON/LGRND/LG(NLG)/IGRND/IG(NIG)/RGRND/RG(NRG)/CGRND/CG(NCG)
LOGICAL LG, QNE
CHARACTER*4 CG
CHARACTER*80 BUFF
SAVE
REAL LEWNO, INLETH
EQUIVALENCE (NZM,IG(1)),(NYM,IG(2)),(IHZ1,IG(3)),(IHZ2,IG(4)),
1 (ILY1,IG(5))
EQUIVALENCE (WFLR,RG(1)),(PAMB,RG(2)),(TAMB,RG(3)),(TWAMB,RG(4)),
1 (WAMB,RG(5)),(HSLDA,RG(6)),(HSN,RG(7)),(VSN,RG(8)),
1 (VSLDA,RG(9)),(HPLDA,RG(10)),(HPN,RG(11)),
1 (VPLDA,RG(12)),(VPN,RG(13)),(HRLDA,RG(14)),(HRN,
1 RG(15)),(VRLDA,RG(16)),(VRN,RG(17)),(TINI,RG(18)),
1 (VELIM,RG(19)),(VSTRT,RG(20)),(FILLH,RG(21)),
1 (RAINH,RG(22)),(INLETH,RG(23)),(ELIMH,RG(24)),
1 (TOWERH,RG(25)),(BTMFR,RG(26)),(TOPFR,RG(27)),
1 (ELIMR,RG(28)),(EXTR,RG(29)),(LEWNO,RG(30))
C
DATA GPI,CFF/3.14,4179./
DATA GCO,GC1/3.148856E6,2.372E3/
DATA GSPECV,GSPECG/1.8003E3,1004.832/
DATA GWH2O,GWAIR,GRCONS/18.02,28.97,8314.0/
DATA GTREF,GPREF,HLAT0/273.,610.0,2.5013E6/
C
IXL=IABS(IXL)
IF(IGR.EQ.13) GO TO 13
IF(IGR.EQ.19) GO TO 19
GO TO (1,2,2,2,2,2,2,2,9,2,11,2,13,2,2,2,2,2,19,2,2,
12,2,2),IGR
C*****************************************************************
1 GO TO (1001,1002),ISC
1001 CONTINUE
CALL WRYT40('Call to spec. Ground TACTGR.F of: 010794')
IF(BFC) CALL PROFIL(NX,NY,NZ)
RETURN
1002 CONTINUE
C .. Calculate ambient conditions CAMB, HAMB and DAMB from
C .. wet and dry-bulb temperatures.
GFACT1=GCO*GWH2O*(1./GTREF-1./(TWAMB+GTREF))/GRCONS
GFACT2=GC1*GWH2O*ALOG((TWAMB+GTREF)/GTREF)/GRCONS
GPSAT=GPREF*EXP(GFACT1-GFACT2)
TWBST=GWH2O*GPSAT/(GWAIR*(PAMB-GPSAT)+GWH2O*GPSAT)
XSATWB=TWBST/(1.-TWBST)
HLATWB=GCO-GC1*(TWAMB+GTREF)
C .. Absolute humidity of air XAIR
XAIR=XSATWB-(950.*(TAMB-TWAMB)/HLATWB)
C .. Moisture content of ambient air CAMB
CAMB=XAIR/(1.+XAIR)
WRITE(BUFF,'(A,1PE10.3)') 'CAMB= ',CAMB
CALL PUT_LINE(BUFF,.true.)
GSPECM=GSPECV*CAMB+(1.-CAMB)*GSPECG
C .. HAMB Enthalpy of ambient air
HAMB=TAMB*GSPECM+CAMB*HLAT0
WRITE(BUFF,'(A,1PE10.3)') 'HAMB= ',HAMB
CALL PUT_LINE(BUFF,.true.)
C .. DAMB - density of ambient air
GWMIX=1./(CAMB/GWH2O+(1.-CAMB)/GWAIR)
DAMB=PAMB*GWMIX/(GRCONS*(TAMB+GTREF))
WRITE(BUFF,'(A,1PE10.3)') 'DAMB= ',DAMB
CALL PUT_LINE(BUFF,.true.)
C
C .. Transfer reference density to the buoyancy term and to PROFIL
C .. for calculating wind mass-inflow
BUOYD=DAMB
BFCA=DAMB
SUMEVP=0.0
RETURN
2 CONTINUE
RETURN
C*****************************************************************
C
C--- GROUP 9. Properties of the medium (or media)
C
9 CONTINUE
GO TO (91,92,92,92,92,96,92,92,92,92,92,92,92),ISC
C*****************************************************************
92 CONTINUE
RETURN
91 CONTINUE
C * ------------------- SECTION 1 ---------------------------
C ** DENS=PAMB*GWMIX/(R*GTEM)
CALL SUB2(L0H,L0F(H1),L0C,L0F(C1))
CALL SUB2(L0RH,L0F(DEN1),L0TA,L0F(LBNAME('TAIR')))
DO 910 IX=1,NX
DO 910 IY=1,NY
ICELL=IY+NY*(IX-1)
GH1=HAMB+F(L0H+ICELL)
GWMIX=1./(F(L0C+ICELL)/GWH2O+(1.-F(L0C+ICELL))/GWAIR)
GSPECM=GSPECV*F(L0C+ICELL)+(1.-F(L0C+ICELL))*GSPECG
F(L0TA+ICELL)=(GH1-F(L0C+ICELL)*HLAT0)/GSPECM
F(L0TA+ICELL)=AMAX1(TAMB,AMIN1(F(L0TA+ICELL),TINI))
F(L0RH+ICELL)=PAMB*GWMIX/(GRCONS*(F(L0TA+ICELL)+GTREF))
910 CONTINUE
CALL FN23(DEN1,DAMB)
RETURN
96 CONTINUE
C * ------------------- SECTION 6 ---------------------------
C For ENUL.LE.GRND--- reference laminar kinematic viscosity.
CALL SUB2(L0VIS,L0F(VISL),L0RH,L0F(DEN1))
if(nx.gt.1) then
crj L0ZD=L0B(70)
crj ICL=NY+NY*NX*(IZ-1)
crj ZDIST=F(L0ZD+ICL)
ZDIST=coord1(1,ny,izstep,3)
C .. Check if position is within boundary layer; (rg(40) is
C .. calculated in PROFIL
IF(ZDIST.GT.RG(40)) THEN
SETVS=1.E-5
ELSE
SETVS=RG(36)
ENDIF
C .. RG(36) calculated in PROFIL as rg(36)=0.35*ustar*hei
C .. IG(6) carries cell-value of shell height
IF(IZ.LT.ig(6)) THEN
DO 971 IX=1,NX
DO 971 IY=NYM+1,NY
ICELL=IY+NY*(IX-1)
971 F(L0VIS+ICELL)=SETVS
DO 970 IX=1,NX
DO 970 IY=1,NYM
ICELL=IY+NY*(IX-1)
970 F(L0VIS+ICELL)=RG(35)
ELSE
DO 972 IX=1,NX
DO 972 IY=1,NY
ICELL=IY+NY*(IX-1)
IF(F(L0RH+ICELL).LT.DAMB) THEN
F(L0VIS+ICELL)=RG(35)
ELSE
F(L0VIS+ICELL)=SETVS
ENDIF
972 CONTINUE
ENDIF
ELSE
DO 973 IX=1,NX
DO 973 IY=NYM+1,NY
ICELL=IY+NY*(IX-1)
973 F(L0VIS+ICELL)=1.E-5
DO 974 IX=1,NX
DO 974 IY=1,NYM
ICELL=IY+NY*(IX-1)
974 F(L0VIS+ICELL)=RG(35)
endif
RETURN
C*****************************************************************
C--- GROUP 11.
11 CONTINUE
IF(NPATCH(1:4).EQ.'IBFT') CALL PROFIL(NX,NY,NZ)
RETURN
C*****************************************************************
C
C--- GROUP 13. Boundary conditions and special sources
C
13 CONTINUE
IF(NPATCH(1:3).EQ.'BFT') CALL PROFIL(NX,NY,NZ)
GO TO (130,131,132,133,134,135,136,137,138,139,1310,
11311,1311,1311,1314,1315,1311,1317,1318,1311,1311,1311),ISC
1311 CONTINUE
RETURN
130 CONTINUE
C------------------- SECTION 1 ------------- coefficient = GRND
C ** Resistance coeff. for V1,U1 & w1 in the spray region
C
IF(NPATCH(1:5).EQ.'SPRAY') THEN
CALL SUB2(L0CO,L0F(CO),L0RH,L0F(DEN1))
CALL SUB2(L0W1,L0F(W1),L0WF,L0F(LBNAME('WFLR')))
CALL SUB2(L0V1,L0F(V1),L0AF,L0F(LBNAME('AIRF')))
IF(SOLVE(U1)) L0U1=L0F(U1)
DO 1363 IX=IXF,IXL
DO 1363 IY=IYF,IYL
ICELL=IY+NY*(IX-1)
if(indvar.eq.u1) vel=f(l0u1+icell)
if(indvar.eq.v1) vel=f(l0v1+icell)
if(indvar.eq.w1) vel=f(l0w1+icell)
F(L0CO+ICELL)=(VSLDA*(F(L0WF+ICELL)/F(L0AF+ICELL))+VSN)*
1 ABS(VEL)*F(L0RH+ICELL)*0.5
1363 CONTINUE
ENDIF
RETURN
131 CONTINUE
C------------------- SECTION 2 ------------- coefficient = GRND1
C ** Heat transfer coeff. in the rain region
C
IF(NPATCH(1:4).EQ.'RAIN') THEN
CALL SUB2(L0CO,L0F(CO),L0HT,L0F(LBNAME('HTCF')))
L0AF=L0F(LBNAME('AIRF'))
CALL SUB2(L0GC,L0F(LBNAME('GCST')),L0BETA,L0F(LBNAME('BETA')))
CALL SUB2(L0C,L0F(C1),L0WF,L0F(LBNAME('WFLR')))
DO 1324 IX=IXF,IXL
DO 1324 IY=IYF,IYL
ICELL=IY+NY*(IX-1)
F(L0BETA+ICELL)=ABS(HRLDA*F(L0WF+ICELL)*
1 (F(L0WF+ICELL)/F(L0AF+ICELL)))**HRN
CE=(F(L0GC+ICELL)-F(L0C+ICELL))/(1-F(L0GC+ICELL))
F(L0CO+ICELL)=F(L0BETA+ICELL)*(LEWNO+CE)/(1.-F(L0C+ICELL))
F(L0HT+ICELL)=F(L0CO+ICELL)
1324 CONTINUE
ENDIF
RETURN
132 CONTINUE
C------------------- SECTION 3 ------------- coefficient = GRND2
C ** Mass transfer coeff. for rain
C
IF(NPATCH(1:4).EQ.'RAIN') THEN
CALL SUB2(L0CO,L0F(CO),L0BT,L0F(LBNAME('BETA')))
CALL SUB2(L0GC,L0F(LBNAME('GCST')),L0MT,L0F(LBNAME('MTCF')))
DO 1348 IX=IXF,IXL
DO 1348 IY=IYF,IYL
ICELL=IY+NY*(IX-1)
F(L0CO+ICELL)=F(L0BT+ICELL)/(1.-F(L0GC+ICELL))
F(L0MT+ICELL)=F(L0CO+ICELL)
1348 CONTINUE
ENDIF
RETURN
133 CONTINUE
C------------------- SECTION 4 ------------- coefficient = GRND3
C ** Resistance coeff. for W1 in the fill region
C
IF(NPATCH(1:4).EQ.'FILL') THEN
CALL SUB2(L0CO,L0F(CO),L0RH,L0F(DEN1))
CALL SUB2(L0W1,L0F(W1),L0WF,L0F(LBNAME('WFLR')))
L0AF=L0F(LBNAME('AIRF'))
DO 1331 IX=IXF,IXL
DO 1331 IY=IYF,IYL
ICELL=IY+NY*(IX-1)
F(L0CO+ICELL)=(VPLDA*(F(L0WF+ICELL)/F(L0AF+ICELL))+VPN)*
1 ABS(F(L0W1+ICELL))*F(L0RH+ICELL)*0.5/FILLH
1331 CONTINUE
ENDIF
RETURN
134 CONTINUE
C------------------- SECTION 5 ------------- coefficient = GRND4
C ** Heat transfer coeff. in the fill region
C
IF(NPATCH(1:4).EQ.'FILL') THEN
CALL SUB2(L0CO,L0F(CO),L0WF,L0F(LBNAME('WFLR')))
CALL SUB2(L0AF,L0F(LBNAME('AIRF')),L0BT,L0F(LBNAME('BETA')))
CALL SUB2(L0GC,L0F(LBNAME('GCST')),L0HT,L0F(LBNAME('HTCF')))
L0C=L0F(C1)
DO 1341 IX=IXF,IXL
DO 1341 IY=IYF,IYL
ICELL=IY+NY*(IX-1)
F(L0BT+ICELL)=ABS(HPLDA*F(L0WF+ICELL)*
1 (F(L0WF+ICELL)/F(L0AF+ICELL)))**HPN
CE=(F(L0GC+ICELL)-F(L0C+ICELL))/(1-F(L0GC+ICELL))
F(L0CO+ICELL)=F(L0BT+ICELL)*(LEWNO+CE)/(1.-F(L0C+ICELL))
F(L0HT+ICELL)=F(L0CO+ICELL)
1341 CONTINUE
ENDIF
RETURN
135 CONTINUE
C------------------- SECTION 6 ------------- coefficient = GRND5
C ** Resistance coeff. for V1,U1 & w1 in the rain region
C
IF(NPATCH(1:4).EQ.'RAIN') THEN
CALL SUB2(L0CO,L0F(CO),L0RH,L0F(DEN1))
CALL SUB2(L0W1,L0F(W1),L0WF,L0F(LBNAME('WFLR')))
CALL SUB2(L0V1,L0F(V1),L0AF,L0F(LBNAME('AIRF')))
IF(SOLVE(U1)) L0U1=L0F(U1)
DO 1366 IX=IXF,IXL
DO 1366 IY=IYF,IYL
ICELL=IY+NY*(IX-1)
if(indvar.eq.u1) vel=f(L0U1+icell)
if(indvar.eq.v1) vel=f(l0V1+icell)
if(indvar.eq.w1) vel=f(l0W1+icell)
F(L0CO+ICELL)=(VRLDA*(F(L0WF+ICELL)/F(L0AF+ICELL))+VRN)*
1 ABS(VEL)*F(L0RH+ICELL)*0.5/RAINH
1366 CONTINUE
ENDIF
RETURN
136 CONTINUE
C------------------- SECTION 7 ------------- coefficient = GRND6
C ** Resistance coeff. for U1 and V1 at the inlet
C
IF(SOLVE(U1)) THEN
C CALL ONLYIF(U1,V1,'ALL')
c IF(NPATCH(1:3).EQ.'ALL') THEN
IF(INDVAR.EQ.U1) THEN
CALL FN21(CO,DEN1,U1,0.0,VSTRT*0.5)
ELSE IF(INDVAR.EQ.V1) THEN
CALL FN21(CO,DEN1,V1,0.0,VSTRT*0.5)
ENDIF
c ENDIF
ELSE
IF(NPATCH(1:5).EQ.'INLET'.AND.INDVAR.EQ.V1) THEN
C CALL ONLYIF(V1,V1,'INLET')
CALL FN21(CO,DEN1,V1,0.0,VSTRT*0.5)
ENDIF
ENDIF
CALL FN40(CO)
RETURN
137 CONTINUE
C------------------- SECTION 8 ------------- coefficient = GRND7
C ** Resistance coeff. for W1 at the eliminator
C
IF(NPATCH(1:4).EQ.'ELIM') THEN
CALL FN21(CO,DEN1,W1,0.0,VELIM*0.5)
CALL FN40(CO)
ENDIF
RETURN
138 CONTINUE
C------------------- SECTION 9 ------------- coefficient = GRND8
C ** Mass transfer coeff. for fill
C
IF(NPATCH(1:4).EQ.'FILL') THEN
CALL SUB2(L0CO,L0F(CO),L0BT,L0F(LBNAME('BETA')))
CALL SUB2(L0GC,L0F(LBNAME('GCST')),L0MT,L0F(LBNAME('MTCF')))
DO 1347 IX=IXF,IXL
DO 1347 IY=IYF,IYL
ICELL=IY+NY*(IX-1)
F(L0CO+ICELL)=F(L0BT+ICELL)/(1.-F(L0GC+ICELL))
F(L0MT+ICELL)=F(L0CO+ICELL)
1347 CONTINUE
ENDIF
RETURN
139 CONTINUE
C------------------- SECTION 10 ------------- coefficient = GRND9
C ** Mass transfer coeff. for spray
C
IF(NPATCH(1:5).EQ.'SPRAY') THEN
CALL SUB2(L0CO,L0F(CO),L0BT,L0F(LBNAME('BETA')))
CALL SUB2(L0GC,L0F(LBNAME('GCST')),L0MT,L0F(LBNAME('MTCF')))
DO 1391 IX=IXF,IXL
DO 1391 IY=IYF,IYL
ICELL=IY+NY*(IX-1)
F(L0CO+ICELL)=F(L0BT+ICELL)/(1.-F(L0GC+ICELL))
F(L0MT+ICELL)=F(L0CO+ICELL)
1391 CONTINUE
ENDIF
RETURN
1310 CONTINUE
C------------------- SECTION 11 ------------- coefficient = GRND10
C ** Heat transfer coeff. in the spray region
C
IF(NPATCH(1:5).EQ.'SPRAY') THEN
CALL SUB2(L0CO,L0F(CO),L0WF,L0F(LBNAME('WFLR')))
CALL SUB2(L0AF,L0F(LBNAME('AIRF')),L0BT,L0F(LBNAME('BETA')))
CALL SUB2(L0GC,L0F(LBNAME('GCST')),L0HT,L0F(LBNAME('HTCF')))
L0C=L0F(C1)
DO 1392 IX=IXF,IXL
DO 1392 IY=IYF,IYL
ICELL=IY+NY*(IX-1)
F(L0BT+ICELL)=ABS(HSLDA*F(L0WF+ICELL)*
1 (F(L0WF+ICELL)/F(L0AF+ICELL)))**HSN
CE=(F(L0GC+ICELL)-F(L0C+ICELL))/(1-F(L0GC+ICELL))
F(L0CO+ICELL)=F(L0BT+ICELL)*(LEWNO+CE)/(1.-F(L0C+ICELL))
F(L0HT+ICELL)=F(L0CO+ICELL)
1392 CONTINUE
ENDIF
RETURN
1314 CONTINUE
C------------------- SECTION 15 -------------------- value = GRND3
IF(INDVAR.EQ.C1) CALL FN1(VAL,CAMB)
RETURN
1315 CONTINUE
C------------------- SECTION 16 -------------------- value = GRND4
C Mass source for moist air flow
IF(NPATCH(1:5).EQ.'MASRC') THEN
CALL SUB2(L0GC,L0F(LBNAME('GCST')),L0VAL,L0F(VAL))
CALL SUB2(L0C,L0F(C1),L0MT,L0F(LBNAME('MTCF')))
DO 1358 IX=IXF,IXL
DO 1358 IY=IYF,IYL
ICELL=IY+NY*(IX-1)
F(L0VAL+ICELL)=F(L0MT+ICELL)/(1.-F(L0C+ICELL))*
1 (F(L0GC+ICELL)-F(L0C+ICELL))
1358 CONTINUE
ENDIF
RETURN
1317 CONTINUE
C------------------- SECTION 18 -------------------- value = GRND6
C CALL ONLYIF(H1,H1,'ALL')
IF(INDVAR.EQ.H1)
1 CALL FN0(VAL,LBNAME('GHFS'))
RETURN
1318 CONTINUE
C------------------- SECTION 19 -------------------- value = GRND7
C CALL ONLYIF(C1,C1,'ALL')
IF(INDVAR.EQ.C1)
1 CALL FN0(VAL,LBNAME('GCST'))
RETURN
C***************************************************************
C
C--- GROUP 19. Special calls to GROUND from EARTH
C
19 GO TO (191,191,193,191,191,196,197,191),ISC
191 CONTINUE
RETURN
193 CONTINUE
C * ------------------- SECTION 3 ---- START OF IZ SLAB.
C ** Calculate saturation pressure by integrating Clausius-Clapeyron
C equation; Calculate saturation enthalpy and concentration
C
IF(ISWEEP.EQ.1.AND.QNE(FIINIT(P1),READFI)) THEN
C ... Set defaults
CALL FN1(DEN1,DAMB)
CALL FN1(C1,CAMB)
CALL FN1(LBNAME('TAIR'),TAMB)
CALL FN1(LBNAME('ACST'),0.)
CALL FN1(LBNAME('MTCF'),0.001)
CALL FN1(LBNAME('BETA'),0.001)
C
IF(IZ.LT.NZM+1) CALL FN1(LBNAME('TWTR'),TINI)
IF(IZ.GT.NZM) CALL FN1(LBNAME('TWTR'),TAMB)
ENDIF
IF(ISWEEP.LT.3.AND.IZ.EQ.1) THEN
L0DIF=L0F(ANYZ(LBNAME('TWTR'),2))-L0F(ANYZ(LBNAME('TWTR'),1))
IJMP=L0DIF
ENDIF
IF(IZ.GT.NZM) RETURN
CALL SUB2(L0V1,L0F(V1),L0W1,L0F(W1))
CALL SUB2(L0RH,L0F(DEN1),L0TA,L0F(LBNAME('TAIR')))
CALL SUB2(L0GH,L0F(LBNAME('GHST')),L0GC,L0F(LBNAME('GCST')))
CALL SUB2(L0AF,L0F(LBNAME('AIRF')),L0TW,L0F(LBNAME('TWTR')))
CALL SUB2(L0C,L0F(C1),L0FS,L0F(LBNAME('GHFS')))
CALL SUB2(L0AC,L0F(LBNAME('ACST')),L0ALFA,L0F(LBNAME('ALFA')))
C
DO 113 IX=1,NX
DO 113 IY=1,NYM
ICELL=IY+NY*(IX-1)
GTC=F(L0TW+ICELL)
GTC=MAX(TAMB,MIN(TINI,GTC))
GFACT1=GCO*GWH2O*(1./GTREF-1./(GTC+GTREF))/GRCONS
GFACT2=GC1*GWH2O*ALOG((GTC+GTREF)/GTREF)/GRCONS
GPSAT=GPREF*EXP(GFACT1-GFACT2)
GC2=GWH2O*GPSAT/(GWAIR*(PAMB-GPSAT)+GWH2O*GPSAT)
GSPECM=GSPECV*GC2+(1.-GC2)*GSPECG
F(L0GH+ICELL)=GSPECM*GTC+GC2*HLAT0-HAMB
C...test vapour saturation for TAIR
GTA=F(L0TA+ICELL)
GFACT1=GCO*GWH2O*(1./GTREF-1./(GTA+GTREF))/GRCONS
GFACT2=GC1*GWH2O*ALOG((GTA+GTREF)/GTREF)/GRCONS
GPSAT=GPREF*EXP(GFACT1-GFACT2)
GSAT=GWH2O*GPSAT/(GWAIR*(PAMB-GPSAT)+GWH2O*GPSAT)
F(L0AC+ICELL)=GSAT
XSAT=GSAT/(1.-GSAT)
C IF(F(L0C+ICELL).GE.GSAT.AND.(ISWEEP.EQ.LSWEEP)) THEN
C WRITE(6,*)'SATURATION AT Iz=',Iz,' XS=',XSAT
C ENDIF
GSPECM=GSPECV*F(L0C+ICELL)+(1.-F(L0C+ICELL))*GSPECG
F(L0GC+ICELL)=GC2
C ... Prevent C becoming greater than GC2 for cases with high
C ... oversaturation
F(L0C+ICELL)=AMIN1(GC2,F(L0C+ICELL))
F(L0ALFA+ICELL)=F(L0BT+ICELL)*GSPECM*LEWNO/
1 (1.-f(l0c+icell))
CE=(GC2-F(L0C+ICELL))/(1.-GC2)
GHSOR=(1.-F(L0C+ICELL))*F(L0GH+ICELL)/(1.-GC2)
DE=(HLAT0+GSPECV*GTC)
F(L0FS+ICELL)=(LEWNO*GHSOR-CE*DE*(LEWNO-1.))/(LEWNO+CE)
IF(SOLVE(U1)) THEN
L0U1=L0F(U1)
GVEL=SQRT(F(L0U1+ICELL)**2+F(L0V1+ICELL)**2+
1 F(L0W1+ICELL)**2)
ELSE
GVEL=SQRT(F(L0V1+ICELL)**2+F(L0W1+ICELL)**2)
ENDIF
F(L0AF+ICELL)=F(L0RH+ICELL)*(1.-F(L0C+ICELL))*GVEL
113 CONTINUE
IF(ISWEEP.GT.1) RETURN
L0WF=L0F(LBNAME('WFLR'))
IF(.NOT.LG(3)) THEN
DO 1003 IX=1,NX
DO 1003 IY=1,NYM
ICELL=IY+NY*(IX-1)
F(L0WF+ICELL)=WFLR
1003 CONTINUE
ELSE
IF(LG(4)) THEN
RYINN=TOPFR/NYM/2.
ELSE
RYINN=BTMFR/(NYM-ILY1)/2.
RYOUT=(TOPFR-BTMFR)/ILY1/2.
ENDIF
RYFIL=0.0
RYFIL=RYFIL+RYINN
DO 1005 IY=1,NYM
DO 1004 IX=1,NX
ICELL=IY+NY*(IX-1)
C --- Linear formula
IF(LG(7)) F(L0WF+ICELL)=0.6*WFLR*(1+RYFIL/TOPFR)
C --- Quadratic distribution
c F(L0FL+ICELL)=0.047365*WFLR*(1+(RYFIL*RYFIL)/TOPFR)
C --- Exponential 1.5
c F(L0FL+ICELL)=0.21618*WFLR*(1+(RYFIL**1.5)/TOPFR)
C --- Exponential 1.3
IF(LG(8)) F(L0WF+ICELL)=0.35252*WFLR*(1+
1 (RYFIL**1.3)/TOPFR)
1004 CONTINUE
IF(LG(4)) THEN
RYFIL=RYFIL+2.*RYINN
ELSE IF(IY.LT.(NYM-ILY1)) THEN
RYFIL=RYFIL+2*RYINN
ELSE IF(IY.EQ.(NYM-ILY1)) THEN
RYFIL=RYFIL+RYINN+RYOUT
ELSE
RYFIL=RYFIL+2.*RYOUT
ENDIF
1005 CONTINUE
ENDIF
RETURN
196 CONTINUE
C * ------------------- SECTION 6 ---- FINISH OF IZ SLAB.
IF(IZ.GT.NZM) RETURN
CALL SUB2(L0H,L0F(H1),L0C,L0F(C1))
CALL SUB2(L0GH,L0F(LBNAME('GHST')),L0GC,L0F(LBNAME('GCST')))
CALL SUB2(L0FS,L0F(LBNAME('GHFS')),L0GQ,L0F(LBNAME('GQ')))
L0HT=L0F(LBNAME('HTCF'))
DO 1961 IX=1,NX
DO 1961 IY=1,NYM
ICELL=IY+NY*(IX-1)
C ... Calculate the amount of heat transfered
F(L0GQ+ICELL)=F(L0HT+ICELL)*(F(L0FS+ICELL)-F(L0H+ICELL))
1961 CONTINUE
IF(ISWEEP.EQ.FSWEEP) RETURN
L0WF=L0F(LBNAME('WFLR'))
C ... Calculate water outlet temperature, air mass flow-rate etc.
CALL SUB2(L0RH,L0F(DEN1),L0W1,L0F(W1))
CALL SUB2(L0H,L0F(H1),L0C,L0F(C1))
L0ME=L0F(LBNAME('WEVP'))
L0AH=L0B(9)
C ... Evaporated water
if(isweep.eq.lsweep-1) then
DO 1968 IX=1,NX
DO 1968 IY=1,NYM
ICELL=IY+NY*(IX-1)
SUMEVP=SUMEVP+F(L0ME+ICELL)*F(L0AH+ICELL)
1968 CONTINUE
endif
IF(IZ.EQ.1) THEN
GSUMA=0.0
GSUMT=0.0
SUMFL=0.0
DO 1962 IX=1,NX
DO 1962 IY=1,NYM
ICELL=IY+NY*(IX-1)
GSUMA=GSUMA+F(L0AH+ICELL)
GSUMT=GSUMT+F(L0TW+ICELL)*F(L0AH+ICELL)
SUMFL=SUMFL+F(L0WF+ICELL)*F(L0AH+ICELL)
1962 CONTINUE
GTOUT=GSUMT/GSUMA
IF(SOLVE(U1)) THEN
SUMFL=SUMFL*2
ELSE
SUMFL=SUMFL*4.0*GPI
ENDIF
ELSE IF(IZ.EQ.NZM) THEN
GSUMM=0.0
GSUMH=0.0
GSUMC=0.0
GSMDAF=0.0
DO 1963 IX=1,NX
DO 1963 IY=1,NYM
ICELL=IY+NY*(IX-1)
GROW=F(L0RH+ICELL)*F(L0W1+ICELL)
GROWH=GROW*F(L0H+ICELL)
GROWC=GROW*F(L0C+ICELL)
GSUMC=GSUMC+GROWC*F(L0AH+ICELL)
GSUMH=GSUMH+GROWH*F(L0AH+ICELL)
GSUMM=GSUMM+GROW*F(L0AH+ICELL)
GSMDAF=GSMDAF+(1.-F(L0C+ICELL))*GROW*F(L0AH+ICELL)
1963 CONTINUE
IF(SOLVE(U1)) THEN
GMAIR=GSUMM*2.
GDELC=GSUMC*2.-GMAIR*CAMB
GHOUT=GSUMH/GSUMM
GQAIR=GSUMH*2.
GSMDAF=GSMDAF*2
GSUMA=GSUMA*2.
ELSE
GMAIR=GSUMM*4.0*GPI
GDELC=GSUMC*4.0*GPI-GMAIR*CAMB
GHOUT=GSUMH/GSUMM
GQAIR=GSUMH*4.0*GPI
GSMDAF=GSMDAF*4.*GPI
GSUMA=GSUMA*4*GPI
ENDIF
GQWAT=CFF*WFLR*gsuma*(TINI-GTOUT)
ENDIF
CEVP=GDELC/SUMFL
CEVP=CEVP*100.
RETURN
197 CONTINUE
C * ------------------- SECTION 7 ---- FINISH OF SWEEP.
C ...Calculate water temperature from heat balance
C
CALL SUB2(L0TW,L0F(LBNAME('TWTR')),L0GQ,L0F(LBNAME('GQ')))
CALL SUB2(L0RH,L0F(DEN1),L0GC,L0F(LBNAME('GCST')))
CALL SUB2(L0H,L0F(H1),L0C,L0F(C1))
CALL SUB2(L0ME,L0F(LBNAME('WEVP')),L0BT,L0F(LBNAME('BETA')))
CALL SUB2(L0WF,L0F(LBNAME('WFLR')),L0AC,L0F(LBNAME('ACST')))
L0ZD=L0B(15)
IZZ=IJMP*(NZM-1)
IZZST=IZZ+IJMP
DO 1972 IZ=NZM,1,-1
DO 1971 IX=1,NX
DO 1971 IY=1,NYM
ICELL=IY+NY*(IX-1)
ICL=IY+1+NY*(IX-1)+NY*NX*(IZ-1)
F(L0WF+IZZST+ICELL)=WFLR
GFL=F(L0WF+IZZ+ICELL)
F(L0TW+IZZST+ICELL)=TINI
C ... Calculate evaporation
DFL=F(L0BT+IZZ+ICELL)/(1.-F(L0C+IZZ+ICELL))/(1.-F(L0GC+
1 IZZ+ICELL))*(F(L0GC+IZZ+ICELL)-F(L0C+IZZ+ICELL))*
1 F(L0ZD+ICL)
C ... LG(6) Merkel - no evapor
IF(LG(6)) THEN
GDT=(F(L0GQ+IZZ+ICELL)*F(L0ZD+ICL))/(GFL*CFF)
DFL=0.0
ELSE
F(L0ME+IZZ+ICELL)=DFL
QEVAP=DFL*CFF*F(L0TW+IZZ+ICELL)
GDT=(F(L0GQ+IZZ+ICELL)*F(L0ZD+ICL)-QEVAP)/(GFL*CFF)
ENDIF
F(L0TW+IZZ+ICELL)=F(L0TW+IZZ+IJMP+ICELL)-GDT
F(L0TW+IZZ+ICELL)=AMAX1(AMIN1(F(L0TW+IZZ+ICELL),TINI),TAMB)
F(L0WF+IZZ+ICELL)=F(L0WF+IZZ+IJMP+ICELL)-DFL
1971 CONTINUE
IZZ=IZZ-IJMP
1972 CONTINUE
IF(MOD(ISWEEP,NPRINT).EQ.0.AND.ISWEEP.EQ.LSWEEP) THEN
C ... aditional calculation
CALL SUB2(L0TA,L0F(LBNAME('TAIR')),L0AF,L0F(LBNAME('AIRF')))
DMW=WFLR*GSUMA-SUMFL
RATIO=WFLR*GSUMA/GSMDAF
C ... Printing
CALL WRITBL
CALL WRIT40('PRINT-OUT OF MAIN TOWER-PERFORMANCE DATA')
CALL WRITBL
CALL WRIT2R('TOUTLET ',GTOUT,'TINLET ',TINI)
CALL WRIT2R('AIRMSFLX',GMAIR,'INTGRWFL',SUMFL*3.6)
CALL WRIT2R('CONLOSS ',GDELC,'WATEVAPR',CEVP)
CALL WRIT2R('AIRWRMUP',GQAIR,'WATCOOL ',GQWAT)
CALL WRIT2R('WATLOSS ',DMW, 'L/GDRY ',RATIO)
CALL WRIT2R('WTRFLRIN',WFLR*GSUMA,'WTRFLROU',SUMFL)
CALL WRIT2R('AREA ',GSUMA,'EVPLOSS ',GDELC*CFF*TINI)
c CALL WRIT2R('DRYAIRMF',GSMDAF,'SUMEVPR ',SUMEVP)
CALL WRIT1R('DRYAIRMF',GSMDAF)
ENDIF
IF(MOD(ISWEEP+1,NPRINT).EQ.0) THEN
IF(LG(1)) THEN
CALL WRITBL
CALL WRIT40('PRINT-OUT OF INITIAL DATA INPUT ')
CALL WRITBL
CALL WRIT2R('WFLR ',WFLR ,'PAMB ',PAMB)
CALL WRIT2R('TAMB ',TAMB ,'TWAMB ',TWAMB)
CALL WRIT2R('X ',XAIR ,'DAMB ',DAMB)
CALL WRIT2R('TINI ',TINI ,'WAMB ',WAMB)
CALL WRIT2R('HSLDA ',HSLDA,'HSN ',HSN)
CALL WRIT2R('VSLDA ',VSLDA,'VSN ',VSN)
CALL WRIT2R('HPLDA ',HPLDA,'HPN ',HPN)
CALL WRIT2R('VPLDA ',VPLDA,'VPN ',VPN)
CALL WRIT2R('HRLDA ',HRLDA,'HRN ',HRN)
CALL WRIT2R('VRLDA ',VRLDA,'VRN ',VRN)
CALL WRIT2R('VELIM ',VELIM,'VSTRT ',VSTRT)
IF(LG(3)) THEN
CALL WRIT40('WATER DISTRIBUTION SET TO NON-UNIFORM ')
ELSE IF(.NOT.LG(3)) THEN
CALL WRIT40('WATER DISTRIBUTION SET TO UNIFORM ')
ENDIF
ENDIF
IF(LG(2)) THEN
CALL WRITBL
CALL WRIT40('PRINT-OUT OF INPUT GEOMETRY DATA ')
CALL WRITBL
CALL WRIT2R('RAINH ',RAINH ,'INLETH ',INLETH)
CALL WRIT2R('FILLH ',FILLH ,'ELIMH ',ELIMH)
CALL WRIT2R('TOWERH ',TOWERH ,'BTMFR ',BTMFR)
CALL WRIT3R('TOPFR ',TOPFR ,'ELIMR ',ELIMR ,
1 'EXITR ',EXTR)
ENDIF
IF(ISWEEP.NE.LSWEEP-1) RETURN
C ** Set water temperature outside tower to eq. air temp. for PHOTON
IZZ=0
DO 1925 IZ=1,NZM
DO 1924 IX=1,NX
DO 1924 IY=NYM+1,NY
ICELL=IZZ+IY+NY*(IX-1)
F(L0TW+ICELL)=TAMB
1924 CONTINUE
IZZ=IZZ+IJMP
1925 CONTINUE
ENDIF
RETURN
END
C*****************************************************************
C**** PROFIL to be used for the Cooling Tower calculation
C ... to allow inlet velocity distribution appropriate
C ... for the atmospheric boundary layer; (single tower, flat-land
C ... surrounding (open terrain)).
C SUBROUTINE PROFIL sets velocity resolutes to represent a uniform
C inflow crossing a curvilinear inlet boundary. It also contains a
C coding sequence which sets an initial condition of uniform flow
C in a curvilinear grid.
C
C A limitation of this subroutine is that only one density is
C possible for all of the BFC inlet PATCHes. This should be
C rectified in the future. It should be noted that BFCA has been
C appropriated to carry the incoming density and that this could
C cause clashes with the subroutines GXLATG and GXRSET called from
C group 19, sections 3 and 6, of GREX3.
C A description of use of this option is supplied in chapter 5 of
C TR 100.
C
C***SUBROUTINE PROFIL is called from TACTGR:-
C group 1 when BFC=T, to set up the patch-wise storage for PROFIL,
C patch names starting BFC and IBFC are changed to BFT and IBFT;
C group 11 when the patch named 'IBFT', to specify uniform flow
C in a curvilinear grid;
C group 13 when the patch name begins with 'BFT', to describe a
C uniform inflow crossing a curvilinear inlet boundary.
C
C...The aerofoil Input Library case 528 and the Mizuki
C Impeller case 524 make use of it.
C
SUBROUTINE PROFIL(NX,NY,NZ)
INCLUDE 'patcmn'
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
DIMENSION EV(3),A(3),B(3),CC(3),D(3),IJKA(3),IJKB(3),IJKC(3),
1 IJKD(3)
COMMON /IDATA/IDFIL1(25),FSWEEP,IDFIL2(44),NUMPAT,IDFIL3(49)
COMMON /GENI/IGNFL1(42),NFTOT,IGNFL2(17)
COMMON /RDATA/RDFIL1(20),RHO1,RDFIL2(4),GRND,RDFIL3(59)
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
INTEGER FSWEEP
LOGICAL DONE1
SAVE
C .. Transfer wind velocity set during menu session
COMMON /RGRND/RG(100)
EQUIVALENCE (WAMB,RG(5))
C .. Next data for correction factor
DATA FCOR/1.1489E-4/
C
NAMSUB = 'PROFIL'
C*******************************************************************
C----Group 1----Section 1-------------------------------------------
C---Allocate PATCH-wise storage for BFC quantities
IF (IGR.EQ.1.AND.ISC.EQ.1) THEN
CALL GXMAKE(L0DUMM,NY*NX,'DUMM')
DO 10 IPAT = 1,NUMPAT
IF (NAMPAT(IPAT) (1:3).EQ.'BFC') THEN
C CALL MAKEPV(PVVMAS,IPAT)
C CALL MAKEPV(PVURES,IPAT)
C CALL MAKEPV(PVVRES,IPAT)
C CALL MAKEPV(PVWRES,IPAT)
NAMPAT(IPAT)(3:3)='T'
ELSEIF (NAMPAT(IPAT)(1:4).EQ.'IBFC') THEN
NAMPAT(IPAT)(4:4)='T'
ENDIF
10 CONTINUE
CALL SUB4(IJKA(2),0,IJKB(2),0,IJKC(2),1,IJKD(2),1)
CALL SUB4(IJKA(3),0,IJKB(3),1,IJKC(3),1,IJKD(3),0)
IF(NX.GT.1) LBNMUC=LBNAME('UCRT')
IF(NY.GT.1) LBNMVC=LBNAME('VCRT')
IF(NZ.GT.1) LBNMWC=LBNAME('WCRT')
ENDIF
C*******************************************************************
C---Make some initial settings for groups 11 and 13
IF (IGR.EQ.11 .OR. IGR.EQ.13) THEN
DONE1 = ISWEEP .GT. FSWEEP
C---Set unit vector for external flow direction
CALL SUB3R(UCRT,0.0,VCRT,0.0,WCRT,0.0)
IF (NX.GT.1) CALL GETCOV(NPATCH,LBNMUC,CCC,UCRT)
IF (NY.GT.1) CALL GETCOV(NPATCH,LBNMVC,CCC,VCRT)
C** Next *********************************
C .. Calculation according to ESDU - Wind eng.)
crj L0ZD=L0B(70)
crj ICL=NY+NY*NX*(IZ-1)
crj ZDIST=F(L0ZD+ICL)
ZDIST=coord1(1,ny,izstep,3)
c .. Iteration loop to calculate u*
if(iz.gt.1) goto 16
adterm=0.0
do 15 i=1,5
USTAR=WAMB/(2.5*(ALOG(10./0.1)+adterm))
C .. Height where velocity reaches constant
HEI=USTAR/6./FCOR
ratio=10./hei
adterm=5.75*RATIO-1.875*RATIO**2.-1.333*RATIO**3.+0.25*RATIO**4.
15 continue
16 continue
RATIO=ZDIST/HEI
C .. Keep value of HEI as free-stream boundary for viscosity calc.
RG(40)=HEI
C .. Effective viscosity
RG(36)=0.35*ustar*hei
C .. Velocity profile
VCRFL=2.5*USTAR*(ALOG(ZDIST/0.1)+5.75*RATIO-1.875*RATIO**2.
1 -1.333*RATIO**3.+0.25*RATIO**4.)
C .. Free stream velocity (ie. maximum velocity)
GRVL=USTAR*2.5*(ALOG(USTAR/FCOR/0.1)+1.)
IF(ZDIST.LE.HEI) THEN
VCRT=-1.*AMIN1(VCRFL,GRVL)
ELSE
VCRT=-1.*GRVL
ENDIF
C ** END ************************
IF (NZ.GT.1) CALL GETCOV(NPATCH,LBNMWC,CCC,WCRT)
IUVW = 3* (INDVAR-U1)/2
IF (IGR.EQ.11) THEN
C*****************************************************************
C----Group 11---Section 2----------------------------(GRND1 coefficient)
IF (ISC.EQ.2) THEN
C---Calculate, and set, resolute
C---IUVW=0,3 or 6 depending upon whether velocity is U1, V1 or W1
C---Calculate initial velocity component over current PATCH
CALL BFCVFX(IUVW+1,IUVW+2,IUVW+3,
1 UCRT,VCRT,WCRT,IZSTEP,VAL,L0DUMM,NX,NY)
END IF
ELSE IF (IGR.EQ.13) THEN
C*****************************************************************
C----Group 13---Section 13-----------------------------(GRND1 value)
IF (ISC.EQ.13) THEN
IMUVW = PVVMAS + (INDVAR-1)/2
IF (DONE1) THEN
C---If not the first sweep, copy mass-inflow rate or velocity
C resolute from IMUVW into VAL
CALL FNPAXY(VAL,IMUVW)
ELSE
IF (INDVAR.EQ.P1) THEN
C---If this is the first sweep of a time-step, calculate mass-
C inflow rate in VAL and store it in the PATCH-wise store PVVMAS
C---For the moment, only one density for all BFC inlets is permitted;
C this is carried in BFCA.
RHO1P = BFCA
IZ1 = IZSTEP
C---Get corner coordinates of the 4 corners of cells at the
C boundaries of the domain, at the current slab. These corners
C are labelled A, B, C and D , and the cartesian
C coordinates XC, YC and ZC for each corner are stored in
C the arrays A(3), B(3), C(3) and D(3).
IF(INTTYP.LT.2 .AND. INTTYP.GT.7) THEN
CALL WRIT40('PROFIL: CALL DISALLOWED. ')
CALL WRIT40('ONLY PATCHES OF TYPES; "EAST",...,"LOW" ')
CALL WRIT40('ARE PERMITTED. ')
CALL SET_ERR(326, 'Error. See result file.',2)
C CALL EAROUT(2)
ENDIF
C---If INTTYP is even, then PATCH-type is east, north or high so add 1
C to the number of the plane in which the PATCH lies. The arrays IJKA
C to IJKD define the 4 points, relative to (I,J,K), used to define
C the vector normal to the cell face in the PATCH.
MIT = MOD(INTTYP+1,2)
CALL SUB4(IJKA(1),MIT,IJKB(1),MIT,IJKC(1),MIT,IJKD(1),
1 MIT)
FLIO = FLOAT(2*MIT-1)
ITD2 = INTTYP/2
CALL SUB3(MITX,MOD(4-ITD2,3)+1,MITY,MOD(5-ITD2,3)+1,
1 MITZ,MOD(6-ITD2,3)+1)
CALL L0F1(VAL,ICELL,IADD,'PROFIL')
DO 20 IX = IXF,IXL
ICELL = ICELL + IADD
DO 30 IY = IYF,IYL
ICELL = ICELL + 1
CALL GETPT(IX+IJKA(MITX),IY+IJKA(MITY),
1 IZ1+IJKA(MITZ),A1,A2,A3)
CALL GETPT(IX+IJKB(MITX),IY+IJKB(MITY),
1 IZ1+IJKB(MITZ),B1,B2,B3)
CALL GETPT(IX+IJKC(MITX),IY+IJKC(MITY),
1 IZ1+IJKC(MITZ),CC1,CC2,CC3)
CALL GETPT(IX+IJKD(MITX),IY+IJKD(MITY),
1 IZ1+IJKD(MITZ),D1,D2,D3)
CALL SUB4R(A(1),A1,B(1),B1,CC(1),CC1,D(1),D1)
CALL SUB4R(A(2),A2,B(2),B2,CC(2),CC2,D(2),D2)
CALL SUB4R(A(3),A3,B(3),B3,CC(3),CC3,D(3),D3)
CALL NORML(A,B,CC,D,EV)
F(ICELL) = FLIO*RHO1P* (UCRT*EV(1)+VCRT*EV(2)+
1 WCRT*EV(3))
30 CONTINUE
20 CONTINUE
ELSEIF(INDVAR.GE.U1.OR.INDVAR.LE.W2) THEN
C.... Calculate initial velocity component over current PATCH
CALL BFCVFX(IUVW+1,IUVW+2,IUVW+3,
1 UCRT,VCRT,WCRT,IZSTEP,VAL,L0DUMM,NX,NY)
ENDIF
C.... Copy mass-inflow rate or velocity resolute, in VAL, into
C PATCH-wise storage IMUVW
CALL FNXYPA(IMUVW,VAL)
ENDIF
ENDIF
ENDIF
ENDIF
NAMSUB = 'profil'
END
C**************************************************