cgxbfgr.for
C.... File name ....................GXBFGR.F.................... 281021
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!-->
cfile-name GXBFGR.htm 281021
c
C**** SUBROUTINE GXBFC 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 SUBROUTINE GXBFC calls UNBFC to set velocity resolutes for the
C CCM-method. The implementation of this option differs from that
C of GXBFC only in the use of GRND2 instead of GRND1 in the
C appropriate COVAL statement.
C
C**** SUBROUTINE GXBFC is called from GREX3:-
C group 1 when BFC=T, to set up the patch-wise storage for GXBFC;
C group 11 when the patch named 'IBFC', to specify uniform flow
C in a curvilinear grid;
C group 13 when the patch name begins with 'BFC', to describe a
C uniform inflow crossing a curvilinear inlet boundary.
C
C The cartesian resolutes of the initial or inlet velocity vector
C are stored in the Q1 file in the last argument of INIT
C statements for UCRT, VCRT and WCRT for phase 1, and U2CR, V2CR
C and W2CR for phase 2.
C
C For constant density flows, the inlet density is passed from Q1
C using BFCA for phase 1, and RSG28 for phase 2. For variable
C density flows, these values would be used for all PATCHes.
C
C If either phase density is variable, and DEN1 or DEN2 is stored,
C the inlet density for individual PATCHes can be passed in the
C last argument of INIT statements for DEN1 or DEN2.
C
C Note that whereas for single phase flows, BFCA or the INIT for
C DEN1 carry the inlet density, for two phase flows it is the
C product of density and volume fraction at the inlet which is
C in question.
C
C... Input-library cases B528 (airfoil) and B524 (Mizuki Impeller)
C make use of it.
C
SUBROUTINE GXBFC
INCLUDE 'patcmn'
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
include 'bfcear'
include 'lbnamer'
DIMENSION EV(3),A(3),B(3),CC(3),D(3),IJKA(3),IJKB(3),IJKC(3),
1 IJKD(3)
COMMON/LDATA/LDUM1(3),ONEPHS,LDUM2(80)
LOGICAL LDUM1,LDUM2,ONEPHS
COMMON /IDATA/NX,NY,NZ,IDFIL1(22),FSWEEP,IDFIL2(44),NUMPAT,
1 IDFIL3(5),DEN1,DEN2,IDFIL4(42)
COMMON /GENI/IGNFL1(42),NFTOT,IGNFL2(17)
COMMON /RDATA/RDFIL1(20),RHO1,RHO2,RDFIL2(63)
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
INTEGER FSWEEP,DEN1,DEN2
LOGICAL DONE1,FLUID1,QNE,LEZ,EQZ
SAVE IJKA,IJKB,IJKC,IJKD,LBNMUC,LBNMVC,LBNMWC,l0DUMM,ISURN
C
NAMSUB = 'GXBFC '
C*******************************************************************
IF(IGR.EQ.1.AND.ISC.EQ.3) THEN
CALL GXMAKE(L0DUMM,NY*NX,'DUMM')
C----Group 1----Section 1-------------------------------------------
C... Allocate PATCH-wise storage for BFC quantities
ELSEIF(IGR.EQ.1.AND.ISC.EQ.1) THEN
DO 10 IPAT = 1,NUMPAT
NPATCH= NAMPAT(IPAT)
IF(NPATCH(1:3).EQ.'BFC' .OR. NPATCH(1:5).EQ.'SCRSF' .OR.
1 NPATCH(1:5).EQ.'SCRSO'.OR.NPATCH(1:4).EQ.'NOCP') THEN
C... Phase 1 variables
CALL MAKEPV(PVVMAS,IPAT)
CALL MAKEPV(PVURES,IPAT)
CALL MAKEPV(PVVRES,IPAT)
CALL MAKEPV(PVWRES,IPAT)
IF(.NOT.ONEPHS) THEN
C... Phase 2 variables
CALL MAKEPV(PVVMS2,IPAT)
CALL MAKEPV(PVURS2,IPAT)
CALL MAKEPV(PVVRS2,IPAT)
CALL MAKEPV(PVWRS2,IPAT)
ENDIF
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)
ENDIF
C*******************************************************************
C... Make some initial settings common to groups 11 and 13
IF(IGR.EQ.11.OR.IGR.EQ.13) THEN
IF(INDVAR==LBSURN.OR.INDVAR==LBSRN2) RETURN
DONE1 = ISWEEP .GT. FSWEEP
IF(.NOT.DONE1) THEN
C... Find pointers for 1st or 2nd phase cartesian components
IF(FLUID1(INDVAR)) THEN
LBNMUC=LBNAME('UCRT')
LBNMVC=LBNAME('VCRT')
LBNMWC=LBNAME('WCRT')
ELSE
LBNMUC=LBNAME('U2CR')
LBNMVC=LBNAME('V2CR')
LBNMWC=LBNAME('W2CR')
ENDIF
C---Set unit vector for external flow direction
CALL SUB3R(UCRT,0.0,VCRT,0.0,WCRT,0.0)
IF(STORE(LBNMUC)) CALL GETCOV(NPATCH,LBNMUC,CCC,UCRT)
IF(STORE(LBNMVC)) CALL GETCOV(NPATCH,LBNMVC,CCC,VCRT)
IF(STORE(LBNMWC)) CALL GETCOV(NPATCH,LBNMWC,CCC,WCRT)
IUVW = 3* ((INDVAR-U1)/2)
ENDIF
IF(IGR.EQ.11) THEN
C*****************************************************************
C----Group 11---Section 2---------------(GRND1 or GRND2 coefficient)
C GRND1 - standard PHOENICS.
C GRND2 - CCM method.
IF(ISC.EQ.2) THEN
C... Calculate, and set, resolute
C IUVW=0,3 or 6 depending upon whether velocity is U1/U2,
C V1/V2 or W1/W2.
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)
c ELSEIF(ISC.EQ.3) THEN
C.... GRND2 coefficient - initialise colocated velocity projection.
c CALL BFCVCM(UCRT,VCRT,WCRT)
ENDIF
ELSEIF(IGR.EQ.13) THEN
C*****************************************************************
C----Group 13---Section 13 or 14--------------(GRND1 or GRND2 value)
C GRND1 - standard PHOENICS.
C GRND2 - CCM method.
IF(ISC.EQ.13.OR.ISC.EQ.14.OR.ISC.EQ.15) THEN
C... Get pointer into PATCH-wise store IMUVW
IF(FLUID1(INDVAR)) THEN
IF(ISC.EQ.13) THEN
INDV = INDVAR
ELSE
IF(INDVAR.EQ.LBNAME('UC1')) THEN
INDV = 3
ELSEIF(INDVAR.EQ.LBNAME('VC1')) THEN
INDV = 5
ELSEIF(INDVAR.EQ.LBNAME('WC1')) THEN
INDV = 7
ELSE
INDV = INDVAR
ENDIF
ENDIF
IMUVW = PVVMAS + (INDV-1)/2
ELSE
IMUVW = PVVMS2 + (INDVAR-1)/2
ENDIF
IF(DONE1) THEN
C... If not the first sweep, copy mass-inflow rate or velocity
C resolute from IMUVW into VAL
CALL GETPTC(NPATCH,GTYPE,JXF,JXL,JYF,JYL,JZF,JZL,JTF,JTL)
JZZ=IZ
IF(IZ.EQ.(JZF-1)) IZ=IZ+1
CALL FNPAXY(VAL,IMUVW)
IZ=JZZ
ELSE
IF(INDVAR.EQ.P1.OR.INDVAR.EQ.P2) 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 or PVVMS2 for the second phase.
C For constant density flows, a single density for all BFC patches
C is carried in BFCA for phase 1, and RSG28 for phase 2. For
C variable density flows, the COVAL for DEN1 or DEN2 can carry the
C density for the current patch. Note that DEN1 or DEN2 must be
C STOREd in the Q1 file.
C Start with constant values:
IF(FLUID1(INDVAR)) THEN
CALL SUB2R(RHOP,BFCA,RHO,RHO1)
CALL SUB2(IRHO,DEN1,IPHS,1)
ELSE
CALL SUB2R(RHOP,RSG28,RHO,RHO2)
CALL SUB2(IRHO,DEN2,IPHS,2)
ENDIF
C... Now check COVALs
IF(LEZ(RHO).AND.STORE(IRHO)) THEN
CALL GETCOV(NPATCH,IRHO,CCC,RHOIN)
IF(QNE(CCC,-999.0)) RHOP=RHOIN
ENDIF
IF( (NPATCH(1:4).EQ.'SCRS'.OR.NPATCH(1:4).EQ.'NOCP')
1 .AND.INDVAR.EQ.P1) THEN
CALL L0F1(VAL,ICELL,IADD,'GXBFC ')
DO 12 IX = IXF,IXL
ICELL = ICELL + IADD
DO 15 IY = IYF,IYL
ICELL = ICELL + 1
RHOP=F(ICELL)
15 CONTINUE
12 CONTINUE
ENDIF
C... If inlet density still 0.0, stop with error message
IF(EQZ(RHOP)) THEN
CALL WRIT40('INLET DENSITY NOT SET IN GXBFC: ')
CALL WRIT1A('PATCH ',NPATCH)
CALL WRIT1I('PHASE ',IPHS)
CALL SET_ERR(507, 'Error. See result file.',2)
C CALL WAYOUT(2)
ENDIF
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('GXBFC: CALL DISALLOWED. ')
CALL WRIT40('ONLY PATCHES OF TYPES; "EAST",...,"LOW" ')
CALL WRIT40('ARE PERMITTED. ')
CALL SET_ERR(508, 'Error. See result file.',2)
C CALL WAYOUT(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,'GXBFC ')
DO 20 IX = IXF,IXL
ICELL = ICELL + IADD
DO 30 IY = IYF,IYL
ICELL = ICELL + 1
C... Get the four corners of the cell
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)
C... Construct the unit vector normal to the cell face, store in EV
CALL NORML(A,B,CC,D,EV)
C... Mass flow is density * (vector velocity).(normal unit vector)
F(ICELL) = FLIO*RHOP* (UCRT*EV(1)+VCRT*EV(2)+
1 WCRT*EV(3))
30 CONTINUE
20 CONTINUE
ELSE
IF( (ISC.EQ.13.OR.ISC.EQ.15) .AND.
1 (INDVAR.GE.U1.AND.INDVAR.LE.W2)) THEN
C... GRND1 value - calculate staggered velocity over inlet PATCH
C CALL BFCVFX(IUVW+UNIVEE,IUVW+UNIVEN,IUVW+UNIVEH,
C 1 UCRT,VCRT,WCRT,IZSTEP,VAL,L0DUMM,NX,NY)
CALL BFCVFX(IUVW+1,IUVW+2,IUVW+3,
1 UCRT,VCRT,WCRT,IZSTEP,VAL,L0DUMM,NX,NY)
ELSEIF(ISC.EQ.14) THEN
C... GRND2 value - calculate colocated velocity over inlet PATCH
CALL BFCVCM(UCRT,VCRT,WCRT)
ENDIF
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='gxbfc'
END
C***********************************************************************
c
C file name gxbfc.ftn
SUBROUTINE BFCVFX(IREFU,IREFV,IREFW,UCRT,VCRT,WCRT,IZSTEP,VAL,
1 L0DUMM,NX,NY)
C.... This subroutine calculates a velocity resolute in each cell
C for the current z-slab of the current PATCH. The velocity
C resolute calculated (U1, V1 or W1) is specified by the
C values of IREFU, IREFV and IREFW. The DUMM array is used as
C work-space.
INCLUDE 'grdbfc'
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
INTEGER VAL
LOGICAL EQZ
C
NAMSUB = 'BFCVFX'
CALL FN1(VAL,0.0)
C.... Get unit vectors for grid direction and project Cartesian
C velocity vector on to the unit vector
IF(.NOT.EQZ(UCRT)) THEN
IF(IREFU.EQ.1) THEN
CALL UNIVCS(-L0DUMM,1,3,IZSTEP)
ELSEIF(IREFU.EQ.4) THEN
CALL UNIVCS(-L0DUMM,1,1,IZSTEP)
ELSEIF(IREFU.EQ.7) THEN
CALL UNIVCS(-L0DUMM,1,5,IZSTEP)
ENDIF
CALL FN34(VAL,-L0DUMM,UCRT)
ENDIF
IF(.NOT.EQZ(VCRT)) THEN
IF(IREFV.EQ.2) THEN
CALL UNIVCS(-L0DUMM,2,3,IZSTEP)
ELSEIF(IREFV.EQ.5) THEN
CALL UNIVCS(-L0DUMM,2,1,IZSTEP)
ELSEIF(IREFV.EQ.8) THEN
CALL UNIVCS(-L0DUMM,2,5,IZSTEP)
ENDIF
CALL FN34(VAL,-L0DUMM,VCRT)
ENDIF
IF(.NOT.EQZ(WCRT)) THEN
IF(IREFW.EQ.3) THEN
CALL UNIVCS(-L0DUMM,3,3,IZSTEP)
ELSEIF(IREFW.EQ.6) THEN
CALL UNIVCS(-L0DUMM,3,1,IZSTEP)
ELSEIF(IREFW.EQ.9) THEN
CALL UNIVCS(-L0DUMM,3,5,IZSTEP)
ENDIF
CALL FN34(VAL,-L0DUMM,WCRT)
ENDIF
NAMSUB = 'bfcvfx'
END
C***********************************************************************
c
SUBROUTINE GXBFGR
C.... This subroutine exemplifies the creation, and modification with
C time, of a BFC grid. It is called from Group 1 section 2 of GREX3
C when SETBFC is true in the satellite, and from Group 19 section 2
C when MOVBFC is also true.
C Three examples are provided, namely:-
C for IBFGR = 1, a box-shaped grid;
C for IBFGR = 2, a spherical grid;
C for IBFGR = 3, a cylindrical grid, with distortion possibility.
c**********************************************************************
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
INCLUDE 'bfcear'
INCLUDE 'bfcwrk'
COMMON/GENI/IFIL1(9),NFM,IFIL2(50)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
CHARACTER*10 CASE
LOGICAL NEZ,dbnow
C
SAVE CASE,XCMIN,YCMIN,ZCMIN,XCMAX,YCMAX,ZCMAX,DXCDZC,DYCDXC,DYCDZC
SAVE TFCXMX,TFCXMN,TFCYMX,TFCYMN,TFCZMX,TFCZMN,RADINN,RADOUT
SAVE TFCRIN,TFCROU,RADIN0,XAXFAC,YAXFAC,ZPOWER,ZLENGT
SAVE TMFCRI,TMFCRO,TMFCXF,TMFCZL,TMFCYZ,TMFCZP,ECCFAC,DRDZ,GAPFAC
SAVE KWWSAV
C
NAMSUB='GXBFGR'
DBNOW=.FALSE.
if(dbnow) call writ4i('igr ',igr,'isc ',isc,
1 'isweep ',isweep,'istep ',istep)
c
IF(IGR.EQ.1.AND.ISC.EQ.2) THEN
C
c.... read special data which may have been provided by spedats in q1,
c initialising each item beforehand in case q1 does not provide it
c
if(dbnow) call writ40('reading of spedat begins')
CASE=' '
CALL GETSDC('GRIDS','CASE',CASE)
IF(CASE.EQ.'BOX') THEN
c.... starting values of box dimensions
CALL SUB3R(XCMIN,0.0,YCMIN,0.0,ZCMIN,0.0)
CALL SUB3R(XCMAX,1.0,YCMAX,1.0,ZCMAX,1.0)
CALL GETSDR('GRIDS', 'XCMIN' , XCMIN)
CALL GETSDR('GRIDS', 'XCMAX' , XCMAX)
CALL GETSDR('GRIDS', 'YCMIN' , YCMIN)
CALL GETSDR('GRIDS', 'YCMAX' , YCMAX)
CALL GETSDR('GRIDS', 'ZCMIN' , ZCMIN)
CALL GETSDR('GRIDS', 'ZCMAX' , ZCMAX)
C
c.... geometric multipliers
CALL SUB3R(DXCDZC,0.0,DYCDXC,0.0,DYCDZC,0.0)
CALL GETSDR('GRIDS', 'DXCDZC', DXCDZC)
CALL GETSDR('GRIDS', 'DYCDXC', DYCDXC)
CALL GETSDR('GRIDS', 'DYCDZC', DYCDZC)
C
c.... time multipliers
CALL SUB2R(TFCXMX,0.0,TFCXMN,0.0)
CALL SUB2R(TFCYMX,0.0,TFCYMN,0.0)
CALL SUB2R(TFCZMX,0.0,TFCZMN,0.0)
CALL GETSDR('GRIDS','TFCXMX', TFCXMX)
CALL GETSDR('GRIDS','TFCXMN', TFCXMN)
CALL GETSDR('GRIDS','TFCYMX', TFCYMX)
CALL GETSDR('GRIDS','TFCYMN', TFCYMN)
CALL GETSDR('GRIDS','TFCZMX', TFCZMX)
CALL GETSDR('GRIDS','TFCZMN', TFCZMN)
C
ELSEIF(CASE.EQ.'SPHERE') THEN
c.... starting inner and outer radii
CALL SUB2R(RADINN,0.0,RADOUT,1.0)
CALL GETSDR('GRIDS','RADINN',RADINN)
CALL GETSDR('GRIDS','RADOUT',RADOUT)
c.... time multipliers
CALL SUB2R(TFCRIN,0.0,TFCROU,0.0)
CALL GETSDR('GRIDS','TFCRIN',TFCRIN)
CALL GETSDR('GRIDS','TFCROU',TFCROU)
C
ELSEIF(CASE.EQ.'CYLINDER') THEN
CALL SUB2R(RADIN0,0.0,RADOUT,1.0)
CALL SUB2R(XAXFAC,0.0,YAXFAC,0.0)
CALL SUB2R(ZPOWER,0.0,ZLENGT,1.0)
CALL GETSDR('GRIDS','RADIN0',RADIN0)
CALL GETSDR('GRIDS','RADOUT',RADOUT)
CALL GETSDR('GRIDS','ZLENGT',ZLENGT)
CALL GETSDR('GRIDS','XAXFAC',XAXFAC)
CALL GETSDR('GRIDS','YAXFAC',YAXFAC)
CALL GETSDR('GRIDS','ZPOWER',ZPOWER)
C
CALL SUB3R(TMFCRI,0.0,TMFCRO,0.0,TMFCXF,0.0)
CALL SUB3R(TMFCZL,0.0,TMFCYZ,0.0,TMFCZP,0.0)
CALL GETSDR('GRIDS','TMFCRI',TMFCRI)
CALL GETSDR('GRIDS','TMFCRO',TMFCRO)
CALL GETSDR('GRIDS','TMFCZL',TMFCZL)
CALL GETSDR('GRIDS','TMFCXF',TMFCXF)
CALL GETSDR('GRIDS','TMFCYZ',TMFCYZ)
CALL GETSDR('GRIDS','TMFCZP',TMFCZP)
C
CALL SUB3R(ECCFAC,0.0,DRDZ,0.0,GAPFAC,0.0)
CALL GETSDR('GRIDS','ECCFAC',ECCFAC)
CALL GETSDR('GRIDS','DRDZ ',DRDZ )
CALL GETSDR('GRIDS','GAPFAC',GAPFAC)
ENDIF
KWWSAV=KWW
ENDIF
IF(MOVBFC) THEN
IF(IGR.EQ.19) THEN
IF(ISC.EQ.1) RETURN
IF(INDVAR.NE.1.AND.INDVAR.NE.9) RETURN
IF(ISWEEP.GT.2) RETURN
IF(ISTEP.GT.FSTEP.AND.ISWEEP.EQ.2) RETURN
ENDIF
ENDIF
IF((IGR.EQ.1.AND.ISC.EQ.2).OR.(IGR.EQ.19.AND.ISC.EQ.1).
1 OR.(MOVBFC.AND.IGR.EQ.19.AND.ISC.EQ.2)) THEN
IF(MOVBFC.AND.IGR.EQ.19.AND.(ISC.EQ.1.OR.ISC.EQ.2)) THEN
IF(.NOT.(ISTEP.EQ.FSTEP.AND.ISWEEP.EQ.1)) THEN
c.... store old corner coordinates for moving-bfc problems
CALL SHINM3(L0B(XCORNO),L0B(XCORNR),L0B(YCORNO),L0B(YCORNR),
1 L0B(ZCORNO),L0B(ZCORNR),(NX+1)*(NY+1)*(NZ+1))
c
C.... Now change geometry by multiplying dimensions by factors
IF(CASE.EQ.'BOX') THEN
XCMIN = XCMIN * TFCXMN
XCMAX = XCMAX * TFCXMX
YCMIN = YCMIN * TFCYMN
YCMAX = YCMAX * TFCYMX
ZCMIN = ZCMIN * TFCZMN
ZCMAX = ZCMAX * TFCZMX
ELSEIF(CASE.EQ.'SPHERE') THEN
RADINN = RADINN * TFCRIN
RADOUT = RADOUT * TFCROU
ELSEIF(CASE.EQ.'CYLINDER') THEN
RADIN0 = RADIN0 * TMFCRI
RADOUT = RADOUT * TMFCRO
if(radout.le.radin0) then
call writ2r('radout ',radout,'radin0 ',radin0)
CALL SET_ERR(247,'Error. See result file.',2)
C call earout(2)
endif
ZLENGT = ZLENGT * TMFCZL
XAXFAC = XAXFAC * TMFCXF
YAXFAC = YAXFAC * TMFCYZ
ZPOWER = ZPOWER * TMFCZP
ENDIF
ENDIF
ENDIF
C............................................ set the corner coordinates
if(dbnow) call writ40('setting corners')
C.... A simple box:
IF(CASE.EQ.'BOX') THEN
C with linear variation of x with z, depending on DXCDZC
C and of y with x and z, depending on DYCDXC and DYCDZC
C
CALL WRIT3R('XCMIN ',XCMIN,'YCMIN ',YCMIN,'ZCMIN ',
1 ZCMIN)
CALL WRIT3R('XCMAX ',XCMAX,'YCMAX ',YCMAX,'ZCMAX ',
1 ZCMAX)
CALL WRIT3R('TFCXMN ',TFCXMN,'TFCYMN ',TFCYMN,'TFCZMN ',
1 TFCZMN)
CALL WRIT3R('TFCXMX ',TFCXMX,'TFCYMX ',TFCYMX,'TFCZMX ',
1 TFCZMX)
CALL WRIT3R('DXCDZC ',DXCDZC,'DYCDXC ',DYCDXC,'DYCDZC ',
1 DYCDZC)
DO 101 K=1,NZ+1
ZC=ZCMIN + (ZCMAX - ZCMIN) * FLOAT(K-1)/FLOAT(NZ)
DO 102 I=1,NX+1
XC=XCMIN + (XCMAX - XCMIN) * FLOAT(I-1)/FLOAT(NX)
C.... Linear variation of x with z
XC=XC + DXCDZC*ZC
DO 102 J=1,NY+1
YC=YCMIN + (YCMAX - YCMIN) * FLOAT(J-1)/FLOAT(NY)
C.... Linear variation of y with x and z
YC=YC + DYCDXC*XC + DYCDZC*ZC
if(dbnow) then
call writ40
1 ('corner coordinates changed in gxbfgr ')
call writ3i(' i',i,' j',j,' k',k)
call writ3r(' xc',xc,' yc',yc,' zc',
1 zc)
endif
CALL SECRNS(XCORNR,YCORNR,ZCORNR,I,J,K,XC,YC,ZC)
102 CONTINUE
101 CONTINUE
c
ELSEIF(CASE.EQ.'SPHERE') THEN
C.... Spherical coordinates: only BFGRA and BFGRB are used
if(dbnow) call writ2r('radout ',radout,'radinn ',radinn)
TWOPI=2.0*3.14159
XANGLE=TWOPI
IF(NX.EQ.1) XANGLE=0.01
C.... J-direction is axial
XMAX=0.
YMAX=0.
ZMAX=0.
DO 201 J=1,NY+1
FACJ=FLOAT(J-1)/FLOAT(NY)
COSFAC =-COS(3.14159*FACJ)
SINFAC = SIN(3.14159*FACJ)
C.... I-direction is circumferential.
DO 202 I=1,NX+1
FACI=FLOAT(I-1)/FLOAT(NX)
FACXC=COS(FACI*XANGLE)
FACYC=SIN(FACI*XANGLE)
C.... K-direction is radial
DO 202 K=1,NZ+1
RAD=RADINN + (RADOUT-RADINN) * FLOAT(K-1)/FLOAT(NZ)
RADJ=RAD*SINFAC
XC=RADJ*FACXC
YC=RADJ*FACYC
ZC=RAD*COSFAC
xmax=amax1(xmax,xc)
ymax=amax1(ymax,yc)
zmax=amax1(ymax,zc)
if(dbnow) then
call writ40
1 ('corner coordinates changed in gxbfgr ')
call writ3i(' i',i,' j',j,' k',k)
call writ3r(' xc',xc,' yc',yc,' zc',
1 zc)
endif
CALL SECRNS(XCORNR,YCORNR,ZCORNR,I,J,K,XC,YC,ZC)
202 CONTINUE
201 CONTINUE
call writ3r(' xmax',xmax,' ymax',ymax,' zmax',
1 zmax)
c
ELSEIF(CASE.EQ.'CYLINDER') THEN
c.... cylindrical coordinates (with possibility of distortion)
if(dbnow) then
call writ40('cylindrical coordinates set')
call writ3r('radin0 ',radin0,'radout ',radout,
1 'zlength ',zlengt)
call writ3r('timfacri',tmfcri,'timfacro',tmfcro,
1 'timfaczl',tmfczl)
call writ1r('zpower ',zpower)
endif
TWOPI=2.0*3.14159
IF(NX.EQ.1) TWOPI=0.01
c.... k-direction is axial
DO 301 K=1,NZ+1
FACK=FLOAT(K-1) / FLOAT(NZ)
ZC=ZLENGT*FACK
FACZ=1.0
IF(NEZ(ZPOWER)) FACZ=(FACK+1.E-10)**ZPOWER
c.... j-direction is radial
DO 301 J=1,NY+1
FACJ=FLOAT(J-1) / FLOAT(NY)
c.... i-direction is circumferential.
DO 301 I=1,NX+1
FACI=FLOAT(I-1) / FLOAT(NX)
c.... eccentricity factor
RADIN=RADIN0 * (1.0 + ECCFAC * COS(2.0*FACI*TWOPI))
c.... variation of radius with z
RADIN=RADIN * (1.0 + DRDZ*FACK)
RADOUT=AMAX1(RADOUT,RADIN*GAPFAC)
RAD=RADIN + (RADOUT-RADIN)*FACJ
ANGLE = FACI*TWOPI-3.14159*0.5
XC=RAD * COS(ANGLE) * XAXFAC
YC=RAD * SIN(ANGLE) * (1.0 -YAXFAC*(1.0-FACZ))
CALL SECRNS(XCORNR,YCORNR,ZCORNR,I,J,K,XC,YC,ZC)
301 CONTINUE
ENDIF
C.... Calculate dependent geometrical quantities
if(dbnow) call writ40('bgeom to be called from ground ')
NODISC=.TRUE.
KWWSTO=KWW
KWW=KWWSAV
CALL BGEOM(1)
CALL BGEOM(2)
KWW=KWWSTO
C.... Put cell volume in 3d store if the variable VOLU exists
NXNY=NX*NY
L0VOLU=L0F(LBNAME('VOLU'))
IF(L0VOLU.NE.0) THEN
DO 1003 IZZ=1,NZ
L0VOL=L0B(VOLUME) + NXNY*(IZZ-1)
CALL SHINXY(L0VOLU + NFM*(IZZ-1),L0VOL)
1003 CONTINUE
ENDIF
IF(IGR.EQ.1.AND.ISC.EQ.2) THEN
C.... Save initial fields and cell-corner cordinates for examination
C via PHOTON
c.... settings and resetting needed to ensure that print-out occurs with
c 0 appended to the name
if(dbnow) call writ40('saving initial fields')
FSTSAV=FSTEP
FSTEP=0
CALL DUMPS(CSG1(1:1),CSG2(1:1),0,1,0,0)
FSTEP=FSTSAV
ENDIF
ENDIF
IF(ISWEEP.EQ.1) THEN
c.... print corner coordinates; 1 means xc, 2 means yc, 3 means zc
c call prxyzc(1)
c call prxyzc(2)
c call prxyzc(3)
c.... calls to st3d which place selected BFC geometry variables in 3D
c storage, usually for the purpose of display via PHOTON
c call st3d(volume,'VOLU')
c call st3d(aprojn,'PRJE')
c call st3d(aproje,'PRJN')
c call st3d(aprojh,'PRJH')
c call st3d(ucarte,'UCRE')
c call st3d(ucartn,'UCRN')
c call st3d(vcarte,'UCRE')
c call st3d(vcartn,'VCRN')
c call st3d(wcartn,'WCRN')
c call st3d(wcarth,'WCRH')
c call st3d(udivye,'UDYE')
c call st3d(udivyn,'UDYN')
c call st3d(udivyh,'UDYH')
c call st3d(vdivxe,'VDYE')
c call st3d(vdivxn,'VDYN')
c call st3d(vdivxh,'VDYH')
c call st3d(mfahph,'MHPH')
c call st3d(dxgpe ,'DP2E')
c call st3d(dygpn ,'DP2N')
c call st3d(mfaepe,'MEPE')
c call st3d(mfaepn,'MEPN')
c call st3d(mfanpn,'MNPN')
c call st3d(mfanph,'MNPH')
c call st3d(mfahpn,'MHPN')
c call st3d(mfahph,'MHPH')
c.... 999 is used to indicate that the variable is not one of those
c in 3D store
call writ40('storing centres')
call st3d(999 ,'XCEN')
call st3d(999 ,'YCEN')
call st3d(999 ,'ZCEN')
ENDIF
NAMSUB='gxbfgr'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
c....This subroutine sets the cartesian coordinates XC, YC, ZC, of the
c south-west-low corners of the cell with indicial coordinates I,J,K
c Its use is illustrated in sub-routine GXBFGR
c
SUBROUTINE SECRNS(XCORNR,YCORNR,ZCORNR,I,J,K,XC,YC,ZC)
INTEGER XCORNR,YCORNR,ZCORNR
CALL SECORN(XCORNR,I,J,K,XC)
CALL SECORN(YCORNR,I,J,K,YC)
CALL SECORN(ZCORNR,I,J,K,ZC)
END
c