cgxbfgr.for

C.... File name ....................GXBFGR.F.................... 090608
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!-->
c

file-name GXBFGR.htm 090608

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' 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 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 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