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 '/phoenics/d_includ/patcmn'
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/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 '/phoenics/d_includ/patcmn'
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoenics/d_includ/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
**************************************************