c

C..FILE.NAME...............GXPARA.FOR ........................020806
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c

c
c
C     SUBROUTINE GXPARA
C     The subroutine is called only when
C        PARAB=T  and NZ > 1 or
C        STEADY=F and NZ = 1
C     and IDISPA > 0,
C     in the Q1 file.
C
C     How to set IDISPA (also IDISPB & IDISPC) :
C     IDISPA = 1  elicits the dump of all the z- (or t-) direction slabs
C                 between IZ (or ISTEP) = IDISPB and IDISPC.
C                 limit the bounds for the first slab and the last slab
C                 to be dumped (defaults: IDISPB=1; IDISPC=NZ).
C     IDISPA > 1  it will 'group' IDISPA cells into one and assign the
C                 value of the central cell to the resulting cell, as
C                 shown in the drawing below. Therefore, the number
C                 of the 'group' cell should be odd number.
C                 Even IDISPA's set by users are automatically converted
C                 to the next odd numbers by the program.
C                 Variables IDISPB and IDISPC set the limits in which the
C                 above grouping will happen. It may be noted that this
C                 trade-off between simplicity and generality will
C                 result in the omission of a couple of slabs, which,
C                 being within the bounds defined by IDISPB and IDISPC,
C                 are not enough to form a new 'group'. The adoption
C                 of this practice is justified by the fact that a
C                 large number of slabs are normally used in a parabolic
C                 simulation of flow.
C
C                  -----------------------------
C                  :        :        :         :
C                  -----------------------------       IDISPA=3
C                  :***************************:
C                  :*  1    :  2-->  :  3     *:       ---- old grid
C                  :***************************:
C                  -----------------------------       **** new cell
C                  :***************************:
C                  :*  1    :  2-->  :  3     *:        : Y
C                   ***************************:        :
C                  -----------------------------        :---> Z (or time)
C                  :        :        :         :
C                  -----------------------------
C
C   a)   If both XULAST and YVLAST are constant and DZ for all
C        IZ is known at the start of the computation (i.e AZDZ=0.0)
C        then only a PARADA (formatted) or PARPHI (direct-access)
C        is written in addition to the original phi/phida file.
C
C   b)   If XULAST and/or YVLAST are not functions of IZ, but
C        AZDZ is set, a scratch file named PARAXX is created
C        during execution and its contents written to a PARADA
C        (formatted) or PARPHI (direct-access) in addition to the
C        original phi/phida file at the end of the run.
C
C   c)   If either XULAST, or YVLAST are functions of IZ then a
C        PARZDA (formatted) or PARXYZ (direct-access) file is
C        written in addition to PARDA or PARPHI.
C        PHOTON treats these files as if results of a BFC case.
C
C.... Formatted or direct-access files are created by GXPARA according
C     to whether PHIDA=T or F in PREFIX. Formatted files are named
C     PARPHI and PARXYZ, direct-access files are PARADA and PARZDA.
C
C.... If you do not wish to dump the values of all variables stored,
C     set the names of variables which you do NOT want as a 4-character
C     name, the last 2 characters of which are '**', in the Q1 file;
C     eg  NAME(R2)=R2**
C
C     (original PHIDA or PHI should be used for restart runs,
C                        NOT  PARADA or PARPHI)
C--------------------------------------------------------------
      SUBROUTINE GXPARA
      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'
      COMMON /PARA1/ MYNAME(150) /PARA2/ MYSTOR(150)
      LOGICAL BFGRID, ZVGRID, MYSTOR, PAIRS, DAFILE, NEZ
      CHARACTER*(4) MYNAME, MYVNAM
      COMMON/LUNITS/LUNIT(60)
      COMMON/NAMFN/NAMFUN,NAMSUB
      COMMON/AUXDAT/ZCCNF
      CHARACTER*6 NAMFUN,NAMSUB
      SAVE DAFILE,IREPHI,IREXYZ,NLUPHI,NLUXYZ,NLUSCR,MYFIZ,MYLIZ,
     1NZEFF,PAIRS,ZVGRID,BFGRID,IFREQ,IFST,IDUMP,LNRC,IREC,IXYZ,
     1ZWLST,ISCR,L0ZT,L0AUX
      SAVE IZVG
C
      NAMFUN='GXPARA'
      NAMSUB=NAMFUN
