cgxintp.for
C SUBROUTINE GXINTP provides for the calculation of the interfacial
C pressure source terms in the momentum equations when the IPSA
C algorithm is in use, as described in PHOENICS Encyclopaedia
C article: Interfacial-Pressure Sources.
C
C Subroutine GXINTP is called from the following Groups of subroutine
C GREX3:-
C
C * Group 1, Section 1:
C in order to allocate storage and assign indices for
C identification of the continuous and dispersed phases.
C
C * Group 13, Section 16:
C in order to compute the interfacial-pressure
C momentum sources for both phases.
C
C * Group 19, Section 11:
C in order to compute Cp and B and also to make some
C once-and-for all settings for the current IZ slab.
C
C-------------------------------------------------------------------
c
SUBROUTINE GXINTP
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
COMMON/GENI/NXNY,NXM1NY,IGFIL1(7),NFM,IGFIL2(50)
COMMON/GENL/LGFIL1(14),DEBGTZ,LGFIL2(45)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB,NPAT5*1
LOGICAL LGFIL1,DEBGTZ,LGFIL2,CONPHS,SWEEP1,GRN,GTZ,LTZ,STOIPU,
1 STOIPV,STOIPW,dbgips
EQUIVALENCE (L0DPCA,L0RDA,L0CP)
SAVE JRC,JRD,JUC,JUD,JVC,JVD,JWC,JWD,JRHOC,JSIPSU,JSIPSV,JSIPSW,
1 L0CP,L0DPC,L0DR,L0SCON,L0SDIS,L0VREL,J0DPC,SWEEP1,RELIPU,
1 RELIPV,RELIPW,STOIPU,STOIPV,STOIPW,dbgips
DATA STOIPU,STOIPV,STOIPW/3*.FALSE./
NAMSUB='GXINTP'
IF(IGR.EQ.1.AND.ISC.EQ.1) THEN
IF(GRN(CPIP).OR.(.NOT.GRN(-CPIP).AND.GTZ(CPIP))) THEN
C phase 1 - continuous phase ; phase 2 - dispersed phase
CALL SUB4(JRC,R1,JRD,R2,JUC,U1,JUD,U2)
CALL SUB4(JVC,V1,JVD,V2,JWC,W1,JWD,W2)
JRHOC=LRHO1
IF(.NOT.NULLPR)
1 CALL WRIT40('phase 1 =continuous:phase 2 = dispersed ')
ELSEIF (GRN(-CPIP).OR.(.NOT.GRN(CPIP).AND.LTZ(CPIP))) THEN
C phase 2 - continuous phase ; phase 1 - dispersed phase
CALL SUB4(JRC,R2,JRD,R1,JUC,U2,JUD,U1)
CALL SUB4(JVC,V2,JVD,V1,JWC,W2,JWD,W1)
JRHOC=LRHO2
IF(.NOT.NULLPR)
1 CALL WRIT40('phase 2 =continuous:phase 1 = dispersed ')
ELSE
CALL WRIT40('CPIP setting incorrect in the Q1 file ')
CALL SET_ERR(245,'Error. CPIP setting incorrect in the Q1.',1)
C CALL EAROUT(2)
ENDIF
C remove any negative signs
IF(GRN(-ABS(CPIP))) THEN
CPIP=-ABS(CPIP)
ELSE
CPIP=ABS(CPIP)
ENDIF
C allocate storage
IF(STORE(U1)) THEN
JSIPSU=LBNAME('IPSU')
STOIPU=STORE(JSIPSU)
IF(STOIPU) RELIPU=ABS(DTFALS(JSIPSU))
ENDIF
IF(STORE(V1)) THEN
JSIPSV=LBNAME('IPSV')
STOIPV=STORE(JSIPSV)
IF(STOIPV) RELIPV=ABS(DTFALS(JSIPSV))
ENDIF
IF(STORE(W1)) THEN
JSIPSW=LBNAME('IPSW')
STOIPW=STORE(JSIPSW)
IF(STOIPW) RELIPW=ABS(DTFALS(JSIPSW))
ENDIF
CALL GXMAKE(L0CP ,NXNY, 'CP ')
MZ=ITWO(1,NZ,PARAB)
CALL GXMAKE(L0DPC ,NXNY*MZ,'DPC ')
CALL GXMAKE(L0DR, NXNY, 'DR ')
CALL GXMAKE(L0SCON,NXNY, 'SCON')
CALL GXMAKE(L0SDIS,NXNY, 'SDIS')
CALL GXMAKE(L0VREL,NXNY, 'VREL')
ELSEIF(IGR.EQ.19.AND.ISC.EQ.11) THEN
dbgips=dbcfip.and.debgtz
SWEEP1=ISWEEP.EQ.FSWEEP
c... calculate interfacial pressure coefficient
IF(GTZ(CPIP)) CALL FN1(-L0CP,CPIP)
IF(GRNDNO(1,CPIP)) CALL FN2(-L0CP,JRC,0.0,CPIPA)
IF(GRNDNO(2,CPIP)) THEN
CALL SUB2(L0RD,L0F(JRD),L0RC,L0F(JRC))
DO 1951 J=1,NXNY
F(L0CP+J)=CPIPA*(1.+F(L0RD+J))*F(L0RC+J)*F(L0RC+J)
1951 CONTINUE
ENDIF
c... calculate Pci-Pc = -cp*rhoc*(ur)**2
CALL FNVSLP(1,L0VREL,CFIPA)
CALL FN51(-L0VREL,2.0)
J0DPC=ITWO(L0DPC,L0DPC+(IZ-1)*NXNY,PARAB)
CALL FN55(-J0DPC,-L0CP,JRHOC,-L0VREL,1.0)
if(dbgips) call prn('dpc ',-J0DPC)
ELSEIF(IGR.EQ.13.AND.ISC.EQ.16) THEN
NPAT5=NPATCH(5:5)
C U-Velocity Interfacial-Pressure Momentum Sources
C ------------------------------------------------
IF(INDVAR.EQ.U1) THEN
C.... calculate sc = 0.5*(DPCP+DPCE)*Ae*(rdE-rdP)
CALL FN1(-L0SCON,0.0)
CALL FNAV(-L0DPCA,-J0DPC,'EAST')
CALL FN103(-L0DR,JRD,3)
CALL FN55(-L0SCON,-L0DPCA,LAEX,-L0DR,1.0)
C.... zero source at ix=nx when xcycle=f
IF(.NOT.XCYCLE) CALL ZERNUM(L0SCON+NXM1NY,NY)
IF(NPAT5.EQ.'L') THEN
C.... calculate sd = 0.5*(rdP+rdE)*Ae*(DPCE-DPCP) { Lahey et al }
CALL FN1(-L0SDIS,0.0)
CALL FNAV(-L0RDA,JRD,'EAST')
CALL FN103(-L0DR,-J0DPC,3)
CALL FN55(-L0SDIS,-L0RDA,LAEX,-L0DR,1.0)
IF(.NOT.XCYCLE) CALL ZERNUM(L0SDIS+NXM1NY,NY)
IF(STOIPU) THEN
CALL FN10(-L0SDIS,-L0SDIS,JSIPSU,0.0,RELIPU,(1.-RELIPU))
CALL FN0(JSIPSU,-L0SDIS)
ENDIF
ELSEIF(NPAT5.EQ.'S') THEN
IF(STOIPU) THEN
CALL FN10(-L0SCON,-L0SCON,JSIPSU,0.0,RELIPU,(1.-RELIPU))
CALL FN0(JSIPSU,-L0SCON)
ENDIF
C.... calculate sd = -sc { Stuhmiller }
CALL FN2(-L0SDIS,-L0SCON,0.0,-1.0)
ENDIF
ENDIF
C V-Velocity Interfacial-Pressure Momentum Sources
C ------------------------------------------------
IF(INDVAR.EQ.V1) THEN
C.... calculate sc = 0.5*(DPCP+DPCN)*An*(rdN-rdP)
CALL FN1(-L0SCON,0.0)
CALL FNAV(-L0DPCA,-J0DPC,'NORTH')
CALL FN103(-L0DR,JRD,1)
CALL FN55(-L0SCON,-L0DPCA,LANY,-L0DR,1.0)
C.... zero source at iy=ny
CALL ZERNY(L0SCON+NY-1)
IF(NPAT5.EQ.'L') THEN
C.... calculate sd = 0.5*(rdP+rdN)*An*(DPCN-DPCP) { Lahey et al }
CALL FN1(-L0SDIS,0.0)
CALL FNAV(-L0RDA,JRD,'NORTH')
CALL FN103(-L0DR,-J0DPC,1)
CALL FN55(-L0SDIS,-L0RDA,LANY,-L0DR,1.0)
CALL ZERNY(L0SDIS+NY-1)
IF(STOIPV) THEN
CALL FN10(-L0SDIS,-L0SDIS,JSIPSV,0.0,RELIPV,(1.-RELIPV))
CALL FN0(JSIPSV,-L0SDIS)
ENDIF
ELSEIF(NPAT5.EQ.'S') THEN
IF(STOIPV) THEN
CALL FN10(-L0SCON,-L0SCON,JSIPSV,0.0,RELIPV,(1.-RELIPV))
CALL FN0(JSIPSV,-L0SCON)
ENDIF
C.... calculate sd = -sc { Stuhmiller }
CALL FN2(-L0SDIS,-L0SCON,0.0,-1.0)
ENDIF
ENDIF
C W-Velocity interphase-pressure Momentum Sources
C ----------------------------------------
IF(INDVAR.EQ.W1) THEN
C.... calculate sc = 0.5*(DPCP+DPCH)*Ah*(rdH-rdP)
CALL FN1(-L0SCON,0.0)
IF(PARAB) THEN
CALL FN1(-L0SDIS,0.0)
ELSE
IF(SWEEP1) THEN
CALL FN0(-L0DPCA,-J0DPC)
ELSE
CALL AVENXY(L0DPCA,J0DPC,J0DPC+NXNY)
ENDIF
CALL FN10(-L0DR,HIGH(JRD),JRD,0.0,1.0,-1.0)
CALL FN55(-L0SCON,-L0DPCA,LAHZ,-L0DR,1.0)
IF(NPAT5.EQ.'L') THEN
C.... calculate sd = 0.5*(rdP+rdH)*Anh(DPCH-DPCP) { Lahey et al }
CALL FN1(-L0SDIS,0.0)
CALL SUB2(L0RD,L0F(JRD),L0RDH,L0F(HIGH(JRD)))
CALL AVENXY(L0RDA,L0RD,L0RDH)
IF(SWEEP1) THEN
CALL FN1(-L0DR,0.0)
ELSE
CALL FN10(-L0DR,-J0DPC,-(J0DPC+NXNY),0.0,-1.0,1.0)
ENDIF
CALL FN55(-L0SDIS,-L0RDA,LAHZ,-L0DR,1.0)
IF(STOIPW) THEN
CALL FN10(-L0SDIS,-L0SDIS,JSIPSW,0.0,RELIPW,(1.-RELIPW))
CALL FN0(JSIPSW,-L0SDIS)
ENDIF
ELSEIF(NPAT5.EQ.'S') THEN
IF(STOIPW) THEN
CALL FN10(-L0SCON,-L0SCON,JSIPSW,0.0,RELIPW,(1.-RELIPW))
CALL FN0(JSIPSW,-L0SCON)
ENDIF
C.... calculate sd = -sc { Stuhmiller }
CALL FN2(-L0SDIS,-L0SCON,0.0,-1.0)
ENDIF
ENDIF
ENDIF
C
CONPHS=INDVAR.EQ.JUC.OR.INDVAR.EQ.JVC.OR.INDVAR.EQ.JWC
IF(CONPHS) THEN
CALL FN0(VAL,-L0SCON)
CALL FN25(VAL,1./TINY)
ELSE
CALL FN0(VAL,-L0SDIS)
CALL FN25(VAL,1./TINY)
ENDIF
if(dbgips) then
call writ3i('indvar ',indvar,'iz ',iz,'isweep ',isweep)
call prn('sipc',-l0scon)
call prn('sipd',-l0sdis)
endif
ENDIF
END
C file-name gxlift.for 071097
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C SUBROUTINE GXLIFT provides for the calculation of the interphase
c
C lift-force terms in the momentum equations when the IPSA
C algorithm is in use, as described in PHOENICS Encyclopaedia
C article: Interfacial-Pressure Sources.
C
C Subroutine GXLIFT is called from the following Groups of subroutine
C GREX3:-
C
C * Group 1, Section 1:
C in order to allocate storage and assign indices for
C identification of the continuous and dispersed phases.
C
C * Group 13, Section 16:
C in order to compute the interphase lift momentum
C sources for both phases.
C
C * Group 19, Section 11:
C in order to compute Cl, curl(Uc), (Ud-Uc) and also to
C make some once-and-for all settings for the current
C IZ slab.
C-------------------------------------------------------------------
SUBROUTINE GXLIFT
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
COMMON/GENI/NXNY,NXM1NY,IGFIL1(7),NFM,IGFIL2(50)
COMMON/GENL/LGFIL1(14),DEBGTZ,LGFIL2(45)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
LOGICAL LGFIL1,DEBGTZ,LGFIL2,CONPHS,SWEEP1,GTZ,STOLIU,STOLIV,
1 STOLIW,dbglft
EQUIVALENCE (L0LFXA,L0LFYA,L0LFZA,L0UDA)
SAVE JCON,JRC,JRD,JUC,JUD,JVC,JVD,JWC,JWD,JRHOC,JSLISU,
1 JSLISV,JSLISW,J0LFZ,L0CL,L0LFX,L0LFY,L0LFZ,dbglft,
1 L0OMGX,L0OMGY,L0OMGZ,L0SCON,L0UCA,L0UDA,L0VRLX,L0VRLY,
1 L0VRLZ,RELLIU,RELLIV,RELLIW,SWEEP1,STOLIU,STOLIV,STOLIW
DATA STOLIU,STOLIV,STOLIW/3*.FALSE./
NAMSUB='GXLIFT'
IF(IGR.EQ.1.AND.ISC.EQ.1) THEN
IF(GTZ(CLIFTB)) THEN
C phase 2 - continuous phase ; phase 1 - dispersed phase
CALL SUB4(JRC,R2,JRD,R1,JUC,U2,JUD,U1)
CALL SUB4(JVC,V2,JVD,V1,JWC,W2,JWD,W1)
CALL SUB2(JRHOC,LRHO2,JCON,1)
IF(.NOT.NULLPR)
1 CALL WRIT40('phase 2 =continuous:phase 1 = dispersed ')
ELSE
C phase 1 - continuous phase ; phase 2 - dispersed phase
CALL SUB4(JRC,R1,JRD,R2,JUC,U1,JUD,U2)
CALL SUB4(JVC,V1,JVD,V2,JWC,W1,JWD,W2)
CALL SUB2(JRHOC,LRHO1,JCON,0)
IF(.NOT.NULLPR)
1 CALL WRIT40('phase 1 =continuous:phase 2 = dispersed ')
ENDIF
C allocate storage
IF(STORE(U1)) THEN
JSLISU=LBNAME('LISU')
STOLIU=STORE(JSLISU)
IF(STOLIU) RELLIU=ABS(DTFALS(JSLISU))
ENDIF
IF(STORE(V1)) THEN
JSLISV=LBNAME('LISV')
STOLIV=STORE(JSLISV)
IF(STOLIV) RELLIV=ABS(DTFALS(JSLISV))
ENDIF
IF(STORE(W1)) THEN
JSLISW=LBNAME('LISW')
STOLIW=STORE(JSLISW)
IF(STOLIW) RELLIW=ABS(DTFALS(JSLISW))
ENDIF
CALL GXMAKE(L0CL ,NXNY,'CL ')
CALL GXMAKE(L0UDA ,NXNY,'UDA ')
CALL GXMAKE(L0UCA ,NXNY,'UCA ')
CALL GXMAKE(L0VRLX,NXNY,'VRLX')
CALL GXMAKE(L0VRLY,NXNY,'VRLY')
CALL GXMAKE(L0VRLZ,NXNY,'VRLZ')
CALL GXMAKE(L0SCON,NXNY,'SCON')
CALL GXMAKE(L0LFX ,NXNY,'LFX ')
CALL GXMAKE(L0LFY ,NXNY,'LFY ')
MZ=ITWO(1,NZ,PARAB)
CALL GXMAKE(L0LFZ ,NXNY*MZ,'LFZ ')
C allocate storage for mean vorticity vector
CALL GXVORT(JCON,PARAB,NX,NY,NZ,L0OMGX,L0OMGY,L0OMGZ)
ELSEIF(IGR.EQ.19.AND.ISC.EQ.11) THEN
SWEEP1=ISWEEP.EQ.FSWEEP
dbglft=dbcfip.and.debgtz
c... calculate lift coefficient
IF(GRNDNO(1,CLIFT)) THEN
CALL FN2(-L0CL,JRC,0.0,CLIFTA)
ELSEIF(GRNDNO(2,CLIFT)) THEN
L0RD=L0F(JRD)
DO 1941 J=1,NXNY
F(L0CL+J)=CLIFTA*(1.-2.78*AMIN1(0.2,F(L0RD+J)))
1941 CONTINUE
ELSE
CALL FN1(-L0CL,CLIFT)
ENDIF
if(dbglft) then
call writst
call writ2i('iz ',iz,'isweep ',isweep)
call prn('cl ',-L0CL)
endif
C... calculate mean vorticity vector at cell centres
CALL GXVORT(JCON,PARAB,NX,NY,NZ,L0OMGX,L0OMGY,L0OMGZ)
if(dbglft) then
call prn('omgx',-L0OMGX)
call prn('omgy',-L0OMGY)
call prn('omgz',-L0OMGZ)
endif
C... calculate slip velocity vector at cell centres
IF(STORE(U1)) THEN
IF(NX.GT.1) THEN
CALL FNAV(-L0UDA,JUD,'WEST')
CALL FNAV(-L0UCA,JUC,'WEST')
CALL FN10(-L0VRLX,-L0UDA,-L0UCA,0.0,1.0,-1.0)
ELSE
CALL FN10(-L0VRLX,JUD,JUC,0.0,1.0,-1.0)
ENDIF
ELSE
CALL FN1(-L0VRLX,0.0)
ENDIF
IF(STORE(V1)) THEN
CALL FNAV(-L0UDA,JVD,'SOUTH')
CALL FNAV(-L0UCA,JVC,'SOUTH')
CALL FN10(-L0VRLY,-L0UDA,-L0UCA,0.0,1.0,-1.0)
ELSE
CALL FN1(-L0VRLY,0.0)
ENDIF
IF(STORE(W1)) THEN
IF(IZ.EQ.1.OR.PARAB) THEN
CALL FN10(-L0VRLZ,JWD,JWC,0.0,1.0,-1.0)
ELSEIF(IZ.EQ.NZ) THEN
CALL FN10(-L0VRLZ,LOW(JWD),LOW(JWC),0.0,1.0,-1.0)
ELSE
CALL FNAV(-L0UDA,JWD,'LOW')
CALL FNAV(-L0UCA,JWC,'LOW')
CALL FN10(-L0VRLZ,-L0UDA,-L0UCA,0.0,1.0,-1.0)
ENDIF
ELSE
CALL FN1(-L0VRLZ,0.0)
ENDIF
if(dbglft) then
call prn('vrlx',-L0VRLX)
call prn('vrly',-L0VRLY)
call prn('vrlz',-L0VRLZ)
endif
C... calculate circulation vector at cell centres
J0LFZ=ITWO(L0LFZ,L0LFZ+(IZ-1)*NXNY,PARAB)
CALL VECPRD(L0LFX,L0LFY,J0LFZ,L0VRLX,L0VRLY,L0VRLZ,
1 L0OMGX,L0OMGY,L0OMGZ)
C... calculate lift force vector per unit volume at cell centres
IF(SOLVE(U1)) THEN
CALL FN55(-L0LFX,-L0CL,JRHOC,-L0LFX,1.0)
CALL FN26(-L0LFX,JRD)
ENDIF
IF(SOLVE(V1)) THEN
CALL FN55(-L0LFY,-L0CL,JRHOC,-L0LFY,1.0)
CALL FN26(-L0LFY,JRD)
ENDIF
IF(SOLVE(W1)) THEN
CALL FN55(-J0LFZ,-L0CL,JRHOC,-J0LFZ,1.0)
CALL FN26(-J0LFZ,JRD)
ENDIF
if(dbglft) then
call writst
call writ2i('iz ',iz,'isweep ',isweep)
call prn('rhoc',JRHOC)
call prn('rd ',JRD)
call prn('lfxp',-L0LFX)
call prn('lfyp',-L0LFY)
call prn('lfzp',-J0LFZ)
endif
ELSEIF(IGR.EQ.13.AND.ISC.EQ.16) THEN
C U-Velocity Interfacial-Lift Momentum Sources
C --------------------------------------------
IF(INDVAR.EQ.U1) THEN
C.... calculate sc = 0.5*(Fx,p + Fx,E)*Ae*dx
CALL FN1(-L0SCON,0.0)
IF(NX.EQ.1) THEN
CALL FN21(-L0SCON,-L0LFX,LVOL,0.0,1.0)
ELSE
CALL FNAV(-L0LFXA,-L0LFX,'EAST')
CALL FN55(-L0SCON,-L0LFXA,LAEX,LXYDXG,1.0)
C.... zero sources at ix=nx when xcycle=f
IF(.NOT.XCYCLE) CALL ZERNUM(L0SCON+NXM1NY,NY)
ENDIF
IF(STOLIU) THEN
CALL FN10(-L0SCON,-L0SCON,JSLISU,0.0,RELLIU,(1.-RELLIU))
CALL FN0(JSLISU,-L0SCON)
ENDIF
ENDIF
C V-Velocity Interfacial-Lift Momentum Sources
C --------------------------------------------
IF(INDVAR.EQ.V1) THEN
C.... calculate sc = 0.5*(Fy,p + Fy,N)*An*dy
CALL FN1(-L0SCON,0.0)
CALL FNAV(-L0LFYA,-L0LFY,'NORTH')
CALL FN55(-L0SCON,-L0LFYA,LANY,LXYDYG,1.0)
C.... zero sources at iy=ny
CALL ZERNY(L0SCON+NY-1)
IF(STOLIV) THEN
CALL FN10(-L0SCON,-L0SCON,JSLISV,0.0,RELLIV,(1.-RELLIV))
CALL FN0(JSLISV,-L0SCON)
ENDIF
ENDIF
C W-Velocity Interfacial-Lift Momentum Sources
C --------------------------------------------
IF(INDVAR.EQ.W1) THEN
C.... calculate sc = 0.5*(Fz,p + Fz,H)*Ah*dz
CALL FN1(-L0SCON,0.0)
IF(SWEEP1.OR.PARAB) THEN
CALL FN0(-L0LFZA,-J0LFZ)
ELSE
CALL AVENXY(L0LFZA,J0LFZ,J0LFZ+NXNY)
ENDIF
IF(PARAB) THEN
CALL FN21(-L0SCON,-L0LFZA,LAHZ,0.0,DZ)
ELSE
CALL FN55(-L0SCON,-L0LFZA,LAHZ,LXYDZG,1.0)
ENDIF
IF(STOLIW) THEN
CALL FN10(-L0SCON,-L0SCON,JSLISW,0.0,RELLIW,(1.-RELLIW))
CALL FN0(JSLISW,-L0SCON)
ENDIF
ENDIF
C
CONPHS=INDVAR.EQ.JUC.OR.INDVAR.EQ.JVC.OR.INDVAR.EQ.JWC
IF(CONPHS) THEN
CALL FN0(VAL,-L0SCON)
CALL FN25(VAL,1./TINY)
ELSE
CALL FN2(VAL,-L0SCON,0.0,-1.0)
CALL FN25(VAL,1./TINY)
ENDIF
ENDIF
NAMSUB='gxlift'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C**** SUBROUTINE GXVORT provides for the calculation of the mean
C vorticity vector. It is called from group 1, section 1 and
C group 19, section 11 of GXLIFT. The coding is restricted to
C cartesian and cylindrical-polar coordinates at present, and
C for parabolic flows the z-direction velocity derivatives are
C taken as zero.
C
C The subroutine may also be used as a general utility if the
C user inserts subroutine calls in the appropriate parts of
C GROUND. When STORE(OMGX,OMGY,OMGZ) appears in the Q1 file, the
C components of the vorticity vector may be printed in the RESULT
C file or viewed via PHOTON and AUTOPLOT.
C
C.... The dummy IPH is the index, which indicates that the calculation
C is provided for the first-phase (IPH=0) or the second-phase
C (IPH=1);
C
SUBROUTINE GXVORT(IPH,PARAB,NX,NY,NZ,L0OMGX,L0OMGY,L0OMGZ)
INCLUDE 'farray'
INCLUDE 'grdear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
COMMON/GENI/NXNY,NXM1NY,IGFIL1(7),NFM,IGFIL2(50)
COMMON/GENL/LGFIL1(14),DEBGTZ,LGFIL2(45)
LOGICAL PARAB,LGFIL1,DEBGTZ,LGFIL2
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
SAVE LBOMGX,LBOMGY,LBOMGZ,L0DWDY,STOMGX,STOMGY,STOMGZ
LOGICAL STOMGX,STOMGY,STOMGZ
EQUIVALENCE (L0DUDY,L0DVDX,L0DUDX,L0DVDY,L0DUDZ,L0DWDX,
1 L0DVDZ,L0DWDY)
C
NAMSUB = 'GXVORT'
IF(IGR.EQ.1.AND.ISC.EQ.1) THEN
CALL GXMAKE(L0DWDY,NXNY,'DWDY')
CALL GXMAKE(L0OMGX,NXNY,'OMGX')
LBOMGX=LBNAME('OMGX')
STOMGX=STORE(LBOMGX)
CALL GXMAKE(L0OMGY,NXNY,'OMGY')
LBOMGY=LBNAME('OMGY')
STOMGY=STORE(LBOMGY)
CALL GXMAKE(L0OMGZ,NXNY,'OMGZ')
LBOMGZ=LBNAME('OMGZ')
STOMGZ=STORE(LBOMGZ)
ELSEIF(IGR.EQ.19.AND.ISC.EQ.11) THEN
C.... omega,x = DW/DR - DV/DZ
CALL FN1(-L0OMGX,0.0)
IF(SOLVE(W1+IPH).AND.NY.GT.1) THEN
CALL FNDWDY(-L0DWDY,W1+IPH)
CALL FN60(-L0OMGX,-L0DWDY)
ENDIF
IF((NZ.GT.1.AND..NOT.PARAB).AND.NY.GT.1) THEN
CALL FNDVDZ(-L0DVDZ,V1+IPH)
CALL FN25(-L0DVDZ,-1.0)
CALL FN60(-L0OMGX,-L0DVDZ)
ENDIF
IF(STOMGX) CALL FN0(LBOMGX,-L0OMGX)
C.... omega,y = DU/DZ - DW/DX
CALL FN1(-L0OMGY,0.0)
IF(SOLVE(U1+IPH).AND.(NZ.GT.1.AND..NOT.PARAB)) THEN
CALL FNDUDZ(-L0DUDZ,U1+IPH)
CALL FN60(-L0OMGY,-L0DUDZ)
ENDIF
IF(SOLVE(W1+IPH).AND.NX.GT.1) THEN
CALL FNDWDX(-L0DWDX,W1+IPH)
CALL FN25(-L0DWDX,-1.0)
CALL FN60(-L0OMGY,-L0DWDX)
ENDIF
IF(STOMGY) CALL FN0(LBOMGY,-L0OMGY)
C.... omega,z = DV/DX - DU/DY
CALL FN1(-L0OMGZ,0.0)
IF(NY.GT.1.AND.NX.GT.1) THEN
CALL FNDVDX(-L0DVDX,V1+IPH)
CALL FN60(-L0OMGZ,-L0DVDX)
ENDIF
IF(SOLVE(U1+IPH).AND.NY.GT.1) THEN
CALL FNDUDY(-L0DUDY,U1+IPH)
CALL FN25(-L0DUDY,-1.0)
CALL FN60(-L0OMGZ,-L0DUDY)
ENDIF
IF(STOMGZ) CALL FN0(LBOMGZ,-L0OMGZ)
ENDIF
NAMSUB = 'gxvort'
END
c