c
C..FILE.NAME...............GXPARA.FOR ........................160414
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 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE '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,BUFF(4)*256
REAL(4) ARY4(LENREC)
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.
C... IFREQ is odd so that the middle slab (of band) saved for plotting purposes
IF(IDISPA>2 .AND. MOD(IDISPA,2).EQ.0) IFREQ=IDISPA+1
NZEFF=NZEFF/IFREQ
IDUMP=0
ELSE
IFREQ=1
IFST=0; IDUMP=0
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
100 CONTINUE
if(debug) call writ1i('open fl',irephi)
CALL OPENZZ(IREPHI,IOS)
IF(IOS/=0) THEN
CALL CLOSZZ(IREPHI) ! Close PHI(DA) file then post query about action
C Likely IOSTAT errors are
IF(IOS==38) THEN ! 38 severe: Error during write
BUFF(1)='Error while writing PHI(DA) file'
ELSEIF(IOS==608) THEN ! 608 severe: No space left on device
BUFF(1)='No space left on device to write PHI(DA) file'
ELSE
CALL IOEMZZ(IOS,BUFF(1))
ENDIF
BUFF(2)='Close any open files or make more space.'
BUFF(3)='Press OK to retry saving file or Cancel to'
BUFF(4)='continue without saving.'
CALL ERRMSG(BUFF,4,222,IEND)
IF(IEND==0) GOTO 100
RETURN
ENDIF
IF(BFGRID) THEN
111 CONTINUE
CALL OPENZZ(IREXYZ,IOS)
IF(IOS/=0) THEN
CALL CLOSZZ(IREXYZ) ! Close XYZ(DA) file then post query about action
C Likely IOSTAT errors are
IF(IOS==38) THEN ! 38 severe: Error during write
BUFF(1)='Error while writing XYZ(DA) file'
ELSEIF(IOS==608) THEN ! 608 severe: No space left on device
BUFF(1)=
1 'No space left on device to write XYZ(DA) file'
ELSE
CALL IOEMZZ(IOS,BUFF(1))
ENDIF
BUFF(2)='Close any open files or make more space.'
BUFF(3)='Press OK to retry saving file or Cancel to'
BUFF(4)='continue without saving.'
CALL ERRMSG(BUFF,4,222,IEND)
IF(IEND==0) GOTO 111
RETURN
ENDIF
ENDIF
IF(ZVGRID) THEN
112 CONTINUE
CALL OPENZZ(33,IOS)
IF(IOS/=0) THEN
CALL CLOSZZ(33) ! Close XYZ(DA) file then post query about action
C Likely IOSTAT errors are
IF(IOS==38) THEN ! 38 severe: Error during write
BUFF(1)='Error while writing XYZ(DA) file'
ELSEIF(IOS==608) THEN ! 608 severe: No space left on device
BUFF(1)=
1 'No space left on device to write XYZ(DA) file'
ELSE
CALL IOEMZZ(IOS,BUFF(1))
ENDIF
BUFF(2)='Close any open files or make more space.'
BUFF(3)='Press OK to retry saving file or Cancel to'
BUFF(4)='continue without saving.'
CALL ERRMSG(BUFF,4,222,IEND)
IF(IEND==0) GOTO 112
RETURN
ENDIF
ENDIF
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) REAL(RRINN,4)
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) (REAL(.9999,4), JJ=I1, I2)
ELSE
IF(KK.EQ.1) THEN
WRITE(NLUPHI,REC=IREC)
1 (REAL(F(L0X+(JJ-1)*NY+1),4),JJ=I1,I2)
ELSEIF(KK.EQ.2) THEN
WRITE(NLUPHI,REC=IREC)(REAL(F(L0Y+JJ),4),JJ=I1,I2)
ELSEIF(.NOT.ZVGRID.AND.KK.EQ.3) THEN
WRITE(NLUPHI,REC=IREC)
1 (REAL(F(L0zt+JJ*IFREQ+MYFIZ-1)*ZWLST,4),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) (REAL(.9999,4), 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 (REAL(F(L0VAR+IJK),4), IJK=I1, I2)
C.... Scratch file
ELSE
ISCR=ISCR+1
WRITE (NLUSCR,REC=ISCR)
1 (REAL(F(L0VAR+IJK),4), 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
C IREC=IZVG+1
IREC=IZVG+I
WRITE(NLUPHI,REC=IREC)
1 (REAL(F(L0AUX+JJ*IFREQ),4),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) (REAL(.9999,4), 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) (ARY4(IJK),IJK=1,I2-I1+1)
DO IJK=I1,I2
F(L0AUX+IJK)=ARY4(IJK-I1+1)
ENDDO
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)
1 (REAL(F(L0AUX+IJK),4),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 '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) (REAL(ZCCNF,4),KK=I1,I2)
ELSEIF(IDIR.EQ.0) THEN
WRITE(NLU,REC=IRXY) (REAL(F(L0AUX+KK),4),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 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE '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