C*****************************************************************
C.... Determine whether z or time is the parabolic dimension
      IF(STEADY) THEN
        KZTFR=KZFR
        NZT=NZSTEP
        IZT=IZ
      ELSE
        IF(NZ.GT.1) THEN
          WRITE(LUPR1,*) 'GXPARA cannot be used for transients if'
          WRITE(LUPR1,*) 'NZ exceeds 1. IDISPA has been set = 0  '
          IDISPA=0
          RETURN
        ENDIF
        KZTFR=KTFR
        NZT=LSTEP
        IZT=ISTEP
      ENDIF
      NZPR=NZT
C.... Group 1 settings
      IF(IGR.EQ.1) THEN
        IFREQ=1
        IF(STEADY) THEN
          ZWLST=ZWLAST
        ELSE
          ZWLST=TLAST
        ENDIF
C.... If PHIDA=T in PREFIX, DAFILE is set true.
        DAFILE=.FALSE.
        IPHI=ITWO(22,23,LUNIT(23).EQ.0)
        IF(IPHI.EQ.23) DAFILE=.TRUE.
C
C.... Logical units for output files
        IF(.NOT.DAFILE) THEN
          IREPHI=31
          IREXYZ=32
        ELSE
          IREPHI=35
          IREXYZ=36
        ENDIF
C
        NLUXYZ=LUNIT(IREXYZ)
        NLUPHI=LUNIT(IREPHI)
        NLUSCR=LUNIT(33)
C.... See first and last slab
        MYFIZ=MAX0(1,IDISPB)
        MYLIZ=NZT
        IF(IDISPC.GT.0.AND.IDISPC.LE.NZT) MYLIZ=IDISPC
        NZEFF=MYLIZ-MYFIZ+1
        CALL GXMAKE(L0ZT,NZT,'ZorT')
        CALL GXMAKE(L0AUX,MAX0((NX+1)*(NY+1),NZEFF),'IAUX')
C.... Grid types
        PAIRS=.FALSE.
        ZVGRID=.FALSE.
        BFGRID=.FALSE.
        IF(NEZ(AZXU).OR.NEZ(AZYV)) THEN ! grid expands, so use BFC treatment
          BFGRID=.TRUE.                 ! for PHOTON display
        ELSEIF(NEZ(AZDZ)) THEN
          ZVGRID=.TRUE.
        ELSEIF(F(KZTFR+1).LT.0.0) THEN  ! use method of pairs  for z grid,
C     store high face coords f-array
          PAIRS=.TRUE.
          IZGR=0
          DO 101 J=1,NZT/2
            IF(IZGR.LT.NZT) THEN
              IZFR=NINT(F(KZTFR+2*J-1))
              IZFR=IABS(IZFR)
              IF(IZGR.EQ.0) THEN
                F(L0ZT+1)=F(KZTFR+2*J)
                IF(IZFR.LE.NZPR) THEN
                  DO 102 I=2,IZFR
                    F(L0ZT+I)=F(KZTFR+2*J)+F(L0ZT+I-1)
 102              CONTINUE
                ENDIF
              ELSE
                IF(IZFR+IZGR.LE.NZPR) THEN
                  DO 103 I=1,IZFR
                    F(L0ZT+I+IZGR)=F(KZTFR+2*J)+F(L0ZT+I+IZGR-1)
 103              CONTINUE
                ENDIF
              ENDIF
              IZGR=IZGR+IZFR
            ENDIF
 101      CONTINUE
          IF(IZGR.GT.NZPR) THEN
            NZT=NZPR+1
            GO TO 210
          ENDIF
        ENDIF
C.... Determine slab frequency for dumping to parphi
        IF(IDISPA.GT.1) THEN
          IFREQ=IDISPA
          IFST=IDISPA/2+1
C.... Check for odd no.
          IF(MOD(IDISPA,2).EQ.0) IFREQ=IDISPA+1
          NZEFF=NZEFF/IFREQ
          IDUMP=0
        ELSE
          IFREQ=1
        ENDIF
  210 CONTINUE
C.... End of Group 1 settings
      ENDIF
C*****************************************************************
      IF(IGR.EQ.19) THEN
C
C.... Not if iz < myfiz or iz > myliz
        IF((IZT.GE.MYFIZ).AND.(IZT.LE.MYLIZ)) THEN
