cgxchki.htm
C element mass fraction and mole fraction
**********************************************************************
c
c
SUBROUTINE GXCHKI(ICELLA,VALOUT)
C-----------------------------------------------------------------
C PHOENICS-CHEMKIN INTERFACE
C-----------------------------------------------------------------
C Function:
C The interface uses CHEMKIN and TRANLIB routines to calculate
C thermodynamic and transport properties for gas mixtures, and
C calculates the diffusion fluxes according to 2 models of
C diffusion. The diffusion models are:
C 0) Standard PHOENICS model - ENULA = 0.0 : IMCOPT = 0 (default)
C 1) Mixture-averaged model - ENULA = GRND9: IMCOPT = 1
C Furthermore, the influence of chemical reactions is treated
C according to 2 different algorithms:
C 1) Provide source terms calculated using CHEMKIN routines
C and solve using the standard PHOENICS algorithms
C 2) Solve the equations using the PBP algorithm; but use
C a Newton-Raphson algorithm operating simultaneously on
C all species and on temperature (enthalpy) to generate
C the corrections
C The source-term options are indicated to the code by way of
C the PATCH name:
C Characters 1-5 must be CHEMK
C Character 6 must be 1 to indicate which of the
C options is to be used
C Characters 7-8 are free to be used by the user
C
C Assumptions:
C Entered only for 1st phase
C Only mixture-averaged model for transport (IMCOPT==0)
C Only solution for TEM1 is implemented (for energy equ'n.)
C Properties are calculated for every sweep
C Non-orthogonality terms have not been included in the enhanced
C diffusion model
C
C Notes: The implementation of the thermal diffusion terms in the
C species diffusion equations was modified on 21/10/03 so
C ensure species mass conservation. The previous treatment
C was discovered to be non-conservative when running library
C cases C206 to C208 inclusive.
C
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C-----------------------------------------------------------------
C
INCLUDE 'chmkin'
INCLUDE 'patcmn'
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'pltcfile'
INCLUDE 'parear'
COMMON/INDAUX/L0ISL,L0IST,L0SL,L0ST,NSTO,NSOL,L0SLRS,L0TTRS,
1 IFL(12)
C
COMMON/GENI/NXNY,IGFL1,NXNYST,IGFL4,KDUMM,IGFL5(4),NFM,NWHOLE,
1 IGSH,IGFL13(30),NFTOT,IGFL44(6),ITEM,IGFL51(5),
1 IPRPS,IGFL57(4)
COMMON/INDUN/INDFL(12),INN,INS,INE,INW,INH,INL
COMMON/IDA7/LOFLUX(150) /IDA8/LODFCO(150)
COMMON/SLDPRP/SOLPRP,PORPRP,VACPRP
COMMON/F01/M0F(600)
C
CHARACTER*256 NKEEP, NTEMP
C
PARAMETER (NCKMI0=84)
PARAMETER (VISSCA=1.E-1,CONSCA=1.E-5,DIFSCA=1.E-4)
C
INTEGER LBCKF, ICKMC0(NCKMI0), PVCKID
LOGICAL LCKERR,GRN,EQZ,VARPRS,QEQ,NEZ,SOLTEM,HRM,
1 DBGLOC,QGE,VARTEM,CONTEM,ERROR,TSTVPO,BLK,BYPASS,
2 BLOCKD,SOLID,FIXCEL,MONPRT,DENSC,CFIRST,TPSEM
CHARACTER*9 CODE
CHARACTER*40 PRECIS
CHARACTER*4 NAMD,CMFILE
C
EQUIVALENCE (PREF,CHSOC),
3 (ICKMC0(1),LBCKF),(PVCKID,LBPVSP(6)),
4 (CMFILE,CSG4)
C
C--- The following section declares the REAL CHEMKIN interface
C buffer to be optionally double- or single-precision, and
C must be configured (using the CHANGE program) to match the
C precision chosen for CHEMKIN
C
C*****double precision
DOUBLE PRECISION
1 RU,RUC,PRES,CPMIX,RHOMIX,VISMIX,
2 CONMIX,VOLUME,ATIM,RTIM,ATOL,RTOL,ABSOL,RELAT,
3 UFAC,DFAC,U0,U,RTMP,QDOT,DTTP,DTMIN,DTMAX,
4 DTEMP0,GRHO
COMMON /CKMCI0/LBCKF,KMCR,KMCI,K0MWT,K0BUF,LINK,LOUTCM,MM,KK,
1 II,L0HK,L0HKH,L0CP,L0CPH,L0H,L0HH,L0DK,L0DKH,
2 L0DJK,L0DJKH,L0VDE,L0VDN,L0VDH,L0CDK,L0TDK,L0HCO,
3 L0HVAL,L0DTF,K0L0,K0SYM,K0ELT,K0SYME,NTDIF,
4 K0ABV,K0BLO,K0TPB,K0TWP,K0IPV,K0SCR,K0VAR,
5 K0VARN,K0RES,K0RESN,K0A,K0SU,K0AP,L0SUK,L0APK,
6 K0H0K,LEVEL(2),LENTWP,ITIND,K0TD,
7 L0VHI,L0DTK,L0HSU,L0HAP,K0WDT,K0ROR,IGRND,
8 IQDOT,IENT,ICON,IMMWT,ICKOP1,IMIXF,IFUEL,IOXYG,
9 IPRO1,IPRO2,IDILU,L0APHI,L0SUHI,ICKOPT,IMCOPT,
1 ICALL1,LEVEL1,LEVEL2,NJAC,ITJAC,IRETIR,NUMDT,
2 NINIT,NRTOT,NITOT,NCTOT,NLTOT,K0EQRW,K0EQIW,K0AWT
2 /CKMCR0/ATIM,RTIM,ATOL,RTOL,UFAC,DFAC,ABSOL,RELAT,
2 RU,RUC,DTMIN,DTMAX,DTEMP0
4 /CKMCL0/VARPRS,CONTEM,SIUNIT
SAVE /CKMCI0/,/CKMCL0/,
1 DODIFF,DOSORC,DOTDIF,SOLTEM,MODDIF,VARTEM,ENTFLX,RELQDT,
2 ALQDOT,STOCMP,BYPASS,DENSC,CFIRST,TPSEM
C
COMMON /CKMIND/ NDCMI,NDCMR,NDCMC,NDCML
LOGICAL DODIFF,DOSORC,DOTDIF,DFN,WHF,MODDIF,ENTFLX,RELQDT,
1 STOCMP,SIUNIT
logical exists
character*(1024) cval, delim*1
C*****double precision
DATA PRECIS /' DOUBLE PRECISION Version is activated '/
C*****END double precision
C
dbgloc=debug.and.izstep>=izdb1 .and.izstep<=izdb2.and.
1 isweep>=iswdb1.and.isweep<=iswdb2
IXL=IABS(IXL)
NXNY=NX*NY
C
IF(IGR==1) CFIRST=.TRUE.
IF(IGR>1) THEN
IF(VARPRS) L0P1=L0F(P1)
IF(ICKOP1>1.AND.STORE(KE).AND.STORE(EP))
1 CALL SUB2(L0KE,L0F(KE),L0EP,L0F(EP))
IF(STORE(ITIND)) THEN
L0TEM=L0F(ITIND)
ELSE
L0TEM=0
VARTEM=.FALSE.
ENDIF
L0DEN=L0F(DEN1)
IF(IQDOT>0) THEN
L0QDOT=L0F(IQDOT)
ALQDOT=DTFALS(IQDOT)
RELQDT=QGE(ALQDOT,-1.0).AND.ALQDOT<0.0
ALQDOT=-DTFALS(IQDOT)
ELSE
L0QDOT=0
ENDIF
IF(IENT>0) THEN
L0ENT=L0F(IENT)
ELSE
L0ENT=0
ENDIF
IF(ICON>0) THEN
L0CON=L0F(ICON)
ELSE
L0CON=0
ENDIF
IF(IMMWT>0) THEN
L0MMWT=L0F(IMMWT)
ELSE
L0MMWT=0
ENDIF
IF(STORE(VPOR)) THEN
L0VPOR=L0F(VPOR)
TSTVPO=BLK(IZSTEP,7)
ELSE
L0VPOR=0
TSTVPO=.FALSE.
ENDIF
IF(IPRPS==0) THEN
L0PRPS=0
ELSE
L0PRPS=L0F(IPRPS)
ENDIF
ENDIF
C
IF(IGR==8.OR.IGR==9.OR.IGR==13.OR.IGR==11
1 .OR.IGR==19) THEN
IF(WHF(LBCKF)) THEN
IZOFF=(IZSTEP-1)*NXNY
ELSE
IZOFF=0
ENDIF
IF(STOCMP) THEN
DO 10002 K=1,KK
ICKMC(K0L0+K)=L0F(LBCKF+K-1)
10002 CONTINUE
ENDIF
ENDIF
GO TO (1,25,25,25,25,25,25,8,9,25,11,25,13,25,25,25,25,25,19,
1 25,25,25,23,25),IGR
25 CONTINUE
go to 999
C-----------------------------------------------------------------
C--- GROUP 1. Run title and other preliminaries
C
1 CONTINUE
C
C User may here change message transmitted to the VDU screen
IF(ISC/=1) GO TO 999
IF(.NOT.NULLPR) THEN
CALL WRITBL
CALL WRITST
CALL WRYT40('GXCHKI file is GXCHKI.HTM of: 240613 ')
CALL WRYT40(' PHOENICS-CHEMKIN Interface v1.0 ')
CALL WRYT40(PRECIS)
CALL WRIT40(' PHOENICS-CHEMKIN Interface v1.0 ')
CALL WRIT40(PRECIS)
CALL WRITST
ENDIF
C
C--- Initialise /CKMCI0/
C
DO 10000 I=1,NCKMI0
ICKMC0(I)=0
10000 CONTINUE
C
LBCKF=16
LBCKF0=NINT(CHSOB)
IF(LBCKF0/=0) LBCKF=LBCKF0
STOCMP=LBCKF>0
DTEMP0=TEMP0
C
IF(STOCMP.AND.GRNDNO(9,CHSOA)) THEN
ICKOPT=1
IF(SELREF) WRITE(LUPR1,*)
1 'SELREF has been set false in GXCHKI because CHSOA=GRND9'
SELREF=.FALSE.
ELSE
ICKOPT=0
ENDIF
C
C--- Use GETSPD to retrieve the Q1 settings which overwrite
C the default parameters
C
CALL GETSPD('CHEM','LEVEL1',2,RDUM,LEVEL1,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','LEVEL2',2,RDUM,LEVEL2,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','NJAC' ,2,RDUM,NJAC ,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','ITJAC' ,2,RDUM,ITJAC ,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','IRETIR',2,RDUM,IRETIR,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','NUMDT' ,2,RDUM,NUMDT ,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','NINIT' ,2,RDUM,NINIT ,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','CALL1' ,2,RDUM,ICALL1,.FALSE.,' ',IERR)
VARPRS=.FALSE.
CALL GETSPD('CHEM','VARPRS',3,RDUM,IDUM,VARPRS,' ',IERR)
CONTEM=.FALSE.
CALL GETSPD('CHEM','CONTEM',3,RDUM,IDUM,CONTEM,' ',IERR)
SIUNIT=.FALSE.
CALL GETSPD('CHEM','SIUNIT',3,RDUM,IDUM,SIUNIT,' ',IERR)
C option to bypass twopnt-solver termination
BYPASS=.FALSE.
CALL GETSPD('CHEM','BYPASS',3,RDUM,IDUM,BYPASS,' ',IERR)
C Provision for density conversion from SI to CGS
DENSC=.FALSE.
CALL GETSPD('CHEM','DENSC',3,RDUM,IDUM,DENSC,' ',IERR)
TPSEM=.TRUE.
CALL GETSPD('CHEM','TPSEM',3,RDUM,IDUM,TPSEM,' ',IERR)
CALL SUB4R(UFAC1,0.0,DFAC1,0.0,DTMIN1,0.0,DTMAX1,0.0)
CALL GETSPD('CHEM','UFAC' ,1,UFAC1 ,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','DFAC' ,1,DFAC1 ,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','DTMIN',1,DTMIN1,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','DTMAX',1,DTMAX1,IDUM,.FALSE.,' ',IERR)
RTOL=1.E-4
ATOL=1.E-8
ATIM=3.*ATOL
RTIM=3.*RTOL
CALL SUB4R(ATOL1,0.0,RTOL1,0.0,ATIM1,0.0,RTIM1,0.0)
CALL GETSPD('CHEM','ATOL',1,ATOL1 ,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','RTOL',1,RTOL1 ,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','ATIN',1,ATIM1,IDUM,.FALSE.,' ',IERR)
CALL GETSPD('CHEM','RTIM',1,RTIM1,IDUM,.FALSE.,' ',IERR)
IF(NEZ(ATOL1)) ATOL=ATOL1
IF(NEZ(RTOL1)) RTOL=RTOL1
IF(NEZ(ATIM1)) ATIM=ATIM1
IF(NEZ(RTIM1)) RTIM=RTIM1
C
C--- Make some debug settings
C
IF(ICKOPT>=1) THEN
CALL SUB2(LEVEL(1),LEVEL1,LEVEL(2),LEVEL2)
LEVEL(1)=MAX(LEVEL(1),LEVEL(2))
ENDIF
C
C--- Are there any source PATCHes, are there any dangerous CHSO
C PATCHes?
C
DOSORC=ICKOPT>0
LCKERR=.FALSE.
DO 10013 IP=1,NUMPAT
IF(.NOT.DOSORC.AND.NAMPAT(IP)(1:5)=='CHEMK') THEN
DOSORC=.TRUE.
READ(NAMPAT(IP)(6:6),'(I1)') ICKOP1
ELSE
LCKERR=LCKERR.OR.NAMPAT(IP)(1:4)=='CHSO'
ENDIF
10013 CONTINUE
C
C--- If CHSO has been set we have to exit
C
IF(LCKERR) THEN
CALL WRIT40(' GXCHKI - CHSO PATCHES CANNOT BE USED ')
CALL WRIT40(' WITH THE CHEMKIN INTERFACE ')
CALL SET_ERR(254, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
C
C--- Set a pointer to various variables
C
IF(STORE(ITEM)) THEN
ITIND=ITEM
ELSEIF(STORE(TEMP1)) THEN
ITIND=TEMP1
ELSE
ITIND=0
ENDIF
IMIXF=LBNAME('MIXF')
IF(.NOT.STORE(IMIXF)) IMIXF=0
SOLTEM=SOLVE(ITIND)
VARPRS=VARPRS.AND.SOLVE(P1)
IQDOT=LBNAME('QDOT')
IF(.NOT.STORE(IQDOT)) IQDOT=0
IENT =LBNAME('ENTH')
IF(.NOT.STORE(IENT)) IENT=0
ICON =LBNAME('COND')
IF(.NOT.STORE(ICON)) ICON=0
IMMWT =LBNAME('MMWT')
IF(.NOT.STORE(IMMWT)) IMMWT=0
C
C--- Initialise CHEMKIN database (name may be constructed from CMFILE)
C
NTEMP=' '
DO 10201 ITRY=0,100,100
JTRY=ITRY
IF(CMFILE/=' ') THEN
I=INDEX(CMFILE,' ')
IF(I==0) THEN
I=4
ELSE
I=I-1
ENDIF
IF(JTRY==0) THEN
NKEEP=NMFIL(55)
NMFIL(55)=CMFILE(1:I)//'ckln'
CALL BLDNAM(NMFIL(55),NTEMP,LNMTMP)
ELSE
I55=IPEXT(55,JTRY/100,0)
IC=INDEX(NMFLEX(I55),'cklink')
IF(IC<=0.OR.IC>40) THEN
IERR=1
GO TO 10202
ENDIF
NKEEP=NMFLEX(I55)
NMFLEX(I55)(IC:)=CMFILE(1:I)//'ckln'
CALL BLDNAM(NMFLEX(I55),NTEMP,LNMTMP)
ENDIF
ENDIF
IERR=-1
IF(MIMD.AND.NPROC>1.AND.JTRY==0)THEN
CVAL=NMFIL(55)
CALL PAROPENFL(90,CVAL,IERR,EXISTS)
IF(EXISTS)THEN
CALL GETSYSID(ISYSID)
DELIM='/'
IF(ISYSID<=2) DELIM=CHAR(92)
CALL COPYBINARYFILE(DELIM,CVAL,IERR,IOS,EXISTS)
IF(IERR/=0)GOTO 10201
ELSE
GOTO 10201
ENDIF
ENDIF
CALL OPENZZ(55+JTRY,IERR)
IF(IERR==0) GO TO 10202
10201 CONTINUE
10202 CONTINUE
IF(IERR/=0) THEN
CALL WRIT40(' GXCHKI - CHEMKIN LINK FILE IS ABSENT ')
CALL SET_ERR(255,
* 'Error. GXCHKI - CHEMKIN LINK FILE IS ABSENT.',1)
C CALL EAROUT(1)
ELSEIF(.NOT.NULLPR) THEN
CALL WRITBL
WRITE(LUPR1,'('' GXCHKI - CHEMKIN LINK FILE: ''/10X,A70)')
1 NTEMP
WRITE(LUPR3,'('' GXCHKI - CHEMKIN LINK FILE: ''/10X,A70)')
1 NTEMP
ENDIF
LINK=LUNIT(55)
LOUTCM=LUPR1
CALL CKLEN(LINK,LOUTCM,LCKI,LCKR,LCKC)
LCKERR=.FALSE.
C.... Allocate CHEMKIN arrays
CALL CHEM_MEM(1,LCKR,LCKI,LCKC,1,NDCMR,NDCMI,NDCMC,NDCML)
IF(LCKR>=NDCMR.OR.LCKI>=NDCMI.OR.LCKC>=NDCMC) THEN
LCKERR=.TRUE.
CALL WRIT40(' GXCHKI - Not enough room for CHEMKIN ')
CALL WRYT40(' GXCHKI - Not enough room for CHEMKIN ')
CALL WRIT40(' Currently: ')
CALL WRIT3I(' NDCMR',NDCMR,' NDCMI',NDCMI,
1 ' NDCMC',NDCMC)
CALL WRIT40(' Require: ')
CALL WRIT3I(' NDCKR', LCKR,' NDCKI', LCKI,
1 ' NDCKC', LCKC)
CALL SET_ERR(256,
* 'Error. GXCHKI - Not enough room for CHEMKIN.',1)
C CALL EAROUT(1)
ENDIF
CALL CKINIT(LCKI,LCKR,LCKC,LINK,LOUTCM,ICKMC,RCKMC,CCKMC1)
CALL PHINDX(1,LINK,LOUTCM,PRECIS,0,0)
CALL CKINDX(ICKMC,RCKMC,MM,KK,II,NFIT)
CALL CLOSFL(55)
IF(JTRY==0) THEN
NMFIL(55)=NKEEP
ELSE
I55=IPEXT(55,JTRY/100,0)
NMFLEX(I55)=NKEEP
ENDIF
IF(.NOT.NULLPR) THEN
CALL WRITBL
CALL WRIT40(' GXCHKI - CHEMKIN initialisation done ')
CALL WRITBL
CALL WRIT40(PRECIS)
ENDIF
C
C--- If MODified DIFfusion law has been selected, set UDIFF, UDIFNE
C and USOURC. BUT the modifications, including the thermal diffusion
C modifications are valid only for truly laminar flow, so we must
C deactivate the modifications if the turbulent diffusivity for the
C species are appreciable
C
IF(GRNDNO(9,ENULA)) THEN
IMCOPT=1
ELSE
IMCOPT=0
ENDIF
MODDIF=IMCOPT==1
ENTFLX=MODDIF.AND.SOLTEM
IF(NEZ(ENUT).AND.STOCMP) THEN
DO 10009 K=1,KK
LK=LBCKF+K-1
MODDIF=MODDIF.AND.QGE(PRT(LK),1.E10)
10009 CONTINUE
ENDIF
IF(MODDIF.OR.ENTFLX) CALL SUB3L(UDIFF,.TRUE.,UDIFNE,.TRUE.,
1 USOURC,.TRUE.)
IF(ICKOPT>0) USOLVE=.TRUE.
ICALL1=MAX(0,ICALL1)
IF(ICALL1==0) THEN
ICALL1=2
ELSE
ICALL1=MIN(2,ICALL1)
ENDIF
C
C--- Initialise TRANLIB if transport properties are required, ie.
C look for activation of diffusion, solution and the specific
C diffusivity option for CHEMKIN variables
C
DODIFF=.FALSE.
DO 10102 I=9,10
DODIFF=DODIFF.OR.DFN(ITIND).AND.GRNDNO(I,-PRNDTL(ITIND))
DO 10103 K=3,7,2
DODIFF=DODIFF.OR.SOLVE(K).AND.GRNDNO(I,ENUL)
10103 CONTINUE
IF(STOCMP) THEN
DO 10101 K=1,KK
LK=LBCKF+K-1
DODIFF=DODIFF.OR.SOLVE(LK).AND.DFN(LK).AND.
1 GRNDNO(I,-PRNDTL(LK))
10101 CONTINUE
ENDIF
10102 CONTINUE
MODDIF=MODDIF.AND.DODIFF
DOTDIF=DOTDIF.AND.DODIFF
IF(DODIFF) THEN
CALL SUB2(KMCR,LCKR+1,KMCI,LCKI+1)
NTEMP=' '
DO 10211 ITRY=0,100,100
IF(CMFILE/=' ') THEN
I=INDEX(CMFILE,' ')
IF(I==0) THEN
I=4
ELSE
I=I-1
ENDIF
JTRY=ITRY
IF(JTRY==0) THEN
NKEEP=NMFIL(56)
NMFIL(56)=CMFILE(1:I)//'mcln'
CALL BLDNAM(NMFIL(56),NTEMP,LNMTMP)
ELSE
I56=IPEXT(56,JTRY/100,0)
NKEEP=NMFLEX(I56)
IC=INDEX(NMFLEX(I56),'tplink')
IF(IC<=0.OR.IC>40) THEN
IERR=1
GO TO 10212
ENDIF
NKEEP=NMFLEX(I56)
NMFLEX(I56)(IC:)=CMFILE(1:I)//'mcln'
CALL BLDNAM(NMFLEX(I56),NTEMP,LNMTMP)
ENDIF
ENDIF
IERR=-1
IF(MIMD.AND.NPROC>1.AND.JTRY==0)THEN
CVAL=NMFIL(56)
CALL PAROPENFL(90,CVAL,IERR,EXISTS)
IF(EXISTS)THEN
CALL GETSYSID(ISYSID)
DELIM='/'
IF(ISYSID<=2) DELIM=CHAR(92)
CALL COPYBINARYFILE(DELIM,CVAL,IERR,IOS,EXISTS)
IF(IERR/=0)GOTO 10211
ELSE
GOTO 10211
ENDIF
ENDIF
CALL OPENZZ(56+ITRY,IERR)
IF(IERR==0) GO TO 10212
10211 CONTINUE
10212 CONTINUE
IF(IERR/=0) THEN
CALL WRIT40(' GXCHKI - TRANLIB LINK FILE IS ABSENT ')
CALL SET_ERR(257,
* 'Error. GXCHKI - TRANLIB LINK FILE IS ABSENT.',1)
C CALL EAROUT(1)
ELSEIF(.NOT.NULLPR) THEN
CALL WRITBL
WRITE(LUPR1,'('' GXCHKI - TRANLIB LINK FILE: ''/10X,A70)')
1 NTEMP
WRITE(LUPR3,'('' GXCHKI - TRANLIB LINK FILE: ''/10X,A70)')
1 NTEMP
ENDIF
LINK=LUNIT(56)
CALL MCLEN(LINK,LOUTCM,LMCI,LMCR)
IF(LMCR>=NDCMR-LCKR.OR.LMCI>=NDCMI-LCKI) THEN
NEEDR=3*(LCKR+LMCR)/2
NEEDI=3*(LCKI+LMCI)/2
CALL CHEM_MEM(2,NEEDR,NEEDI,NDCMC,NDCML,NDCMR,NDCMI,NDCMC,
1 NDCML)
ENDIF
IF(LMCR>=NDCMR-LCKR.OR.LMCI>=NDCMI-LCKI) THEN
LCKERR=.TRUE.
CALL WRIT40(' GXCHKI - Not enough room for TRANLIB ')
CALL WRYT40(' GXCHKI - Not enough room for TRANLIB ')
CALL WRIT40(' Currently: ')
CALL WRIT3I(' NDCMR',NDCMR,' NDCMI',NDCMI,
1 ' NDCMC',NDCMC)
CALL WRIT40(' Require: ')
CALL WRIT2I(' NDCKR', LCKR+LMCR,' NDCKI', LCKI+LMCI)
CALL SET_ERR(258,
* 'Error. GXCHKI - Not enough room for TRANLIB.',1)
C CALL EAROUT(1)
ENDIF
CALL MCINIT(LINK,LOUTCM,LMCI,LMCR,ICKMC(KMCI),RCKMC(KMCR))
CALL PHINDX(2,LINK,LOUTCM,PRECIS,NTDIF,K0TD)
K0TD=K0TD+LCKI
IF(.NOT.NULLPR) THEN
CALL WRITBL
CALL WRIT40(' GXCHKI - TRANLIB initialisation done ')
ENDIF
CALL CLOSFL(56)
IF(JTRY==0) THEN
NMFIL(56)=NKEEP
ELSE
I56=IPEXT(56,JTRY/100,0)
NMFLEX(I56)=NKEEP
ENDIF
ELSE
CALL SUB3(IMCOPT,0,LMCR,0,LMCI,0)
ENDIF
C
C--- Need an interface buffer and atomic weights
C
IF(IMCOPT<2) THEN
NBUF=MAX(KK*2,II)
ELSE
NBUF=KK*(1+KK)
ENDIF
KK1=KK+1
LENTWP=7*KK1+7*KK1+2
K0MWT =LCKR+LMCR
C Introduce atomic-weight storage
K0AWT= K0MWT+KK1
K0VAR =K0AWT+KK1
C
K0BUF =K0VAR+KK1
K0H0K =K0BUF+NBUF+KK1
NRTOT =K0H0K+KK1
IF(ICKOPT==1) THEN
K0ABV =NRTOT
K0BLO =K0ABV +KK1
K0TPB =K0BLO +KK1
K0TWP =K0TPB +KK1
K0SCR =K0TWP +LENTWP
K0VARN=K0SCR +8*MAX(II,KK1)
K0RES =K0VARN+KK1
K0RESN=K0RES +KK1
K0A =K0RESN+KK1
K0SU =K0A +KK1*KK1
K0AP =K0SU +KK1
NRTOT =K0AP +KK1
ENDIF
IF(NRTOT>NDCMR) THEN
NEEDR=NRTOT+1
CALL CHEM_MEM(2,NEEDR,NDCMI,NDCMC,NDCML,NDCMR,NDCMI,NDCMC,NDCML)
ENDIF
IF(NRTOT>NDCMR) THEN
CALL WRYT40(' GXCHKI - RCKMC too small for buffers ')
CALL WRIT40(' GXCHKI - RCKMC too small for buffers ')
CALL WRIT1I(' Add no.',NRTOT-NDCMR)
CALL SET_ERR(259,
* 'Error. GXCHKI - RCKMC too small for buffers.',1)
C CALL EAROUT(1)
ENDIF
C
C--- Storage for indices to species, and elemental composition
C
K0L0 =LCKI +LMCI
K0ELT =K0L0 +KK
K0IPV =K0ELT+MM*KK
K0WDT =K0IPV+KK
K0ROR =K0WDT+2*KK
K0EBUS=K0ROR+2*II
NITOT =K0EBUS+KK*II
IF(NITOT>NDCMI) THEN
NEEDI=NITOT+1
CALL CHEM_MEM(2,NDCMR,NEEDI,NDCMC,NDCML,NDCMR,NDCMI,NDCMC,NDCML)
ENDIF
IF(NITOT>NDCMI) THEN
CALL WRYT40(' GXCHKI - ICKMC too small for buffers ')
CALL WRIT40(' GXCHKI - ICKMC too small for buffers ')
CALL WRIT1I(' Add no.',NITOT-NDCMI)
CALL SET_ERR(260,
* 'Error. GXCHKI - ICKMC too small for buffers.',1)
C CALL EAROUT(1)
ENDIF
C
C--- Storage for the chemical species and element names
C
K0SYM =LCKC
K0SYME=K0SYM+KK
c K0AWT=K0SYME+MM
c NCTOT=K0AWT+MM
NCTOT=K0SYME+MM
IF(NCTOT>NDCMC) THEN
NEEDC=NCTOT+1
CALL CHEM_MEM(2,NDCMR,NDCMI,NEEDC,NDCML,NDCMR,NDCMI,NDCMC,NDCML)
ENDIF
IF(NCTOT>NDCMC) THEN
CALL WRYT40(' GXCHKI - CCKMC too small for buffers ')
CALL WRIT40(' GXCHKI - CCKMC too small for buffers ')
CALL WRIT1I(' Add no.',NCTOT-NDCMC)
CALL SET_ERR(261,
* 'Error. GXCHKI - CCKMC too small for buffers.',1)
C CALL EAROUT(1)
ENDIF
C
C--- Storage for the ACTIVE flags for TWOPNT
C
IF(ICKOPT==1) THEN
NLTOT=KK
IF(NLTOT>NDCML) THEN
NEEDL=NLTOT+1
CALL CHEM_MEM(2,NDCMR,NDCMI,NDCMC,NEEDL,NDCMR,NDCMI,NDCMC,
1 NDCML)
ENDIF
IF(NLTOT>NDCML) THEN
CALL WRYT40(' GXCHKI - LCKMC too small for buffers ')
CALL WRIT40(' GXCHKI - LCKMC too small for buffers ')
CALL WRIT1I(' Add no.',NLTOT-NDCML)
CALL SET_ERR(262,
* 'Error. GXCHKI - LCKMC too small for buffers.',1)
C CALL EAROUT(1)
ELSE
DO 10005 K=1,KK
LCKMC(K)=.FALSE.
10005 CONTINUE
ENDIF
ENDIF
C
C--- Set value of PRESS0 from CHEMKIN data; get no. of elements etc.;
C get molecular weights; species names. If SIUNIT, then convert
C pressure and RU to SI units.
C
IF(EQZ(PREF)) PREF=1.0
CALL CKRP(ICKMC,RCKMC,RU,RUC,PRES)
IF(SIUNIT) THEN
RU=RU*1.E-4
PRES=PRES*1.E-1
ENDIF
PRESS0=PRES*PREF
CALL CKWT(ICKMC,RCKMC,RCKMC(K0MWT+1))
C atomic-weight storage
CALL CKAWT(ICKMC,RCKMC,RCKMC(K0AWT+1))
CALL CKSYMS(CCKMC1,LCMOUT,CCKMC2(K0SYM+1),LCKERR)
CALL CKSYME(CCKMC1,LCMOUT,CCKMC2(K0SYME+1),LCKERR)
VARTEM=.NOT.CONTEM
RTMP=273.0
CALL CKHMS(RTMP,ICKMC,RCKMC,RCKMC(K0H0K+1))
CALL CKNCF(MM,ICKMC,RCKMC,ICKMC(K0ELT+1))
IF(LCKERR) THEN
CALL WRYT40(' GXCHKI - Failed to load species names ')
CALL WRIT40(' GXCHKI - Failed to load species names ')
CALL SET_ERR(263,
* 'Error. GXCHKI - Failed to load species names.',1)
C CALL EAROUT(1)
ENDIF
C
C--- Check that ITEM comes after species
C
IF(STORE(ITIND).AND.ITIND species ')
CALL WRIT40(' indices; please reorder the variables')
CALL SET_ERR(264,
* 'Error. GXCHKI - TEM1 index must be > species indices.',1)
C CALL EAROUT(1)
ELSE
LCKERR=.FALSE.
IF(STOCMP) THEN
K1=LBCKF
DO 10014 K=1,KK-1
LCKERR=LCKERR.OR..NOT.STORE(K1)
K1=K1+1
10014 CONTINUE
LCKERR=LCKERR.OR..NOT.STORE(LBCKF+KK-1)
IF(LCKERR) THEN
CALL WRIT40(' GXCHKI - too few variables for species ')
CALL WRIT40(' increase the no. of variables ')
CALL SET_ERR(265,
* 'Error. GXCHKI - too few variables for species.',1)
C CALL EAROUT(1)
ENDIF
IF(MOD(ISLN(LBCKF+KK-1),3)==0) THEN
ISLN(LBCKF+KK-1)=ISLN(LBCKF+KK-1)/3
F(L0ISL+LBCKF+KK-1)=-F(L0ISL+LBCKF+KK-1)
ENDIF
ENDIF
ENDIF
C
C--- Declare storage for properties
C
C--- Settings to ensure PBP flags are set so that Ap is set,
C and to ensure whole-field solver storage is deactivated
C
IF(ICKOPT==1) THEN
IF(STOCMP) THEN
K1=LBCKF
DO 10021 K=1,KK
IF(MOD(ISLN(K1),7)/=0) ISLN(K1)=ISLN(K1)*7
IF(MOD(ISLN(K1),5)==0) THEN
ISLN(K1)=ISLN(K1)/5
NWHOLE=NWHOLE-1
NWHOLE=MAX0(1,NWHOLE)
ENDIF
K1=K1+1
10021 CONTINUE
ENDIF
IF(SOLTEM) THEN
IF(MOD(ISLN(ITIND),7)/=0) ISLN(ITIND)=ISLN(ITIND)*7
IF(MOD(ISLN(ITIND),5)==0) THEN
ISLN(ITIND)=ISLN(ITIND)/5
NWHOLE=NWHOLE-1
NWHOLE=MAX0(1,NWHOLE)
ENDIF
ENDIF
ENDIF
NFTOT0=NFTOT
CALL GXMAKE(L0DTF,KK1,'DTF ')
NXYZ=NXNY
NXYZ2=NXNY
IF(WHF(LBCKF)) NXYZ=NXYZ*NZ
IF(NZ>1) NXYZ2=2*NXNY
CALL GXMAKE(L0HK ,NXYZ2*KK,'HK ')
CALL GXMAKE(L0H ,NXYZ,'H ')
CALL GXMAKE(L0CP ,NXYZ,'CP ')
IF(DODIFF) THEN
IF(IMCOPT<2) CALL GXMAKE(L0DK ,MAX(NXYZ,NXYZ2)*KK,'DK ')
IF(IMCOPT>=1) THEN
IF(IMCOPT>=2) THEN
CALL GXMAKE(L0DJK ,NXNY*KK*KK,'DJK ')
IF(NZ>1) CALL GXMAKE(L0DJKH,NXNY*KK*KK,'DJKH')
ENDIF
IF(NX>1) CALL GXMAKE(L0VDE,NXYZ*KK,'VDFE')
IF(NY>1) CALL GXMAKE(L0VDN,NXYZ*KK,'VDFN')
IF(NZ>1) CALL GXMAKE(L0VDH,NXYZ*KK,'VDFH')
ENDIF
C
C--- Determine whether or not we require the thermal diffusivities
C on the basis of are MODifying the DIFfusion law, and whether
C any of the species undergo thermal diffusion
C
DOTDIF=GRNDNO(9,ENULB).AND.NTDIF>0.AND.MODDIF
IF(DOTDIF) CALL GXMAKE(L0DTK,MAX(NXYZ,NXYZ2)*KK,'DTK ')
IF(MODDIF.OR.ENTFLX) THEN
CALL GXMAKE(L0HSU,MAX(NXYZ,NXYZ2),'HSU ')
CALL GXMAKE(L0HAP,MAX(NXYZ,NXYZ2),'HAP ')
IF(NZ>1) THEN
CALL GXMAKE(L0SUHI,NXNY,'SUHI')
CALL GXMAKE(L0APHI,NXNY,'APHI')
ENDIF
ENDIF
ENDIF
IF(DOSORC) THEN
CALL GXMAKE(L0CDK,NXYZ*KK,'CDOT')
CALL GXMAKE(L0TDK,NXYZ*KK,'TAUD')
IF(ICKOPT==0) THEN
CALL GXMAKE(L0HCO ,NXYZ,'HCO ')
CALL GXMAKE(L0HVAL,NXYZ,'HVAL')
ELSEIF(ICKOPT==1) THEN
L0SUK=L0CDK
L0APK=L0TDK
cccc USOLVE=.TRUE.
USOLVE=.TRUE.
ENDIF
ENDIF
C call writ2i('l0suk ',l0suk,'l0apk ',l0apk)
C
C--- Look for Inlet Diffusion PATCHes and create PATCH-wise
C variables
C
IF(SOLTEM) THEN
DO 10022 IP=1,NUMPAT
IF(NAMPAT(IP)(1:6)=='CHEMID') THEN
CALL MAKEPV(PVCKID,IP)
ENDIF
10022 CONTINUE
ENDIF
C
C--- Print memory requirements into the RESULT file
C
IF(.NOT.NULLPR) THEN
CALL WRITBL
CALL WRIT40(' GXCHKI - The following memory has been ')
CALL WRIT40(' taken for CHEMKIN & TRANLIB ')
CALL WRITBL
CALL WRIT2I('REAL:NO.',NRTOT,'INT.:NO.',NITOT)
CALL WRIT2I('CHR.:NO.',NCTOT,'In F:NO.',NFTOT-NFTOT0)
IF(ICKOPT==1) CALL WRIT1I('LOG.:NO.',NLTOT)
ENDIF
C
C--- Set up indices for storing WDoTs and Rates-Of-Reaction
C
IF(STOCMP) THEN
IGRND=NINT(GRND)
ICKMC(K0WDT+1)=IGRND
DO 10031 K=2,2*KK
ICKMC(K0WDT+K)=0
10031 CONTINUE
ICKMC(K0ROR+1)=IGRND
DO 10032 K=2,2*II
ICKMC(K0ROR+K)=0
10032 CONTINUE
DO 10025 K=16,NPHI
NAMD=' '
K0=0
IF(NAME(K)(4:4)=='+') THEN
NAMD=NAME(K)(1:3)//' '
I=LBNAME(NAMD)
K0=K0WDT
ELSEIF(NAME(K)(3:3)=='+') THEN
NAMD=NAME(K)(1:2)//' '
I=LBNAME(NAMD)
K0=K0WDT
ELSEIF(NAME(K)(2:2)=='+') THEN
NAMD=NAME(K)(1:1)//' '
I=LBNAME(NAMD)
K0=K0WDT
ELSEIF(NAME(K)(4:4)=='&') THEN
READ(NAME(K),'(I3)') I
K0=K0ROR
ELSEIF(NAME(K)(3:3)=='&') THEN
READ(NAME(K),'(I2)') I
K0=K0ROR
ELSEIF(NAME(K)(2:2)=='&') THEN
READ(NAME(K),'(I1)') I
K0=K0ROR
ENDIF
IF(K0>0) THEN
IF(K0==K0WDT) THEN
IF(ILBCKF+KK-1) THEN
CALL WRIT40(' GXCHKI - Failed to get WDT var. name ')
CALL SET_ERR(266,
* 'Error. GXCHKI - Failed to get WDT var. name.',1)
C CALL EAROUT(1)
ELSE
I=I-LBCKF+1
ENDIF
ELSEIF(K0==K0ROR) THEN
IF(I<1.OR.I>II) THEN
CALL WRIT40(' GXCHKI - Failed to get ROR var. name ')
CALL SET_ERR(266,
* 'Error. GXCHKI - Failed to get ROR var. name.',1)
C CALL EAROUT(1)
ENDIF
ENDIF
IF(ICKMC(K0+1)==IGRND) ICKMC(K0+1)=0
ICKMC(K0+I)=K
ENDIF
10025 CONTINUE
ENDIF
C
IF(.NOT.NULLPR) THEN
CALL WRITBL
CALL WRIT40(' The chemical species treated are: ')
WRITE(LUPR1, '('' No. Name '',
1 '' Mol.Wt. (Status)'')')
IF(STOCMP) THEN
CODE='(Solved) '
ELSE
CODE='(Virtual)'
ENDIF
ENDIF
DO 10012 K=1,KK
IF(STOCMP) THEN
IF(.NOT.SOLVE(LBCKF+K-1)) CODE='(Derived)'
ENDIF
IL=INDEX(CCKMC2(K0SYM+K),' ')
IF(IL>0.AND.IL<=5) THEN
IF(STOCMP) NAME(LBCKF+K-1)=CCKMC2(K0SYM+K)
IF(.NOT.NULLPR) WRITE(LUPR1,'(2X,I3,4X,A16,8X,F6.2,3X,A9)')
1 K,CCKMC2(K0SYM+K),RCKMC(K0MWT+K),CODE
ELSE
IF(.NOT.NULLPR.AND.STOCMP) WRITE(LUPR1,'(2X,I3,4X,A16,
1 '' ('',A4,'') '',F6.2,3X,A9)') K,CCKMC2(K0SYM+K),
2 NAME(LBCKF+K-1),RCKMC(K0MWT+K),CODE
ENDIF
10012 CONTINUE
IF(.NOT.NULLPR) THEN
CALL WRITBL
CALL WRIT1R(' PRESS0 ',PRESS0)
CALL WRITBL
IF(CONTEM) THEN
CALL WRIT40(' CONTEM=T constant temperature option ')
CALL WRITBL
CALL WRIT1R(' TEMP0 ',TEMP0)
CALL WRITBL
ENDIF
IF(.NOT.SIUNIT) THEN
CALL WRIT40(' CGS units: cm, g, erg, mol, dynes, s, K')
CALL WRITST
CALL WRITBL
ENDIF
ENDIF
C
C--- Copy the DTFALS values for the CHEMKIN species
C
IF(STOCMP) THEN
K1=LBCKF
DO 10016 K=1,KK
F(L0DTF+K)=DTFALS(K1)
K1=K1+1
10016 CONTINUE
ENDIF
IF(SOLTEM) F(L0DTF+KK+1)=DTFALS(ITIND)
C
C--- Special initialisations for the P-B-P Newton-Raphson solution
C
IF(ICKOPT==1) THEN
C
C--- Settings for the TWOPNT driver
C
IF(NJAC==0) NJAC =20
IF(ITJAC==0) ITJAC =0
IF(IRETIR==0) IRETIR=50
IF(NUMDT==0) NUMDT =50
NUMDT=MAX(0,NUMDT)
IF(NINIT<0) NINIT =50
NINIT=MAX(0,NINIT)
UFAC=UFAC1
DFAC=DFAC1
IF(EQZ(UFAC1)) UFAC=2.
IF(EQZ(DFAC1)) DFAC=4.1
UFAC=MAX(REAL(UFAC),1.01)
DFAC=MAX(REAL(DFAC),1.01)
DTMIN=DTMIN1
DTMAX=DTMAX1
IF(EQZ(DTMIN1)) DTMIN=1.E-10
IF(EQZ(DTMAX1)) DTMAX=0.1
C
C--- Determine U, the smallest increment to 1.0 which the computer
C can detect, ie. the computer's resolution
C
U = 1.0
10017 CONTINUE
U = U*0.5
U0 = 1.0 + U
c!!!! do NOT replace the next line by IF(QNE ....)) !!!!
IF(U0 /= 1.0) GO TO 10017
ABSOL = SQRT(2.0*U)
RELAT = ABSOL
IF(SOLTEM) THEN
NY0=1
RCKMC(K0ABV+1)=VARMAX(ITIND)
RCKMC(K0BLO+1)=VARMIN(ITIND)
ELSE
NY0=0
ENDIF
K1=LBCKF
ATOL=ENDIT(LBCKF)
IF(STOCMP) THEN
DO 10018 K=1,KK-1
ATOL=MIN(REAL(ATOL),ENDIT(K1))
RCKMC(K0ABV+NY0+K)=VARMAX(K1)
RCKMC(K0BLO+NY0+K)=VARMIN(K1)
K1=K1+1
10018 CONTINUE
ENDIF
C
IF(.NOT.STOCMP.AND.IMIXF>0)
1 CALL CKNU(KK,ICKMC,RCKMC,ICKMC(K0EBUS+1))
C
ENDIF
C
GO TO 999
C-----------------------------------------------------------------
C * Make changes for this group only in group 19.
C--- GROUP 7. Variables stored, solved & named
C-----------------------------------------------------------------
C
C--- GROUP 8. Terms (in differential equations) & devices
C
8 GO TO (81,82,83,84,85,86,87,88,89,810,811,812,813,814,815)
1,ISC
81 CONTINUE
C * ------------------- SECTION 1 ---------------------------
C For U1AD<=GRND--- phase 1 additional velocity. Index VELAD
GO TO 999
82 CONTINUE
C * ------------------- SECTION 2 ---------------------------
C For U2AD<=GRND--- phase 2 additional velocity. Index VELAD
GO TO 999
83 CONTINUE
C * ------------------- SECTION 3 ---------------------------
C For V1AD<=GRND--- phase 1 additional velocity. Index VELAD
GO TO 999
84 CONTINUE
C * ------------------- SECTION 4 ---------------------------
C For V2AD<=GRND--- phase 2 additional velocity. Index VELAD
GO TO 999
85 CONTINUE
C * ------------------- SECTION 5 ---------------------------
C For W1AD<=GRND--- phase 1 additional velocity. Index VELAD
GO TO 999
86 CONTINUE
C * ------------------- SECTION 6 ---------------------------
C For W2AD<=GRND--- phase 2 additional velocity. Index VELAD
GO TO 999
87 CONTINUE
C * ------------------- SECTION 7 ---- Volumetric source for gala
GO TO 999
88 CONTINUE
C * ------------------- SECTION 8 ---- Convection fluxes
GO TO 999
89 CONTINUE
C * ------------------- SECTION 9 ---- Diffusion coefficients
C--- Entered when UDIFF =.TRUE.; block-location indices are LAE
C for east, LAW for west, LAN for north, LAS for
C south, LD11 for high, and LD11 for low.
C User should provide INDVAR and NDIREC IF's as above.
C EARTH will apply the DIFCUT and GP12 modifications after the user
C has made his settings.
C
C--- Modify the diffusion coefficients to give derivatives of mole-
C fractions. See the section for diffusion neighbours.
C
IF(MODDIF.AND.INDVAR>=LBCKF.AND.
1 INDVAR<=LBCKF+KK-2) THEN
IF(NDIREC==1) THEN
CALL SUB3(KDP,L0F(LAN),KDM,L0F(LAS),L0VDI,L0VDN)
CALL SUB2(INI,INN,L0DX,L0F(LXYDY))
ELSEIF(NDIREC==3) THEN
CALL SUB3(KDP,L0F(LAE),KDM,L0F(LAW),L0VDI,L0VDE)
CALL SUB2(INI,INE,L0DX,L0F(LXYDX))
ELSEIF(NDIREC==5) THEN
CALL SUB3(KDP,L0F(LD11),KDM,0,L0VDI,L0VDH)
CALL SUB3(INI,INH,L0DX,L0F(LXYDZ),L0DXH,L0F(LXYDZG))
ENDIF
PRES=PRESS0
RTMP=TEMP0
TMPNEI=TEMP0
PRSNEI=PRESS0
K=INDVAR-LBCKF+1
K0IV=ICKMC(K0L0+K)
IVOFF =IZOFF*KK+(INDVAR-LBCKF)*NXNY
IVTOFF=IZOFF*kk+(INDVAR-LBCKF)*NXNY
K0IV=ICKMC(K0L0+K)
DO 80911 ICELL=1,NXNY
IF(VARTEM) RTMP=TEMP0+F(L0TEM+ICELL)
IF(VARPRS) PRES=PRESS0+F(L0P1+ICELL)
WTMEAN=0.0
DO 80905 KA=1,KK
WTMEAN=WTMEAN+F(ICKMC(K0L0+KA)+ICELL)/RCKMC(K0MWT+KA)
80905 CONTINUE
WTMEAN=1./WTMEAN
c write(lupr1,'('' species mass fractions:''/1p5e12.5)')
c 1 (f(ickmc(k0l0+ka)+icell),ka=1,kk)
c call writ1r(' wtmean ',wtmean)
FLUXTD=F(KDP+ICELL)
IF(KDP>0) F(KDP+ICELL)=FLUXTD*WTMEAN
IF(KDM>0) F(KDM+ICELL)=F(KDM+ICELL)*WTMEAN
IF(DOTDIF) THEN
C
C--- For thermal diffusion, we must first calculate an equivalent
C mean diffusion coefficient
C
INEI=NINT(F(INI+ICELL))
IF(NDIREC==5) THEN
INEI1=ICELL+KK*NXNY
INEI2=ICELL+L0DXH-L0DX
DXNEI=2.*F(L0DXH+ICELL)-F(L0DX+ICELL)
INEI3=ITWO(ICELL,INEI,.NOT.STORE(DEN1).AND.RHO1>=0.)
ELSE
INEI1=INEI
DXNEI =F(L0DX+INEI)
INEI3=INEI
ENDIF
RHO =F(L0DEN+ICELL)
RHONEI=F(L0DEN+INEI3)
IF(VARPRS) PRSNEI=PRESS0+F(L0P1+INEI)
IF(VARTEM) TMPNEI=TEMP0+F(L0TEM+INEI)
DTK =F(L0DTK+IVTOFF+ICELL)+TINY
DTKNEI=F(L0DTK+IVTOFF+INEI1)+TINY
DCK =F(L0DK+IVOFF+ICELL)+TINY
DCKNEI=F(L0DK+IVOFF+INEI1)+TINY
DX =F(L0DX+ICELL)
IF(HRM(INDVAR)) THEN
AVEDTK=(DX/(RHO*DCK)+DXNEI/(RHONEI*DCKNEI))/
1 (DX/(RHO*DTK*DCK)+DXNEI/(RHONEI*DTKNEI*DCKNEI))
ELSE
AVEDTK=(RHO*DTK*DCK+RHONEI*DTKNEI*DCKNEI)/
1 (RHO*DCK+RHONEI*DCKNEI)
ENDIF
GRAT=TMPNEI/RTMP
FLUXTD=-AVEDTK*FLUXTD*WTMEAN*LOG(TMPNEI/RTMP)
F(L0VDI+IVTOFF+ICELL)=FLUXTD
ENDIF
80911 CONTINUE
C
ENDIF
GO TO 999
810 CONTINUE
C * ------------------- SECTION 10 --- Convection neighbours
GO TO 999
811 CONTINUE
C * ------------------- SECTION 11 --- Diffusion neighbours
C
C--- Modify the diffusion neighbours to take account of mole-
C fraction derivatives and thermal diffusion
C
IF((ENTFLX.OR.MODDIF).AND.INDVAR>=LBCKF.AND.
1 INDVAR<=LBCKF+KK-2) THEN
IF(WHF(LBCKF)) THEN
CALL SUB2(KSU,-LBZ(L3SU,IZ),KAP,-LBZ(L3AP,IZ))
ELSE
CALL SUB2(KSU,L0F(LSU),KAP,L0F(LAP))
ENDIF
K=INDVAR-LBCKF+1
IF(INDVAR==LBCKF.AND.SOLTEM) THEN
CALL FN1(-L0HAP,0.0)
CALL FN1(-L0HSU,0.0)
IF(NZSTEP>1) THEN
CALL FN1(-L0HAP-NXNY,0.0)
CALL FN1(-L0HSU-NXNY,0.0)
ENDIF
ENDIF
IF(NDIREC==1) THEN
CALL SUB4(KDP,L0F(LAN),KDM,L0F(LAS),L0VDI,L0VDN,
1 KDNEI,L0F(LD8))
CALL SUB2(INI,INN,L0DX,L0F(LXYDY))
ELSEIF(NDIREC==3) THEN
CALL SUB4(KDP,L0F(LAE),KDM,L0F(LAW),L0VDI,L0VDE,
1 KDNEI,L0F(LD8))
CALL SUB2(INI,INE,L0DX,L0F(LXYDX))
ELSEIF(NDIREC==5) THEN
C Must use LAH here not LD11
CALL SUB4(KDP,L0F(LAH),KDM,0,L0VDI,L0VDH,KDNEI,L0F(LD7))
CALL SUB4(INI,INH,L0DX,L0F(LXYDZ),K0IHH,
1 L0HK+(K-1)*NXNY+KK*NXNY,K0HHKK,L0HK+(KK-1)*NXNY+KK*NXNY)
CALL FN1(-L0APHI,0.0)
CALL FN1(-L0SUHI,0.0)
ENDIF
PRES=PRESS0
PRESP=PRESS0
RTMP=TEMP0
RTMP=TEMP0
CALL SUB4(K0IV,ICKMC(K0L0+K),K0IH,L0HK+(K-1)*NXNY,
1 K0IHKK,L0HK+(KK-1)*NXNY,IVTOFF,IZOFF*KK+(K-1)*NXNY)
DO 81111 ICELL=1,NXNY
INEI=NINT(F(INI+ICELL))
IF(MODDIF.OR.ENTFLX) THEN
IF(VARPRS) THEN
PRES =PRESS0+F(L0P1+ICELL)
PRESP=PRESS0+F(L0P1+INEI)
ENDIF
IF(VARTEM) THEN
RTMP =TEMP0+F(L0TEM+ICELL)
TMPNEI=TEMP0+F(L0TEM+INEI)
ENDIF
ENDIF
IF(MODDIF) THEN
WTMEAN=0.0
WTMNEI=0.0
DO 81105 KA=1,KK
WTMEAN=WTMEAN+F(ICKMC(K0L0+KA)+ICELL)/RCKMC(K0MWT+KA)
WTMNEI=WTMNEI+F(ICKMC(K0L0+KA)+INEI )/RCKMC(K0MWT+KA)
81105 CONTINUE
WTMEAN=1./WTMEAN
WTMNEI=1./WTMNEI
WTMRAT=WTMEAN/WTMNEI
ELSE
WTMRAT=1.0
ENDIF
FLUX=0.0
IF(DOTDIF) THEN
FLUXTD=F(L0VDI+IVTOFF+ICELL)
IF(NDIREC==5) THEN
KAPNEI=L0APHI+ICELL
KSUNEI=L0SUHI+ICELL
ELSE
KAPNEI=KAP+INEI
KSUNEI=KSU+INEI
ENDIF
IF(FLUXTD<0.0) THEN
F(KAP+ICELL)=F(KAP+ICELL)-FLUXTD/F(K0IV+ICELL)
ELSEIF(FLUXTD>0.0) THEN
F(KAPNEI) =F(KAPNEI)+FLUXTD/F(KDNEI+ICELL)
ENDIF
F(KSU+ICELL)=F(KSU+ICELL) - FLUXTD
F(KSUNEI) =F(KSUNEI) + FLUXTD
FLUX=FLUX+FLUXTD
ENDIF
IF(MODDIF) F(KDNEI+ICELL)=F(KDNEI+ICELL)/WTMRAT
FLUX=FLUX+F(KDP+ICELL)*(F(K0IV+ICELL)-F(KDNEI+ICELL))
C--- Accumulate the enthalpy flux that results from mass-diffusion
C (note that the arrays for SU and AP terms must be zeroed on
C entry for the first CHEMKIN species)
C
IF(ENTFLX) THEN
C
C--- The enthalpy fluxes must be upwinded
C
ENTI =F(K0IH+ICELL) -F(K0IHKK+ICELL)
IF(NDIREC==5) THEN
ENTNEI=F(K0IHH+ICELL) -F(K0HHKK+ICELL)
INEIH=NXNY+ICELL
ELSE
ENTNEI=F(K0IH+INEI) -F(K0IHKK+INEI)
INEIH=INEI
ENDIF
IF(FLUX<0.0) THEN
F(L0HAP+ICELL)=F(L0HAP+ICELL)-FLUX*ENTI/RTMP
F(L0HSU+ICELL)=F(L0HSU+ICELL)-FLUX*(ENTNEI-ENTI)
ELSEIF(FLUX>0.0) THEN
F(L0HAP+INEIH)=F(L0HAP+INEIH)+FLUX*ENTNEI/TMPNEI
F(L0HSU+INEIH)=F(L0HSU+INEIH)+FLUX*(ENTNEI-ENTI)
ENDIF
ENDIF
81111 CONTINUE
ENDIF
C
GO TO 999
812 CONTINUE
C * ------------------- SECTION 12 --- Linearised sources
C
C--- Add the enthalpy flux resulting from mass-diffusion to the
C SU for the energy conservation equation, and make the
C corresponding additions to AP
C
IF((ENTFLX.OR.MODDIF).AND.SOLTEM.AND.INDVAR==ITIND) THEN
IF(WHF(LBCKF)) THEN
CALL SUB2(KAP,LBZ(L3AP,IZ),KSU,LBZ(L3SU,IZ))
ELSE
CALL SUB2(KAP,-L0F(LAP),KSU,-L0F(LSU))
ENDIF
CALL FN60(KAP,-L0HAP)
CALL FN60(KSU,-L0HSU)
ENDIF
C
GO TO 999
813 CONTINUE
C * ------------------- SECTION 13 --- Correction coefficients
GO TO 999
814 CONTINUE
C * ------------------- SECTION 14 --- User's own solver
C
C--- The CHEMKIN variables are solved on a point-by-point (PBP)
C basis using the Sandia TWOPNT solver (Newton-Raphson iteration)
C If this is not a CHEMKIN variable (assumed to be the default)
C then we want to go into the PHOENICS solver so we set USER=.F.
C For CHEMKIN variables we must initialise LD7, into which the
C PHOENICS solver would put the corrections, to zero.
C
USER=.FALSE.
C
C--- Set up indices for storage of Rates-Of-Reaction and WDoTs
C
IF(STOCMP) THEN
IF(ICKMC(K0WDT+1)/=IGRND) THEN
DO 81401 K=1,KK
I=ICKMC(K0WDT+K)
IF(I>0) ICKMC(K0WDT+KK+K)=L0F(I)
81401 CONTINUE
ENDIF
IF(ICKMC(K0ROR+1)/=IGRND) THEN
DO 81402 K=1,II
I=ICKMC(K0ROR+K)
IF(I>0) ICKMC(K0ROR+II+K)=L0F(I)
81402 CONTINUE
ENDIF
ENDIF
C
IF((SOLTEM.AND.INDVAR==ITIND.OR.
1 .NOT.SOLTEM.AND.INDVAR==LBCKF+KK-2)
2 .AND.ICKOPT==1) THEN
C
C--- Initialise ENUFSW for later setting on the basis of residual
C on entry to PHOTWO, set USER to cut out the PHOENICS solver,
C and ensure that the dummy array (LD7) that usually contains
C the corrections is set to zero
C
ENUFSW=.TRUE.
DO 81411 K=1,KK
LK=LBCKF+K-1
IF(ISL(LK)/=0) F(L0SLRS+IABS(ISL(LK)))=0.0
81411 CONTINUE
IF(SOLTEM) F(L0SLRS+ISL(ITIND))=0.0
IF(IZSTEP==1) THEN
DO 81412 K=1,KK
LK=LBCKF+K-1
IF(ISL(LK)/=0) F(L0TTRS+IABS(ISL(LK)))=0.0
81412 CONTINUE
IF(SOLTEM) F(L0TTRS+ISL(ITIND))=0.0
ENDIF
USER=.TRUE.
CALL FN1(LD7,0.0)
C
C--- Final CHEMKIN variable so now loop over cells solving on a point-
C by-point basis: first add in the extra diffusion sources for
C enthalpy and set up indices dependent on whether or not the
C temperature variable is being solved
C
CALL SUB3(L0SU,L0F(LSU),L0AP,L0F(LAP),L0VOL,L0F(VOL))
IF(INDVAR==ITIND) THEN
NATJ=KK
NY0=1
ELSE
NATJ=KK-1
NY0=0
ENDIF
I0=L0SUK
ISA=L0APK-L0SUK
PRES=PRESS0
RTMP=TEMP0
DTTP=DTFALS(LBCKF)/(UFAC-1.)
C
C--- Now loop over cells calling the solver-driver PHOTWO to generate
C the solutions
C
DO 81441 ICELL=1,NXNY
C
C--- Before doing anything else look to see whether the cell is
C blocked. If blocked we do nothing for the chemical species
C but still do a PBP solution for the TEM1 variable
C
CALL SUB2L(BLOCKD,.FALSE.,SOLID,.FALSE.)
IF(IPRPS>0) THEN
SOLID =F(L0PRPS+ICELL)>=SOLPRP
BLOCKD=F(L0PRPS+ICELL)>=PORPRP
ENDIF
C IF(TSTVPO) BLOCKD=BLOCKD.OR.F(L0VPOR+ICELL)<=1.E-10
IF(SOLID.AND.INDVAR==ITIND) THEN
C
C--- If the cell is a solid then we solve only for temperature, and
C if temperature is solved we skip over all solution
C
RESID=ABS(F(L0SU+ICELL))
DELVAR=F(L0SU+ICELL)/F(L0AP+ICELL)
F(L0TEM+ICELL)=F(L0TEM+ICELL)+DELVAR
GO TO 81427
ELSEIF(BLOCKD.OR.SOLID) THEN
C
C--- If the cell is completely blocked, do nothing
C
GO TO 81439
ENDIF
C--- If FIXVAL cell then (a) fix temperature and species; and
C (b) skip to end of loop bypassing call to photwo
FIXCEL=.FALSE.
IF(NY0>0) THEN
IF(F(L0AP+ICELL)>=FIXVAL) THEN
FIXCEL=.TRUE.
DELVAR=F(L0SU+ICELL)/F(L0AP+ICELL)
F(L0TEM+ICELL)=F(L0TEM+ICELL)+DELVAR
ENDIF
ENDIF
DO 81499 K=1,NATJ-1
IF(F(I0+ISA+K)>=FIXVAL) THEN
FIXCEL=.TRUE.
DELVAR=F(I0+K)/F(I0+ISA+K)
F(ICKMC(K0L0+K)+ICELL)=F(ICKMC(K0L0+K)+ICELL)+DELVAR
ENDIF
81499 CONTINUE
IF(FIXCEL) GO TO 81439
C---
C
C--- First copy the sources, Su, the coefficients, Ap, and the CHEMKIN
C variables solved into buffer stores, set the pressure, if local
C pressure is to be used in CHEMKIN, and get the cell volume for
C calculating the reaction source and internal false time-step
C terms
C
IF(VARPRS) PRES=PRESS0+F(L0P1+ICELL)
IF(VARTEM) RTMP=F(L0TEM+ICELL)
DO 81405 K=1,KK
RCKMC(K0SU+NY0+K) =F(I0+K)
RCKMC(K0AP+NY0+K) =F(I0+ISA+K)
RCKMC(K0VAR+NY0+K) =F(ICKMC(K0L0+K)+ICELL)
RCKMC(K0VARN+NY0+K)=F(ICKMC(K0L0+K)+ICELL)
81405 CONTINUE
IF(INDVAR==ITIND) THEN
RCKMC(K0VAR+1)=F(L0TEM+ICELL)
RCKMC(K0VARN+1)=F(L0TEM+ICELL)
RCKMC(K0SU+1)=F(L0SU+ICELL)
RCKMC(K0AP+1)=F(L0AP+ICELL)
ELSE
RCKMC(K0SU+NATJ)=F(L0SU+ICELL)
RCKMC(K0AP+NATJ)=F(L0AP+ICELL)
ENDIF
VOLUME=F(L0VOL+ICELL)
if(dbgloc) call writ1i(' icell',icell)
C
C--- Call the PHOenics TWOpnt driver routine PHOTWO to generate the
C solution for the current cell
C
CALL PHOTWO(LOUTCM,KK,NATJ,1,1,SOLTEM,LENTWP,LEVEL,
1 LCKMC,ICKMC,RCKMC,VOLUME,RCKMC(K0MWT+1),
2 ATIM,RTIM,ATOL,RTOL,ABSOL,RELAT,RCKMC(K0ABV+1),
3 RCKMC(K0BLO+1),RCKMC(K0TPB+1),RCKMC(K0TWP+1),
4 ICKMC(K0IPV+1),RCKMC(K0SCR+1),RCKMC(K0VAR+1),
5 RCKMC(K0VARN+1),RCKMC(K0RES+1),RCKMC(K0RESN+1),
6 RCKMC(K0A+1),RCKMC(K0SU+1),RCKMC(K0AP+1),
7 RCKMC(K0H0K+1),UFAC,DFAC,PRES,RTMP,NINIT,NJAC,
8 ITJAC,IRETIR,NUMDT,DTTP,DTMIN,DTMAX,
9 QDOT,ICALL1,DTEMP0,DBGLOC.OR.FLAG,ERROR,TPSEM)
C
C--- First check for errors reported by PHOTWO
C
IF(ERROR) THEN
IF(BYPASS) THEN
IF(TPSEM) THEN
CALL WRITBL
CALL WRITST
CALL WRIT40(' GXCHKI - PHOTWO failed to converge at ')
CALL WRIT4I(' ISWEEP',ISWEEP,' ITHYD',ITHYD,
1 ' IZSTEP',IZSTEP,' ICELL',ICELL)
ENDIF
ELSE
CALL WRITBL
CALL WRITST
CALL WRIT40(' GXCHKI - PHOTWO reports a fatal error ')
CALL WRIT1I('INDVAR ',INDVAR)
CALL WRIT4I(' ISWEEP',ISWEEP,' ITHYD',ITHYD,
1 ' IZSTEP',IZSTEP,' ICELL',ICELL)
CALL WRIT40(' Solution for this cell is: ')
WRITE(LUPR1,'('' Spec. Value VARMIN'',
1 '' VARMAX Ap Res'')')
IF(NY0>0) WRITE(LUPR1,'(1X,A4,3X,1P5E12.3)') 'TEM1',
1 RCKMC(K0VAR+1), VARMIN(ITIND), VARMAX(ITIND),
2 RCKMC(K0AP+1), RCKMC(K0SCR+2*KK+1)
DO 81491 K=NY0+1,NATJ+1
KIV=LBCKF+K-NY0-1
IF(SOLVE(KIV)) THEN
WRITE(LUPR1,'(1X,A4,3X,1P5E12.3)') NAME(KIV),
1 RCKMC(K0VAR+K), VARMIN(KIV), VARMAX(KIV),
2 RCKMC(K0AP+K), RCKMC(K0SCR+2*KK+K)
ELSE
WRITE(LUPR1,'(1X,A4,3X,1P5E12.3)') NAME(KIV),
1 RCKMC(K0VAR+K), VARMIN(KIV), VARMAX(KIV)
ENDIF
81491 CONTINUE
CALL WRITBL
CALL WRITST
CALL WRITBL
CALL SET_ERR(267, 'Error. See result file.',1)
C CALL EAROUT(1)
ENDIF
ENDIF
C
C--- Copy the new solution back into the PHOENICS variable stores,
C calculating the mass-fraction of the last species form the sum
C of the other species mass-fractions
C
SUMY=0.0
DO 81431 K=1,KK-1
SUMY=SUMY+RCKMC(K0VAR+NY0+K)
RCKMC(K0BUF+K)=RCKMC(K0VAR+NY0+K)
F(ICKMC(K0L0+K)+ICELL)=RCKMC(K0VAR+NY0+K)
81431 CONTINUE
RCKMC(K0BUF+KK)=1.-SUMY
F(ICKMC(K0L0+KK)+ICELL)=RCKMC(K0BUF+KK)
IF(SOLTEM) F(L0TEM+ICELL)=RCKMC(K0VAR+1)
C
C--- Load the print-out arrays for Rates-Of-Reaction, WDoTs &
C Note that the first KK elements of the scratch array (K0SCR+...)
C contain the rates-of-production (WDoTs)
C
IF(L0QDOT>0) F(L0QDOT+ICELL)=QDOT
IF(ICKMC(K0WDT+1)/=IGRND) THEN
DO 81421 K=1,KK
I=ICKMC(K0WDT+KK+K)
IF(I>0) F(I+ICELL)=RCKMC(K0SCR+K)*RCKMC(K0MWT+K)
81421 CONTINUE
ENDIF
IF(ICKMC(K0ROR+1)/=IGRND) THEN
CALL CKQYP(PRES,RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,
1 RCKMC(K0BUF+KK+1))
DO 81422 K=1,II
I=ICKMC(K0ROR+II+K)
IF(I>0) F(I+ICELL)=RCKMC(K0BUF+KK+K)
81422 CONTINUE
ENDIF
C
C--- Set ENUFSW on the basis of residual for the current cell on entry
C to the solution procedure, and on the basis of the change in the
C variables solved using PHOTWO for this iteration. If both are
C small enough when compared with the RESREF and ENDIT variables
C signal that ENUFSW need not be .FALSE.
C
DO 81425 K=NY0+1,NATJ
KIV=LBCKF+K-NY0-1
IF(RESREF(KIV)<0) ENUFSW=ENUFSW.AND.
1 ABS(RCKMC(K0SCR+2*KK+K))/(TINY+ABS(RESREF(KIV)))<1.0
DELVAR=RCKMC(K0VAR+K)-RCKMC(K0VARN+K)
ENUFSW=ENUFSW.AND.ABS(DELVAR)=LBCKF.AND.INDVAR<=LBCKF+KK-2.AND.
1 ICKOPT==1) THEN
C
C--- If this is a CHEMKIN variable, but not the last, then just
C save the SU and AP, zero the correction store, and set USER
C to indicate that the PHOENICS solver should do nothing
C
USER=.TRUE.
CALL FN1(LD7,0.0)
I=L0SUK+INDVAR-LBCKF+1
ISA=L0APK-L0SUK
L0SU=L0F(LSU)
L0AP=L0F(LAP)
DO 81451 ICELL=1,NXNY
C
CALL SUB2L(BLOCKD,.FALSE.,SOLID,.FALSE.)
IF(IPRPS>0) THEN
SOLID =F(L0PRPS+ICELL)>=SOLPRP
BLOCKD=F(L0PRPS+ICELL)>=PORPRP
ENDIF
C IF(TSTVPO) BLOCKD=BLOCKD.OR.F(L0VPOR+ICELL)<=1.E-10
IF(.NOT.(BLOCKD.OR.SOLID)) THEN
C
C--- If the cell is completely blocked, do nothing
C
F(I)=F(L0SU+ICELL)
F(I+ISA)=F(L0AP+ICELL)
ELSE
F(I)=0.0
F(I+ISA)=1.E10
ENDIF
I=I+KK
81451 CONTINUE
ENDIF
GO TO 999
815 CONTINUE
C * ------------------- SECTION 15 --- Change solution
GO TO 999
C
C * See the equivalent section in GREX for the indices to be
C used in sections 7 - 15
C
C * Make all other group-8 changes in GROUP 19.
C-----------------------------------------------------------------
C
C--- GROUP 9. Properties of the medium (or media)
C
C The sections in this group are arranged sequentially in their
C order of calling from EARTH. Thus, as can be seen from below,
C the temperature sections (10 and 11) precede the density
C sections (1 and 3); so, density formulae can refer to
C temperature stores already set.
9 CONTINUE
C-----------------------------------------------------------------
C
C--- The next section is done for all sections
C
ICELL=ICELLA
PRES=PRESS0
RTMP=TEMP0
IF(VARPRS) PRES=PRES+F(L0P1+ICELL)
IF(VARTEM) RTMP=RTMP+F(L0TEM+ICELL)
SUMY=0.0
SUMW=0.0
C IF(ICKOP1/=2.OR.IMIXF==0) THEN
IF(STOCMP) THEN
DO 9001 K=1,KK-1
RCKMC(K0BUF+K)=F(ICKMC(K0L0+K)+ICELL)
SUMW=SUMW+RCKMC(K0BUF+K)/RCKMC(K0MWT+K)
SUMY=SUMY+RCKMC(K0BUF+K)
9001 CONTINUE
ENDIF
C
RCKMC(K0BUF+KK)=1.-SUMY
SUMW=SUMW+RCKMC(K0BUF+KK)/RCKMC(K0MWT+KK)
IF(IMMWT>0) F(L0MMWT+ICELL)=1./SUMW
IF(IGSH>0) IZOFF=IZOFF+NXNY
C--- Check whether we are in a solid (special actions taken)
C probably done in GXPRPS
C
IF(ISC==14) THEN
C * ------------------- SECTION 14 ---------------------------
C For SOLVE(TEMP1)-------------------------- phase-1 specific heat
C
C....
CALL SUB2(L0HZ,L0H+IZOFF,L0CPZ,L0CP+IZOFF)
CALL CKHMS(RTMP,ICKMC,RCKMC,RCKMC(K0BUF+KK+1))
HMIX=0.0
IF(IGSH>0) THEN
CALL SUB2(K0L,L0HK,K0,L0HK+KK*NXNY)
KADD=0
DO 90001 K=1,KK
C CALL FN0(-K0L-KADD,-K0-KADD)
F(K0L+KADD+ICELL)=F(K0+KADD+ICELL)
KADD=KADD+NXNY
90001 CONTINUE
ELSE
CALL SUB2(K0,L0HK,K0H,L0HK+KK*NXNY)
ENDIF
C
C--- Calculate the sensible enthalpies, store them and calculate
C the mixture enthalpy
C
DO 9002 K=1,KK
HSENS=RCKMC(K0BUF+KK+K)-RCKMC(K0H0K+K)
F(K0+ICELL)=HSENS
HMIX=HMIX+RCKMC(K0BUF+K)*HSENS
K0=K0+NXNY
9002 CONTINUE
c VALOUT=HMIX/RTMP
VALOUT=HMIX/f(l0tem+icell)
F(L0HZ+ICELL)=HMIX
IF(L0ENT>0) THEN
HMIXT=0.0
DO 9021 K=1,KK
HMIXT=HMIXT+RCKMC(K0BUF+K)*RCKMC(K0BUF+KK+K)
9021 CONTINUE
F(L0ENT+ICELL)=F(L0DEN+ICELL)*HMIXT
ENDIF
CALL CKCPBS(RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,CPMIX)
F(L0CPZ+ICELL)=CPMIX
IF(IGSH==0.AND.NZ>1) THEN
CALL SUB3(K0,L0HK,K0H,L0HK+KK*NXNY,KADD,0)
DO 90011 K=1,KK
C CALL FN0(-K0H-KADD,-K0-KADD)
F(K0H+KADD+ICELL)=F(K0+KADD+ICELL)
KADD=KADD+NXNY
90011 CONTINUE
ENDIF
C
if(dbgloc.and.dbt) then
call writ3i(' icell',icell,' isweep',isweep,
1 ' izstep',izstep)
call writ2r(' hmix',real(hmix),' cpmix',real(cpmix))
write(lupr1,'('' individual enthalpies are:''/1p5e12.5)')
1 (rckmc(k0buf+kk+k),k=1,kk)
endif
C
ELSEIF(ISC==1) THEN
C * ------------------- SECTION 1 ---------------------------
C For RHO1<=GRND--- density for phase 1 Index DEN1
C
CALL CKRHOY(PRES,RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,RHOMIX)
VALOUT=RHOMIX
C
if(dbgloc.and.dbrho) then
call writ3i(' icell',icell,' isweep',isweep,
1 ' izstep',izstep)
call writ1r(' rho',valout)
endif
C
ELSEIF(ISC==6) THEN
C * ------------------- SECTION 6 ---------------------------
C For ENUL<=GRND--- reference laminar kinematic viscosity
C Index VISL
C
VALOUT=0.0
IF(DODIFF) THEN
VISMIX=0.0
CALL CKYTX(RCKMC(K0BUF+1),ICKMC,RCKMC,RCKMC(K0BUF+1))
CALL MCAVIS(RTMP,RCKMC(K0BUF+1),RCKMC(KMCR),VISMIX)
VALOUT=VISMIX/F(L0DEN+ICELL)
IF(SIUNIT) VALOUT=VALOUT*VISSCA
C
if(dbgloc.and.dbemu) then
call writ3i(' icell',icell,' isweep',isweep,
1 ' izstep',izstep)
call writ1r(' enu',valout)
endif
ENDIF
C
ELSEIF(ISC==7) THEN
C * ------------------- SECTION 7 ---------------------------
C For PRNDTL( )<=GRND--- laminar PRANDTL nos., or diffusivity
C Index LAMPR
IF(DODIFF) THEN
CALL CKYTX(RCKMC(K0BUF+1),ICKMC,RCKMC,RCKMC(K0BUF+1))
IF(INDVAR==ITIND) THEN
CALL MCACON(RTMP,RCKMC(K0BUF+1),RCKMC(KMCR),CONMIX)
VALOUT=CONMIX
IF(SIUNIT) VALOUT=VALOUT*CONSCA
IF(L0CON>0) F(L0CON+ICELL)=VALOUT
C
if(dbgloc.and.dbgam) then
call writ3i(' icell',icell,' isweep',isweep,
1 ' izstep',izstep)
call writ1r(' cond',valout)
endif
C
ELSEIF(IMCOPT<2) THEN
IF(INDVAR==LBCKF) THEN
if(dbgloc.and.dbgam) then
call writ3i(' icell',icell,' isweep',isweep,
1 ' izstep',izstep)
write(lupr1,'('' mole fractions are:''/1p5e12.5)')
1 (rckmc(k0buf+k),k=1,kk)
endif
CALL MCADIF(PRES,RTMP,RCKMC(K0BUF+1),RCKMC(KMCR),
1 RCKMC(K0BUF+KK+1))
K0=L0DK+IZOFF*KK+ICELL
DO 9003 K=1,KK
F(K0)=RCKMC(K0BUF+KK+K)
K0=K0+NXNY
9003 CONTINUE
IF(SIUNIT) THEN
K0=L0DK+IZOFF*KK+ICELL
DO 90035 K=1,KK
F(K0)=F(K0)*DIFSCA
K0=K0+NXNY
90035 CONTINUE
ENDIF
ENDIF
K0=L0DK+IZOFF*KK
KIV=INDVAR-LBCKF
VALOUT=F(K0+KIV*NXNY+ICELL)
IF(MODDIF) THEN
VALOUT=VALOUT*SUMW
ENDIF
C
if(dbgloc.and.dbgam) then
call writ3i(' icell',icell,' isweep',isweep,
1 ' izstep',izstep)
write(lupr1,'('' diffusion coefficients:''/1p5e12.5)')
1 (rckmc(k0buf+kk+k),k=1,kk)
call writ1l(' moddif ',moddif)
cccc call writ2r(' sumw ',sumw,' pres ',pres)
call writ1r(' valout ',valout)
endif
C
IF(DOTDIF) THEN
C
C--- Calculate the thermal diffusion ratios and store in L0DTK
C
IF(INDVAR==LBCKF) THEN
CALL MCATDR(RTMP,RCKMC(K0BUF+1),ICKMC(KMCI),
1 RCKMC(KMCR),RCKMC(K0BUF+KK+1))
K0 =L0DTK+IZOFF*KK+ICELL
DO 9004 K=1,KK
F(K0)=RCKMC(K0BUF+KK+K)
K0=K0+NXNY
9004 CONTINUE
C
if(dbgloc.and.dbgam) then
call writ3i(' icell',icell,' isweep',isweep,
1 ' izstep',izstep)
write(lupr1,'('' thermal diffusion ratios:''
1 /1p5e12.5)') (rckmc(k0buf+kk+k),k=1,kk)
endif
ENDIF
ENDIF
ENDIF
ENDIF
ENDIF
GO TO 999
C-----------------------------------------------------------------
C
C--- GROUP 11. Initialization of variable or porosity fields
C Index VAL
11 CONTINUE
C For RHO1<=GRND--- density for phase 1
C
IF(NPATCH(1:6)=='ICHEMK') THEN
L0VAL=L0F(VAL)
IF(INDVAR==DEN1) THEN
PRES=PRESS0
RTMP=TEMP0
DO 11011 ICELL=1,NXNY
IF(VARPRS) PRES=PRESS0+F(L0P1+ICELL)
IF(VARTEM) RTMP=TEMP0+F(L0TEM+ICELL)
DO 11001 K=1,KK
RCKMC(K0BUF+K)=F(ICKMC(K0L0+K)+ICELL)
11001 CONTINUE
CALL CKRHOY(PRES,RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,RHOMIX)
F(L0VAL+ICELL)=RHOMIX
11011 CONTINUE
ELSEIF(INDVAR==W1) THEN
CALL POISSL(KVEL,IMUL)
ENDIF
ENDIF
C
GO TO 999
C-----------------------------------------------------------------
C
C--- GROUP 12. Convection and diffusion adjustments
C
c 12 CONTINUE
c GO TO 999
C-----------------------------------------------------------------
C
C--- GROUP 13. Boundary conditions and special sources
C Index for Coefficient - CO
C Index for Value - VAL
13 CONTINUE
C
IF(NPATCH(1:5)/='CHEMK'.AND.NPATCH(1:4)/='NOCP'.AND.
2 NPATCH(1:6)/='CHEMID'.AND.NPATCH(1:5)/='CKWAL') THEN
CALL WRIT40(' GXCHKI - Should not be here for PATCH ')
WRITE(LUPR1,'('' NAME = '',A8)') NPATCH
CALL SET_ERR(268,
* 'Error. GXCHKI - Should not be here for PATCH.',1)
C CALL EAROUT(1)
ENDIF
C
C!!!! Note to users: The next statement indicates how the Eddy-Break-Up
C model MIGHT be implemented; but it has not been tested. Users
C wishing to work with it should "uncomment".
c IF(NPATCH(1:6)=='CHEMK1'.OR.
c 1 NPATCH(1:6)=='CHEMK2'.OR.
c 1 NPATCH(1:6)=='CHEMK3') THEN
C
IF(NPATCH(1:6)=='CHEMK1') THEN
C
C--- Non-stiff / slightly-stiff sources
C CKCTYP - calculates rates of creation and characteristic
C destruction time-scales
C CKWYP - calculates rates of formation of species
C CKCTYR - calculates rates of creation and characteristic
C timescales for destruction of species
C
CALL SUB2(L0HVZ,L0HVAL+IZOFF*KK,L0HCZ,L0HCO+IZOFF*KK)
IF(ISC<11) THEN
IF(INDVAR==LBCKF) THEN
DO 1301 K=1,KK
F(L0DTF+K)=DTFALS(LBCKF+K-1)
1301 CONTINUE
IF(SOLTEM) F(L0DTF+KK+1)=DTFALS(ITIND)
CALL L0F1(CO,ICELL,IADD,'GXCHKI')
ICELL=ICELL+CO
PRES=PRESS0
RTMP=TEMP0
DO 1305 IX=IXF,IXL
ICELL=ICELL+IADD
DO 1304 IY=IYF,IYL
ICELL=ICELL+1
BLOCKD=.FALSE.
IF(IPRPS>0) BLOCKD=F(L0PRPS+ICELL)>=SOLPRP
C IF(TSTVPO) BLOCKD=BLOCKD.OR.F(L0VPOR+ICELL)<=1.E-10
IF(BLOCKD) THEN
K0C=L0CDK+IZOFF*KK
K0T=L0TDK+IZOFF*KK
DO 1391 K=1,KK
F(K0C+ICELL)=0.0
F(K0T+ICELL)=0.0
K0C=K0C+NXNY
K0T=K0T+NXNY
1391 CONTINUE
F(L0HVZ+ICELL)=0.0
F(L0HCZ+ICELL)=0.0
GO TO 1304
ENDIF
IF(VARPRS) PRES=PRES+F(L0P1+ICELL)
IF(VARTEM) RTMP=TEMP0+F(L0TEM+ICELL)
SUMY=0.0
SUMW=0.0
DO 1302 K=1,KK-1
RCKMC(K0VAR+K)=F(ICKMC(K0L0+K)+ICELL)
RCKMC(K0BUF+K)=0.0
RCKMC(K0BUF+KK+K)=0.0
SUMY=SUMY+RCKMC(K0VAR+K)
SUMW=SUMW+RCKMC(K0VAR+K)/RCKMC(K0MWT+K)
1302 CONTINUE
RCKMC(K0VAR+KK)=1.-SUMY
SUMW=SUMW+RCKMC(K0VAR+KK)/RCKMC(K0MWT+KK)
IF(SOLTEM) THEN
IF(ICKOP1==1.OR.ICKOP1==3) THEN
c---- calculate wdotk, i.e. molar production rates of species k
CALL CKWYP(PRES,1.001*RTMP,RCKMC(K0VAR+1),
1 ICKMC,RCKMC,RCKMC(K0BUF+KK+1))
ENDIF
C
HSORP=0.0
C---- calculate qdot = sum(wdotk*wtk*h0k) over all k
DO 1306 K=1,KK
HSORP=HSORP+RCKMC(K0BUF+KK+K)*RCKMC(K0MWT+K)*
1 RCKMC(K0H0K+K)
1306 CONTINUE
IF(RELQDT) HSORP=ALQDOT*HSORP-
1 (1.-ALQDOT)*F(L0QDOT+ICELL)
ENDIF
IF(ICKOP1==1.OR.ICKOP1==3) THEN
GRHO=F(L0DEN+ICELL)
CALL CKCTYR(GRHO,RTMP,RCKMC(K0VAR+1),ICKMC,RCKMC,
1 RCKMC(K0BUF+1),RCKMC(K0BUF+KK+1))
C CALL CKCTYP(PRES,RTMP,RCKMC(K0VAR+1),ICKMC,RCKMC,
C 1 RCKMC(K0BUF+1),RCKMC(K0BUF+KK+1))
ENDIF
RHO=F(L0DEN+ICELL)
HSOR=0.0
K0C=L0CDK+IZOFF*KK
K0T=L0TDK+IZOFF*KK
DTFT=F(L0DTF+KK+1)
DO 1303 K=1,KK
CDOTK=RCKMC(K0BUF+K)
DTF=F(L0DTF+K)
IF(RCKMC(K0BUF+KK+K)>GREAT.OR.
1 RCKMC(K0BUF+KK+K)<=0.0) THEN
TAUK = GREAT
ELSE
TAUK =RCKMC(K0BUF+KK+K)
ENDIF
IF(TAUK>1.E-20.OR.CDOTK>1.E-20) THEN
RHODT=RHO/TAUK
YCONR=RHO/(SUMW*RCKMC(K0MWT+K))
TAUCK=YCONR*(RCKMC(K0VAR+K)+1.E-4)/(TINY+CDOTK)
DTF=MIN(DTF,TAUCK)
DTF=MIN(DTF,TAUK)
F(L0DTF+K)=DTF
F(K0C+ICELL)=RHODT
F(K0T+ICELL)=CDOTK*RCKMC(K0MWT+K)/RHODT
IF(SOLTEM) HSOR=HSOR+RHODT*RCKMC(K0H0K+K)*
1 (F(K0T+ICELL)-F(ICKMC(K0L0+K)+ICELL))
ELSE
F(K0C+ICELL)=0.0
F(K0T+ICELL)=0.0
ENDIF
K0C=K0C+NXNY
K0T=K0T+NXNY
1303 CONTINUE
IF(SOLTEM) THEN
IF(RELQDT) HSOR=ALQDOT*HSOR-
1 (1.-ALQDOT)*F(L0QDOT+ICELL)
DHSBDT=(HSORP-HSOR)/(0.001*RTMP)
if(hsor>0.0) then
IF(DHSBDT<0.0) THEN
F(L0HVZ+ICELL)=RTMP-HSOR/DHSBDT
F(L0HCZ+ICELL)=-DHSBDT
ELSE
F(L0HVZ+ICELL)=-HSOR/FIXFLU
F(L0HCZ+ICELL)=FIXFLU
ENDIF
else
IF(DHSBDT>0.0) THEN
F(L0HVZ+ICELL)=RTMP-HSOR/DHSBDT
F(L0HCZ+ICELL)=DHSBDT
ELSE
F(L0HVZ+ICELL)=-HSOR/FIXFLU
F(L0HCZ+ICELL)=FIXFLU
ENDIF
endif
ENDIF
IF(L0QDOT>0) F(L0QDOT+ICELL)=-HSOR
1304 CONTINUE
1305 CONTINUE
ENDIF
IF(INDVAR==ITEM) THEN
CALL FN0(CO,-L0HCZ)
ELSEIF(INDVAR>=LBCKF.AND.INDVAR0.0) CALL FN34(LAP,LM1,RECT)
ENDIF
ELSE
IF(INDVAR==ITEM) THEN
CALL FN0(VAL,-L0HVZ)
ELSEIF(INDVAR>=LBCKF.AND.INDVAR=U1.AND.(EQZ(VCO)
1 .OR.VCO>=0.8*FIXVAL).OR.
1 INDVAR==P1.AND.QEQ(PCO,FIXFLU)) THEN
C
C--- Setting of inlet velocity: get incoming mass-flux
C
IF(GRN(PCO)) THEN
CALL WRIT40(' GXCHKI - Cannot do INLET with GRND P-CO')
CALL SET_ERR(269,
* 'Error. GXCHKI - Cannot do INLET with GRND P-CO.',1)
C CALL EAROUT(1)
ELSEIF(GRN(RHO1)) THEN
C
C--- Get the incoming gas composition and temperature to calculate the
C density of the incoming gas: note that the incoming temperature
C for GRND set enthalpy fluxes comes form TINLET
C
CALL GETCOV(NPATCH,ITIND,TCO,TVAL)
IF(GRN(TVAL)) THEN
CALL GETSPD(NPATCH,'TINLET',1,TINLET,IDUM,.FALSE.,
1 ' ',IERR)
IF(IERR==1) TINLET=TMP1A
RTMP=TEMP0+TINLET
ELSE
RTMP=TEMP0+TVAL
ENDIF
LCKERR=.FALSE.
SUMY=0.0
DO 13211 K=1,KK-1
CALL GETCOV(NPATCH,LBCKF+K-1,YCO,YVAL)
LCKERR=LCKERR.OR.GRN(YVAL)
RCKMC(K0BUF+K)=YVAL
SUMY=SUMY+YVAL
13211 CONTINUE
RCKMC(K0BUF+KK)=1.-SUMY
IF(LCKERR) THEN
CALL WRIT40(' GXCHKI - Cannot do INLET with GRND Y-CO')
CALL SET_ERR(270,
* 'Error. GXCHKI - Cannot do INLET with GRND Y-CO.',1)
C CALL EAROUT(1)
ELSE
PRES=PRESS0
C inlet species defined in mol fractions
IF(NPATCH(8:8)=='X') THEN
CALL CKRHOX(PRES,RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,RHOMIX)
ELSE
CALL CKRHOY(PRES,RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,RHOMIX)
ENDIF
ENDIF
RHOIN=REAL(RHOMIX)
ELSE
RHOIN=RHO1
ENDIF
C
C--- Now look to see if we have a Poisseulle flow PATCH; this is
C indicated by the 6th letter in the PATCH name, by having
C a non-zero value of the relevant PROFx variable, and by
C having the right PATCH-type etc. It is assumed that if the
C velocity is set as Poisseulle flow, then so also is the
C mass-flux (through P1). The work is all done in the subroutine
C POISSL
C
IF(BFC.AND.ISC==13) THEN
IF(INDVAR==P1.OR.(INDVAR>=U1.AND.INDVAR<=W1)) THEN
CALL FN1(VAL,RHOIN)
CALL GXBFC
ENDIF
ELSE
CALL POISSL(KVEL,IMUL)
C
IF(INDVAR==P1) THEN
C
C--- If this is the pressure, we are determining the mass-flux from
C the magnitude of the velocity multiplied by the density of the
C inflowing gas (NB.: the calculation is incorrect for BFC grids)
C
IF(KVEL/=W1) THEN
CALL FN1(VAL,0.0)
DO 13221 IVEL=U1,W1,2
CALL GETCOV(NPATCH,IVEL,VCO,VVALA)
IF(GRN(VVALA)) THEN
ELSEIF(NEZ(VVALA)) THEN
CALL FN33(VAL,VVALA*VVALA)
ENDIF
13221 CONTINUE
CALL FN30(VAL)
ENDIF
CALL FN25(VAL,RHOIN)
ELSEIF(.NOT.GRN(PVAL)) THEN
C
C--- If this is a velocity, we are determining the inflow-velocity
C from the mass-flux and the inflow density (note that this assumes
C only 1 non-zero inflow velocity)
C
CALL FN1(VAL,PVAL/RHOIN)
ELSE
CALL FN25(VAL,FLOAT(IMUL))
ENDIF
ENDIF
ELSEIF(INDVAR==ITIND) THEN
C
C--- Set the value of the incoming enthalpy NOCP PATCHes; get the
C incoming temperature form TINLET (always have GRND set enthapy
C flux for this option)
C
LCKERR=.FALSE.
CALL GETSPD(NPATCH,'TINLET',1,TINLET,IDUM,.FALSE.,
1 ' ',IERR)
IF(IERR==1) TINLET=TMP1A
RTMP=TEMP0+TINLET
SUMY=0.0
DO 13212 K=1,KK-1
CALL GETCOV(NPATCH,LBCKF+K-1,YCO,YVAL)
LCKERR=LCKERR.OR.GRN(YVAL)
RCKMC(K0BUF+K)=YVAL
SUMY=SUMY+YVAL
13212 CONTINUE
RCKMC(K0BUF+KK)=1.-SUMY
IF(LCKERR) THEN
CALL WRIT40(' GXCHKI - Cannot do INLET with GRND Y-CO')
CALL SET_ERR(271,
* 'Error. GXCHKI - Cannot do INLET with GRND Y-CO.',1)
C CALL EAROUT(1)
ELSE
CALL CKHBMS(RTMP,RCKMC(K0BUF+1),ICKMC,RCKMC,CPMIX)
C
C--- Must subtract off the heat of formation part of the enthalpy
C (Note: CPMIX is used to store the enthalpy here)
C
DO 13213 K=1,KK
CPMIX=CPMIX-RCKMC(K0BUF+K)*RCKMC(K0H0K+K)
13213 CONTINUE
CALL FN1(VAL,REAL(CPMIX))
ENDIF
ENDIF
C
ELSEIF(NPATCH(1:6)=='CHEMID') THEN
C
C--- Inlet diffusion PATCHes - the enthalpy flux due to mass-diffusion
C is summed and the added into the enthalpy source
C
IF(INDVAR>=LBCKF.AND.INDVAR<=LBCKF+KK-2.AND.
1 ISC==10) THEN
IPV=L0PVAR(PVCKID,IREG,IZSTEP)
IF(INDVAR==LBCKF) CALL FNSCPA(-IPV,0.0)
CALL GETCOV(NPATCH,INDVAR,YCO,YVAL)
CALL GETSPD(NPATCH,'TINLET',1,TINLET,IDUM,.FALSE.,
1 ' ',IERR)
IF(IERR==1) TINLET=TMP1A
CALL L0F2(CO,INDVAR,ICO,IFIC,IADD,'GXCHKI')
DO 13301 IX=IXF,IXL
ICO=ICO+IADD
DO 13301 IY=IYF,IYL
ICO=ICO+1
IPV=IPV+1
RTMP=TEMP0+TINLET
CALL CKHMS(RTMP,ICKMC,RCKMC,RCKMC(K0BUF+1))
F(IPV)=F(IPV)+MAX(F(ICO)*(YVAL-F(ICO+IFIC)),0.0)
1 *(RCKMC(K0BUF+K)-RCKMC(K0BUF+KK))
13301 CONTINUE
ELSEIF(INDVAR==ITIND.AND.ISC==21) THEN
CALL L0F1(VAL,IVAL,IADD,'GXCHKI')
IPV=L0PVAR(PVCKID,IREG,IZSTEP)
DO 13302 IX=IXF,IXL
IVAL=IVAL+IADD
DO 13302 IY=IYF,IYL
IVAL=IVAL+1
IPV=IPV+1
F(IVAL)=F(IPV)
13302 CONTINUE
ENDIF
C
ELSEIF(MODDIF.AND.NPATCH(1:5)=='CKWAL'.AND.ISC==10.AND.
1 INDVAR>=LBCKF.AND.INDVAR=LBCKF.AND.INDVAR
SUBROUTINE PHINDX(NGO,LINK,LOUT,PRECIS,NTDA,K0TD)
C
C---------------------------------------------------------------
C PHINDX - Checks that the precision of GXCHKI matches that of
C the CHEMKIN and TRANLIB data-bases; also sets various
C indices for GXCHKI
C
C returns .TRUE. if there are thermal diffusion
C coeff.s
C---------------------------------------------------------------
C
C*****double precision
DOUBLE PRECISION
C*****END double precision
C*****single precision
C REAL
C*****END single precision
1 RU, PATMOS, SMALL
CHARACTER*40 PRECIS
CHARACTER*16 PREC, VERS, PRECMC, VERSMC
LOGICAL ERROR, KERR, KERRMC
COMMON /MCCONS/ VERSMC, PRECMC, KERR, LENI, LENR
COMMON /CKCONS/ VERS, PREC, KERRMC, LENIMC, LENRMC
C
COMMON /MCMCMC/ RU, PATMOS, SMALL, NKK, NO, NLITE, INLIN, IKTDIF,
1 IPVT, NWT, NEPS, NSIG, NDIP, NPOL, NZROT, NLAM,
2 NETA, NDIF, NTDIF, NXX, NVIS, NXI, NCP,
3 NCROT, NCINT, NPARK, NBIND, NEOK, NSGM,
4 NAST, NBST, NCST, NXL, NR, NWRK, K3
C
IF(NGO==1) THEN
C
C--- Check CHEMKIN data-base
C
ERROR = PREC(1:6)/=PRECIS(2:7)
C
ELSEIF(NGO==2) THEN
C
C--- Check the TRANLIB data-base; get the number of species for which
C there are thermal diffusion ratios/coefficients
C
ERROR = PREC(1:6)/=PRECIS(2:7)
C
NTDA=NLITE
K0TD=IKTDIF-1
C
ENDIF
C
END
C
C===============================================================
C
SUBROUTINE POISSL(KVEL,IMUL)
C
C-----------------------------------------------------------------
C SUBROUTINE POISSL - puts a POISSL flow into the VAL slot for
C group 13 or group 11 PATCHes
C-----------------------------------------------------------------
C Arguments:
C KVEL - returns a code that indicates whether or not a
C setting has been made
C-----------------------------------------------------------------
C
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
C
LOGICAL GRN
C
IMUL=0
KVEL=0
IF((INDVAR==P1.AND.IGR==13).OR.INDVAR==W1) THEN
VPOISS=GRND
IF(NPATCH(7:7)=='A') THEN
VPOISS=PROFA
ELSEIF(NPATCH(7:7)=='B') THEN
VPOISS=PROFB
ELSEIF(NPATCH(7:7)=='C') THEN
VPOISS=PROFC
ENDIF
KVEL=0
IF(.NOT.GRN(VPOISS)) THEN
IF(IGR==13.AND.(INTTYP==6.OR.INTTYP==7).OR.
1 IGR==11.AND.IPATNF(2,IREG)==26) THEN
KVEL=7
ENDIF
IF(IGR==13) THEN
IGEO=PATGEO
ELSEIF(IGR==11) THEN
IGEO=ANORTH
ENDIF
IF(VPOISS>0.0) THEN
IDVP=YG2D
IF(IYF==1) THEN
XYF=0.0
ELSE
XYF=F(L0F(YV2D)+IYF-1)
ENDIF
XYL=F(L0F(YV2D)+IYL)
IF(.NOT.CARTES) THEN
XYF=XYF+RINNER
XYL=XYL+RINNER
ENDIF
ELSE
VPOISS=-VPOISS
IDVP=XG2D
IF(IXF==1) THEN
XYF=0.0
ELSE
XYF=F(L0F(XU2D)+(IXF-2)*NY)
ENDIF
XYL=F(L0F(XU2D)+(IXL-1)*NY)
ENDIF
IF(IGR==13) IMUL=2*MOD(INTTYP,2)-1
SUM=0.0
ASUM=0.0
CALL L0F3(VAL,IDVP,IGEO,I,IDMI,IGMI,IADD,'GXCHKI')
IF(XYF>0.0) RLRAT=(XYL**2-XYF**2)/LOG(XYL/XYF)
DO 13215 IX=IXF,IXL
I=I+IADD
DO 13215 IY=IYF,IYL
I=I+1
F(I)=XYL**2-F(I+IDMI)**2
IF(XYF>0.0) F(I)=F(I)+RLRAT*LOG(F(I+IDMI)/XYL)
SUM=SUM+F(I)*F(I+IGMI)
ASUM=ASUM+F(I+IGMI)
13215 CONTINUE
CALL FN25(VAL,VPOISS*ASUM/SUM)
ENDIF
ENDIF
END
c