C
C.... F-array zero-location indices of geometrical quantities
          L0X=L0F(XU2D)
          L0Y=L0F(YV2D)
          IF(.NOT.CARTES) L0R=L0F(RV2D)
          IF(BFGRID .AND. .NOT.CARTES .AND. NX.GT.1) THEN
            L0SP=L0F(EASP2)
            L0XG=L0F(XG2D)
          ENDIF
C
C.... Part 1: ....................................................
C     First-slab practices for xyz and phi
C
          IF(IZT.EQ.MYFIZ) THEN
C
C.... 1a. open files
C
            if(debug) call writ1i('open  fl',irephi)
            CALL OPENFL(IREPHI)
            IF(BFGRID) CALL OPENFL(IREXYZ)
            IF(ZVGRID) CALL OPENFL(33)
C
C.... Set record length for PHI or PHIDA files
            LNRC=LENREC
            IF(.NOT.DAFILE) LNRC=6
C
            RRINN=RINNER
            IF(BFGRID) RRINN=0.0
            IF(.NOT.DAFILE) THEN
C
C.... 1b. nx, ny, nz to xyz
C
              IF(BFGRID) WRITE (NLUXYZ, '(3I5)' ) NX+1, NY+1, NZEFF+1
C
C.... 1c. title and other preliminaries to PHI
C
              WRITE(NLUPHI,201) MESS(1:40)
  201         FORMAT(1X,A40,9A4)
              WRITE (NLUPHI,'(1X,79L1)') CARTES,ONEPHS,BFGRID,XCYCLE
              WRITE (NLUPHI,'(1X,7I10)') NX,NY,NZEFF,NPHI,DEN1,DEN2,
     1                              EPOR,NPOR,HPOR,VPOR,LENREC
              WRITE (NLUPHI,'(1PE13.6)') RRINN
            ELSE
C
C.... 1b. nx, ny, nz to xyz
C
              IF(BFGRID) WRITE (NLUXYZ, REC=1 ) INT4(NX+1), INT4(NY+1),
     1           INT4(NZEFF+1)
C
C.... 1c. title and other preliminaries to phi
C
              WRITE (NLUPHI,REC=1) MESS(1:40)
              WRITE (NLUPHI,REC=2) LOGICAL(CARTES,4),
     1           LOGICAL(ONEPHS,4),LOGICAL(BFGRID,4),LOGICAL(XCYCLE,4)
              WRITE (NLUPHI,REC=3) INT4(NX),INT4(NY),INT4(NZEFF),
     1           INT4(NPHI),INT4(DEN1),INT4(DEN2),INT4(EPOR),
     1           INT4(NPOR),INT4(HPOR),INT4(VPOR),INT4(LENREC)
              WRITE (NLUPHI,REC=4) RRINN
            ENDIF
C
C.... 1d. variable names to phi
C           NB: name of velocities overwritten if bfgrid
C           (This is because the default vector components
C           used by PHOTON are the cartesian components.)
C
            DO 10 II=1, NPHI
              MYNAME (II) = NAME (II)
  10        CONTINUE
            DO 110 JJ=1, NPHI
              MYSTOR(JJ)=.FALSE.
              MYVNAM=MYNAME(JJ)
              IF(STORE(JJ)) then
                IF(MYVNAM(3:4).NE.'**') THEN
                  IF(MOD(IPRN(JJ),5).EQ.0) MYSTOR(JJ)=.TRUE.
                ENDIF
              ENDIF
  110       CONTINUE
            IF(BFGRID) THEN
              IF(MYSTOR(U1)) MYNAME (U1) = 'UCRT'
              IF(MYSTOR(V1)) MYNAME (V1) = 'VCRT'
              IF(MYSTOR(W1)) MYNAME (W1) = 'WCRT'
              IF(MYSTOR(U2)) MYNAME (U2) = 'U2CR'
              IF(MYSTOR(V2)) MYNAME (V2) = 'V2CR'
              IF(MYSTOR(W2)) MYNAME (W2) = 'W2CR'
            ENDIF
            IF(.NOT.DAFILE) THEN
              NREC=NUMREC (NPHI,19)
              I1=-18
              I2=0
              DO 11 I=1, NREC
                I1=I1+19
                I2=MIN0 (I2+19, NPHI)
                WRITE (NLUPHI,'(1X,19A4)') (MYNAME (J), J=I1, I2)
  11          CONTINUE
            ELSE
              WRITE (NLUPHI,REC=5) (MYNAME (J), J=1,NPHI)
            ENDIF
C
C.... 1e. cell faces (first x, then y, then z) to PHI
C     NB dummy values written if bfgrid since geometry file (xyz) is used
C
            IF(DAFILE) IREC=5
            DO 12 KK=1, 3
              IF(KK.EQ.1) THEN
                NN=NX
              ELSEIF(KK.EQ.2) THEN
                NN=NY
              ELSEIF(KK.EQ.3) THEN
                NN=NZEFF
              ENDIF
              NREC=NUMREC (NN,LNRC)
              I1=1-LNRC
              I2=0
              IF(KK.EQ.3.AND.ZVGRID) IZVG=IREC
              DO 13 I=1, NREC
                I1=I1+LNRC
                I2=MIN0 (I2+LNRC, NN)
                IF(.NOT.BFGRID.AND..NOT.ZVGRID.AND.KK.EQ.3) THEN
                  IF(.NOT.PAIRS.AND.MYFIZ.GT.1) F(L0ZT+MYFIZ-1)=
     1                 F(KZTFR+MYFIZ-1)
                  DO 15 J=I1,I2
                    IF(.NOT.PAIRS) F(L0ZT+J*IFREQ+MYFIZ-1)=F(KZTFR+
     1                      J*IFREQ+MYFIZ-1)
                    IF(MYFIZ.GT.1) F(L0ZT+J*IFREQ+MYFIZ-1)=
     1                      F(L0ZT+J*IFREQ+MYFIZ-1)-F(L0ZT+MYFIZ-1)
 15               CONTINUE
                ENDIF
C.... Formatted files
                IF(.NOT.DAFILE) THEN
                  IF(BFGRID) THEN
                    WRITE (NLUPHI,1962) (.9999, JJ=I1, I2)
                  ELSE
                    IF(KK.EQ.1) WRITE (NLUPHI,1962)
     1                              (F(L0X+(JJ-1)*NY+1),JJ=I1,I2)
                    IF(KK.EQ.2) WRITE (NLUPHI,1962)
     1                   (F(L0Y+JJ),JJ=I1,I2)
                    IF((.NOT.ZVGRID).AND.(KK.EQ.3)) WRITE(NLUPHI,1962)
     1                 (F(L0zt+JJ*IFREQ+MYFIZ-1)*ZWLST,JJ=I1,I2)
                  ENDIF
                ELSE
C.... Direct-access files
                  IREC=IREC+1
                  IF(BFGRID) THEN
                    WRITE (NLUPHI,REC=IREC) (.9999, JJ=I1, I2)
                  ELSE
                    IF(KK.EQ.1) THEN
                      WRITE (NLUPHI,REC=IREC)
     1                 (F(L0X+(JJ-1)*NY+1),JJ=I1,I2)
                    ELSEIF(KK.EQ.2) THEN
                      WRITE (NLUPHI,REC=IREC) (F(L0Y+JJ),JJ=I1,I2)
                    ELSEIF(.NOT.ZVGRID.AND.KK.EQ.3) THEN
                      WRITE(NLUPHI,REC=IREC)
     1                 (F(L0zt+JJ*IFREQ+MYFIZ-1)*ZWLST,JJ=I1,I2)
                    ENDIF
                  ENDIF
                ENDIF
  13          CONTINUE
  12        CONTINUE
C
C.... 1f. pressure correction at each slab to phi
C         (the file will never be used for restarts, so
C         dummy values are written.)
C
            IF(.NOT.ZVGRID) THEN
              NREC=NUMREC (NZEFF,LNRC)
              I1=1-LNRC
              I2=0
              DO 14 I=1, NREC
                I1=I1+LNRC
                I2=MIN0 (I2+LNRC, NZEFF)
                IF(.NOT.DAFILE) THEN
                  WRITE (NLUPHI, 1962) (.9999, JJ=I1, I2)
                ELSE
                  IREC=IREC+1
                  WRITE (NLUPHI, REC=IREC) (.9999, JJ=I1, I2)
                ENDIF
 14           CONTINUE
            ENDIF
C
C.... 1g. and store(phii) to phi
C
            IF(.NOT.ZVGRID) THEN
              IF(.NOT.DAFILE) THEN
                WRITE (NLUPHI,'(1X,79L1)') (MYSTOR(I), I=1, NPHI)
              ELSE
                IREC=IREC+1
                WRITE(NLUPHI,REC=IREC) (LOGICAL(MYSTOR(I),4),I=1,NPHI)
              ENDIF
            ENDIF
          ENDIF
C
C.... Part 2: ....................................................
C     2a.  slab values of each stored var to phi
C          NB: u and v components are converted to cartesian
C          components if bfgrid.and..not.cartes.and.(nx.gt.1)
C
C.... Variables only dumped every ifreq slabs
C
          IWRIT=IFST+IDUMP*IFREQ+MYFIZ-1
          IF(ZVGRID.AND.IZT.EQ.MYFIZ) ISCR=0
          IF((IFREQ.EQ.1).OR.(IWRIT.EQ.IZT)) THEN
            DO 20 KK=1, NPHI
              IF(MYSTOR(KK)) THEN
                L0VAR=L0F(KK)
C
C.... This is the conversion mentioned above; the cartesian
C     components are stored in the spare array EASP2.
C
                IF(BFGRID.AND.(.NOT.CARTES).AND.(NX.GT.1).AND
     1            .((KK.EQ.U1).OR.(KK.EQ.V1).OR.(KK.EQ.U2).OR.
     1            (KK.EQ.V2))) THEN
                  IF((KK.EQ.U1).OR.(KK.EQ.U2)) THEN
                    L0ANG=L0X
                  ELSE
                    L0ANG=L0XG
                  ENDIF
                  DO 21 ICELL=1, NX*NY
                    F(L0SP+ICELL)=F(L0VAR+ICELL)*COS(F(L0ANG+ICELL))
 21               CONTINUE
                  L0VAR=L0SP
                ENDIF
                LENPHI=LNRC
                IF(ZVGRID) LENPHI=LENREC
                NREC=NUMREC (NX*NY,LENPHI)
                I1=1-LENPHI
                I2=0
                DO 22 I=1, NREC
                  I1=I1+LENPHI
                  I2=MIN0 (I2+LENPHI, NX*NY)
C.... Formatted files
                  IF(.NOT.DAFILE.AND..NOT.ZVGRID) THEN
                    WRITE (NLUPHI, 1962) (F(L0VAR+IJK), IJK=I1, I2)
                  ELSE
C.... PHIDA file
                    IF(.NOT.ZVGRID) THEN
                      IREC=IREC+1
                      WRITE (NLUPHI,REC=IREC)
     1                  (F(L0VAR+IJK), IJK=I1, I2)
C.... Scratch file
                    ELSE
                      ISCR=ISCR+1
                      WRITE (NLUSCR,REC=ISCR)
     1                  (F(L0VAR+IJK), IJK=I1, I2)
                    ENDIF
                  ENDIF
 22             CONTINUE
              ENDIF
 20         CONTINUE
            IDUMP=IDUMP+1
          ENDIF
C
C.... 2b. If zvgrid, store the high-face coordinate
C
          IF(ZVGRID) THEN
            IF(IZT.EQ.MYFIZ) THEN
              F(L0AUX+1)=DZ
            ELSE
              F(L0AUX+IZT-MYFIZ+1)=F(L0AUX+IZT-MYFIZ)+DZ
            ENDIF
          ENDIF
C
C.... Part 3: ....................................................
C     cell-corner coordinates to xyz if bfgrid
C
          IF(BFGRID) THEN
            IF(DAFILE.AND.IZT.EQ.MYFIZ) IXYZ=1
            IGEOM=(IDUMP)*IFREQ+MYFIZ
            IF((IFREQ.EQ.1).OR.(IGEOM.EQ.IZT)) THEN
C
C.... 3a. x coordinates
C
              IA=1
              DO 30 II=1, NY+1
                F(L0AUX+IA)=0.0
                IA=IA+1
 30           CONTINUE
              IF(CARTES) THEN
                DO 31 JJ=1,NX
                  DO 131 II=1,NY+1
                    F(L0AUX+IA)=F(L0X+NY*(JJ-1)+1)
                    IA=IA+1
 131              CONTINUE
 31             CONTINUE
              ELSE
                DO 32 II=1, NX
                  F(L0AUX+IA)=RINNER*SIN(F(L0X+NY*(II-1)+1))
                  IA=IA+1
                  DO 33 JJ=1, NY
                    IC=(II-1)*NY+JJ
                    F(L0AUX+IA)=F(L0R+IC)*SIN(F(L0X+IC))
                    IA=IA+1
 33               CONTINUE
 32             CONTINUE
              ENDIF
              IF(.NOT.DAFILE) THEN
                WRITE (NLUXYZ, 1961) (F(L0AUX+II),II=1,(NX+1)*(NY+1))
              ELSE
                NREC=NUMREC((NY+1)*(NX+1),LENREC)
                IDIR=0
                IF(IZT.NE.MYFIZ) IXYZ=IXYZ+NREC-1
                CALL WRTBFG(LENREC,NREC,IXYZ,IDIR,NLUXYZ,(NX+1)*(NY+1),
     1                      L0AUX)
              ENDIF
C
C.... 3b. y coordinates
C
              IA=1
              IF(CARTES) THEN
                DO 38 II=1,NX+1
                  F(L0AUX+IA)=0.0
                  IA=1+IA
                  DO 34 JJ=1, NY
                    F(L0AUX+IA)=F(L0Y+JJ)
                    IA=IA+1
 34               CONTINUE
 38             CONTINUE
              ELSE
                F(L0AUX+1)= RINNER
                IA=IA+1
                DO 35 JJ=1, NY
                  F(L0AUX+IA) = F(L0R+JJ)
                  IA=IA+1
 35             CONTINUE
                DO 36 II=1, NX
                  F(L0AUX+IA)=RINNER*COS(F(L0X+NY*(II-1)+1))
                  IA=IA+1
                  DO 37 JJ=1, NY
                    IC=(II-1)*NY+JJ
                    F(L0AUX+IA)=F(L0R+IC)*COS(F(L0X+IC))
                    IA=IA+1
 37               CONTINUE
 36             CONTINUE
              ENDIF
              IF(.NOT.DAFILE) THEN
                WRITE (NLUXYZ, 1961) (F(L0AUX+II),II=1,(NX+1)*(NY+1))
              ELSE
                NREC=NUMREC((NY+1)*(NX+1),LENREC)
                IDIR=0
                IXYZ=IXYZ+NREC-1
                CALL WRTBFG(LENREC,NREC,IXYZ,IDIR,NLUXYZ,(NX+1)*(NY+1),
     1                      L0AUX)
              ENDIF
            ENDIF
C
C.... 3c. z coordinates
C
            IF(IZT.EQ.MYFIZ) THEN
              ZCCNF=0.0
            ELSE
              ZCCNF=ZCCNF+DZ
            ENDIF
            IF((IFREQ.EQ.1).OR.(IGEOM.EQ.IZT)) THEN
              IF(.NOT.DAFILE) THEN
                WRITE (NLUXYZ,1961) (ZCCNF, II=1, (NX+1)*(NY+1))
              ELSE
                NREC=NUMREC((NY+1)*(NX+1),LENREC)
                IDIR=1
                IXYZ=IXYZ+NREC-1
                CALL WRTBFG(LENREC,NREC,IXYZ,IDIR,NLUXYZ,(NX+1)*(NY+1),
     1                      L0AUX)
              ENDIF
            ENDIF
          ENDIF
C
C.... Part 4: ....................................................
C     last-slab practices.
C
C.... 4a. If zvgrid, the phi file is completed here, and the
C         scratch file is closed and deleted.
C
          IF(ZVGRID.AND.(IZT.EQ.MYLIZ)) THEN
            NREC=NUMREC (NZEFF,LNRC)
            I1=1-LNRC
            I2=0
            DO 40 I=1, NREC
              I1=I1+LNRC
              I2=MIN0 (I2+LNRC, NZEFF)
              IF(.NOT.DAFILE) THEN
                WRITE (NLUPHI,1962) (F(L0AUX+JJ*IFREQ), JJ=I1,I2)
              ELSE
                IREC=IZVG+1
                WRITE(NLUPHI,REC=IREC) (F(L0AUX+JJ*IFREQ),JJ=I1,I2)
              ENDIF
 40         CONTINUE
C
C.... Now pressure correction
C
            NREC=NUMREC (NZEFF,LNRC)
            I1=1-LNRC
            I2=0
            DO 41 I=1, NREC
              I1=I1+LNRC
              I2=MIN0 (I2+LNRC, NZEFF)
              IF(.NOT.DAFILE) THEN
                WRITE (NLUPHI, 1962) (.9999, JJ=I1, I2)
              ELSE
                IREC=IREC+1
                WRITE (NLUPHI,REC=IREC) (.9999, JJ=I1,I2)
              ENDIF
 41         CONTINUE
C
C.... Now store (phi)
C
            IF(.NOT.DAFILE) THEN
              WRITE (NLUPHI,'(1X,79L1)') (MYSTOR(I), I=1, NPHI)
            ELSE
              IREC=IREC+1
              WRITE(NLUPHI,REC=IREC) (LOGICAL(MYSTOR(I),4),I=1,NPHI)
            ENDIF
C
C.... Now the scratch file (da file)
C
c            CALL CLOSFL(33)
c            CALL OPENFL(34)
            ISCR=0
            DO 42 IIZZ=1, NZEFF
              DO 43 KK=1, NPHI
                IF(MYSTOR(KK)) THEN
                  NREC=NUMREC (NX*NY,LENREC)
                  I1=1-LENREC
                  I2=0
                  DO 44 I=1, NREC
                    I1=I1+LENREC
                    I2=MIN0 (I2+LENREC, NX*NY)
                    ISCR=ISCR+1
                    READ(NLUSCR,REC=ISCR) (F(L0AUX+IJK), IJK=I1,I2)
 44               CONTINUE
                  NREC=NUMREC (NX*NY,LNRC)
                  I3=1-LNRC
                  I4=0
                  DO 144 I=1, NREC
                    I3=I3+LNRC
                    I4=MIN0 (I4+LNRC, NX*NY)
                    IF(.NOT.DAFILE) THEN
                      WRITE (NLUPHI,1962) (F(L0AUX+IJK), IJK=I3,I4)
                    ELSE
                      IREC=IREC+1
                      WRITE(NLUPHI,REC=IREC) (F(L0AUX+IJK),IJK=I3,I4)
                    ENDIF
 144              CONTINUE
                ENDIF
 43           CONTINUE
 42         CONTINUE
C
cc            CALL CLOSFL(34)
            CALL CLOSFL(33)
          ENDIF
C
C.... 4b. write the high face to xyz.
C
          IF((IZT.EQ.MYLIZ).AND.BFGRID) THEN
C
C.... x coordinates
            IA=1
            DO 45 II=1, NY+1
              F(L0AUX+IA)=0.0
              IA=IA+1
 45         CONTINUE
            IF(CARTES) THEN
              DO 46 JJ=1,NX
                DO 146 II=1,NY+1
                  F(L0AUX+IA)=F(L0X+NY*(JJ-1)+1)
                  IA=IA+1
 146            CONTINUE
 46           CONTINUE
            ELSE
              DO 47 II=1, NX
                F(L0AUX+IA)=RINNER*SIN(F(L0X+NY*(II-1)+1))
                IA=IA+1
                DO 48 JJ=1, NY
                  IC=(II-1)*NY+JJ
                  F(L0AUX+IA)=F(L0R+IC)*SIN(F(L0X+IC))
                  IA=IA+1
 48             CONTINUE
 47           CONTINUE
            ENDIF
            IF(.NOT.DAFILE) THEN
              WRITE (NLUXYZ, 1961) (F(L0AUX+II),II=1,(NX+1)*(NY+1))
            ELSE
              NREC=NUMREC((NY+1)*(NX+1),LENREC)
              IDIR=0
              IXYZ=IXYZ+NREC-1
              CALL WRTBFG(LENREC,NREC,IXYZ,IDIR,NLUXYZ,(NX+1)*(NY+1),
     1                    L0AUX)
            ENDIF
C
C.... y coordinates
C
            IA=1
            IF(CARTES) THEN
              DO 49 II=1,NX+1
                F(L0AUX+IA)=0.0
                IA=1+IA
                DO 149 JJ=1, NY
                  F(L0AUX+IA)=F(L0Y+JJ)
                  IA=IA+1
 149            CONTINUE
 49           CONTINUE
            ELSE
              F(L0AUX+1)= RINNER
              IA=IA+1
              DO 410 JJ=1, NY
                F(L0AUX+IA) = F(L0R+JJ)
                IA=IA+1
 410          CONTINUE
              DO 411 II=1, NX
                F(L0AUX+IA)=RINNER*COS(F(L0X+NY*(II-1)+1))
                IA=IA+1
                DO 412 JJ=1, NY
                  IC=(II-1)*NY+JJ
                  F(L0AUX+IA)=F(L0R+IC)*COS(F(L0X+IC))
                  IA=IA+1
 412            CONTINUE
 411          CONTINUE
            ENDIF
            IF(.NOT.DAFILE) THEN
              WRITE (NLUXYZ, 1961) (F(L0AUX+II),II=1,(NX+1)*(NY+1))
            ELSE
              NREC=NUMREC((NY+1)*(NX+1),LENREC)
              IDIR=0
              IXYZ=IXYZ+NREC-1
              CALL WRTBFG(LENREC,NREC,IXYZ,IDIR,NLUXYZ,(NX+1)*(NY+1),
     1                    L0AUX)
            ENDIF
C
C.... z coordinates
C
            ZCCNF=ZCCNF+DZ
            IF(.NOT.DAFILE) THEN
              WRITE (NLUXYZ,1961) (ZCCNF,II=1,(NX+1)*(NY+1))
            ELSE
              NREC=NUMREC((NY+1)*(NX+1),LENREC)
              IDIR=1
              IXYZ=IXYZ+NREC-1
              CALL WRTBFG(LENREC,NREC,IXYZ,IDIR,NLUXYZ,(NX+1)*(NY+1),
     1                    L0AUX)
            ENDIF
C
C.... 4c. close PHI and XYZ
C
            CALL CLOSFL(IREXYZ)
          ENDIF
          IF(IZT.EQ.MYLIZ) then
            if(debug) call writ1i('close fl',irephi)
            CALL CLOSFL(IREPHI)
          endif
        ENDIF
      ENDIF
C
 1961 FORMAT (5(1PE13.6))
 1962 FORMAT (6(1PE13.6))
      NAMFUN='gxpara'
      NAMSUB=NAMFUN
      END
C-----------------------------------------------------------------------
c
      SUBROUTINE WRTBFG(LENREC,N,IRXY,IDIR,NLU,NXNY1,L0AUX)
      INCLUDE '/phoenics/d_includ/farray'
      COMMON/AUXDAT/ZCCNF
      I1=1-LENREC
      I2=0
      DO 10 I=1,N
        IRXY=IRXY+1
        I1=I1+LENREC
        I2=MIN0(I2+LENREC,NXNY1)
        IF(IDIR.EQ.1) THEN
          WRITE (NLU,REC=IRXY) (ZCCNF,KK=I1,I2)
        ELSEIF(IDIR.EQ.0) THEN
          WRITE (NLU,REC=IRXY) (F(L0AUX+KK),KK=I1,I2)
        ENDIF
 10   CONTINUE
      END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
C     SUBROUTINE GXPVAD
C     The subroutine is called only when
C        PARAB=T  and U1AD, U2AD, V1AD or V2AD = GRND1
C     in the Q1 file.
C
C.... The subroutine adds to the cross stream velocity (V1 say)
C     a component equal  -W1*TAN(ALF) where ALF is the grid angle
C.    at the velocity location. It is needed for parabolic
C     contracting and expanding grids, i.e. grids for which
C     constant-y or constant-x lines are not precisely at right
C     angles to constant-z lines. It is essential when using such
C     grids with the IPARAB=4 (wholly supersonic flow) or IPARAB=5
C     (supersonic flow with allowance for a subsonic free stream)
C     options.
C
C-----------------------------------------------------------------------
      SUBROUTINE GXPVAD(XLAST,XR,XU,VELZ,CDIR)
      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'
C
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
      CHARACTER*1 CDIR
      INTEGER VELZ,XU
C
      NAMFUN='GXPVAD'
      NAMSUB=NAMFUN
c
c.... IYL is set to NY-1 (and later re-set to NY) to reflect the fact that
c     there is no flow through the outer boundary
      IF(CDIR.EQ.'X') THEN
        IXL=NX-1
      ELSE
        IYL=NY-1
      ENDIF
      GDXDZ=XLAST*(1.-1./XR)/DZ
      CALL FN2(VELAD,XU,0.0, -GDXDZ/XLAST)
      CALL FN26(VELAD,VELZ)
      IF(CDIR.EQ.'X') THEN
        IXL=NX
      ELSE
        IYL=NY
      ENDIF
      NAMFUN='gxpvad'
      NAMSUB=NAMFUN
      END
c