cc
C.... File Name ................ GREX3.FOR ...................... 280623
c
C This subroutine is called by EARTH within the solution cycle.
C
C Subroutine GREX3, with the subsidiary subroutines GXxxxx,
C allows the PHOENICS EARTH program to access a wide range of
C special physical models, algorithm adjustments and print-out
C features.
C
C 'EX' stands for 'example', signifying that the coding practices
C exhibited here may be used as patterns by PHOENICS users who
C wish to introduce their own models, algorithms, etc.
C
C The GX subroutines that are called in group 13 are activated
C by settings in the SATELLITE, mainly by way of PATCH names.
C Many of them, like GREX3 itself, call on 'utility' subroutines'
C named FNxxx.
C
C The FNxxx subroutines perform various arithmetic operations on
C segments of the F-array corresponding to variables of like kind,
C indicated by the arguments. The FUNCT entry of the Encyclopaedia
C provides full explanations.
C
c
SUBROUTINE GREX3
C--------------------------------------- INCLUDE COMMON BLOCK files
INCLUDE 'patcmn'
INCLUDE 'objnam'
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
INCLUDE 'parear'
INCLUDE 'lbnamer'
CHARACTER*48 COMMAND
COMMON/IMAIN2/IMAIN1(9),ISWSTP,IMAIN2(33)
COMMON/F0/KF01(70),KAP,KF02(233)
1 /MRKCMN/LBMARK,L0MARK
COMMON/LBFC/STORSA(6),STORWD(6),LBFCSP /ARRT/ARRIVT
LOGICAL STORSA,STORWD,LBFCSP,LGTEMP,SPCLGR,ARRIVT
COMMON/RDLOG/LOGQ(10)/RDINT/INTQ(10)/RDREAL/REALQ(10)
COMMON/RAKEEP/RADIAT/DVMOD/IDVCGR
LOGICAL RADIAT,NEZ,dbgloc
COMMON /LRNTM3/L0UTAU,NMWALL,L0WALL,L0DSKN,IVPRST
LOGICAL NOCALC,DRAGON,DRG1ST,PBFC,DFAIL,CALLSL
COMMON/DRAGCM/NOCALC,DRAGON,DRG1ST,PBFC,DFAIL
COMMON /LUNITS/LUNIT(60)
1 /GENL/LFIL1(3),TIMPRN,LGF1(10),debgtz,SWPRIN,LGF2(40),
1 TURB,LGF3(3)
LOGICAL SLIPVL,LG,KEBUOY,LGF1,debgtz,TIMPRN,SWPRIN,LGF2,TURB,LGF3,
1 HOCS,ESCRS,FSQSOR,RADCVD,BLIN,LGXIO,LSUN,LWAVE,SURFT
SAVE SLIPVL,NUMSPPR,KEBUOY,SPCLGR,ESCRS,FSQSOR,HOCS,PSIVEL,
1 BLIN,LWAVE,SURFT
COMMON/GRXFLAGS/LGXIO,LSUN
COMMON/LBINGREX/LBSLU,LBSLV,LBSLW,LBWAVE,LBGENK,LBGNK2,LBMACH,
1 LBPTOT,LBVABS,LBVLSQ,LBPTO2,LBVAB2,LBMAC2,LBFSQ,
1 LBARRT,LBWCRT,LBPSIV,LBVOXY,LBVOYZ,LBVOZX
COMMON/LRNTM1/L0WDIS,L0FMU,L0FONE,L0FTWO,L0REYN,L0REYT,L0UD1,L0UD2
1 ,L0UD3,L0UD4
COMMON/PASQLF/ITPRO,PASQBUOY,BUOSSG,MOLEN
LOGICAL PASQBUOY,BUOSSG
COMMON/EXPLOSION/EXPLOS
LOGICAL EXPLOS,PSIVEL
C.... Provision for the use of GRAPH
DIMENSION TIME(100),OR(100,10),LABOR(10),IOSC(10)
CHARACTER*4 LABOR,CG
PARAMETER (NLG=100, NIG=100, NRG=200, NCG=100)
COMMON/LGRND/LG(NLG)/IGRND/IG(NIG)/RGRND/RG(NRG)/CGRND/CG(NCG)
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
COMMON/GENI/NXNY,IFIL1(8),NFM,IGFIL1(21),IPRL,IBTAU,ILTLS,
1 IGFIL(15),ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,
1 IRADZ,IVFOL
PARAMETER (MXXLET=20)
COMMON /LETL0C/ L0C(MXXLET)
CHARACTER*(MXXLET) CINDEX, WREF*2
COMMON /LETBLK/ CINDEX
COMMON/IELAPS/IREC,ITIME1,NUMSOL
COMMON/GRX/NUMSEC
COMMON/RHILO/HI3D,RLO3D
LOGICAL DONE192,GXMON,SLD,DOIT
COMMON/GXMON2/PLTCLSX; LOGICAL PLTCLSX
SAVE DONE192
SAVE TIME,OR,LABOR
DIMENSION DIGIT(10)
CHARACTER*2 DIGIT
DATA DIGIT/'1','2','3','4','5','6','7','8','9','10'/
DATA DONE192 /.FALSE./
c
dbgloc=debgtz.and.dbgrnd
if(dbgloc) then
call writbl
call banner(1,'grex ',120121)
write(14,*)'entry to grex for: igr,isc,indvar ',igr,isc,indvar
endif
if(test) call fcheck(0,'grex3 ')
c----------------------------------------------------------------------
NAMSUB='GREX3'
C.... GO TO the group indicated by IGR, then to section ISC
IF(IGR==19) THEN
GO TO 19
ELSEIF(IGR==13) THEN
GO TO 13
ELSE
GO TO (1,2,3,4,5,6,25,8,9,10,11,12,13,14,25,25,25,25,19,
1 20,21,25,25,24) IGR
ENDIF
C****************************************************************
C
C---- GROUP 1. Run title and other preliminaries
C
1 GO TO (1001,1002,1003),ISC
C
C.... group 1 sections are called in the order 1, 2, 3.
C****************************************************************
C
C * -----------GROUP 1 SECTION 1 ------ Initializations and
C allocations of storage for variables which it is essential to
C dump to the PHI (or PHIDA) file for restart purposes
C
1001 IX=1 ; IY=1
C----------------------------------- Identification of GREX3 to VDU
IF(.NOT.NULLPR.AND.UWATCH)
1 CALL WRYT40('GREX3 of 280623 has been called')
C
IF(READQ1) CALL GXRDQ1 ! Read Q1EAR file when READQ1 = T in Q1
c
c.... Set absolute value below which zeroes will appear in field-value
c print-out (unless debug=t, which sets print0 to 0.0).
c Users wishing to change print0 without re-compiling should insert
c PRT0BEGIN
c PRINT0 [desired value]
c PRT0END
c at the top of their Q1 file
c
IF(.NOT.DEBUG) PRINT0=1.E-11
CALL RQ1R('PRT0','PRINT0',PRINT0)
C.... Set up PATCH-wise storage, if necessary, for GXBFC
IF(BFC) CALL GXBFC
PSIVEL=LBPSIV/=0
C------------------------ PHI and XYZ files for parabolic cases
C
IF(PARAB.OR.(.NOT.STEADY.AND.PNAM==' ')) THEN
IF(IDISPA>0) CALL GXPARA
ENDIF
NUMSEC=0
C.... Initializations dependent on the existence of patch names
RADIAT=.FALSE. ; KEBUOY=.FALSE. ; HOCS=.FALSE. ; FSQSOR=.FALSE.
BLIN=.FALSE. ;LGXIO=.FALSE.; LSUN=.FALSE.; LWAVE=.FALSE.
SURFT=.FALSE. ; PASQBUOY=.FALSE.
DO 10011 I = 1,NUMPAT
IF(NAMPAT(I)(1:6)=='KEBUOY') THEN
KEBUOY=.TRUE.
ELSEIF(NAMPAT(I)(1:4)=='RADI') THEN
RADIAT=.TRUE.
CALL GXRADI(CARTES,TEMP1,DEN1,RHO1) ! Radiation initialization
ELSEIF(NAMPAT(I)(1:4)=='HOCS') THEN
HOCS=.TRUE.
ELSEIF(NAMPAT(I)(1:4)=='STEN') THEN
SURFT=.TRUE.
ELSEIF(NAMPAT(I)(1:5)=='FSQSO') THEN
FSQSOR=.TRUE.
ELSEIF(NAMPAT(I)(1:4)=='BLIN'.OR.NAMPAT(I)(1:4)=='GXBL')THEN
BLIN=.TRUE. ! Boundary-layer inlet profile facility
ELSEIF(.NOT.LWAVE.AND.NAMPAT(I)(1:4)=='WAVE'.AND.IVFOL>0) THEN
CALL GETCOV(NAMPAT(I),IVFOL,GCO,GVAL)
IF(GCO/=-999.) LWAVE=.TRUE.
ELSEIF(.NOT.BFC) THEN
COMMAND=' '
CALL GETSDC(OBJNAM(I),'INLET',COMMAND)
IF(COMMAND==' ') CALL GETSDC(OBJNAM(I),'OUTLET',COMMAND) ! test for
! angled-out object
IF(COMMAND=='GXI'.OR.COMMAND=='GXO') LGXIO=.TRUE.
ENDIF
10011 CONTINUE
LGXIO=LGXIO.OR.(MIMD.AND.NPROC>1)
COMMAND=' '
CALL GETSDC('SUNLIGHT','TIME',COMMAND)
IF(COMMAND/=' ') LSUN=.TRUE.
C---------------------------------- SLIPVL setting
IF(.NOT.ONEPHS) THEN
SLIPVL = LBNAME('SLPU')/=0.OR.LBNAME('SLPV')/=0.OR.
1 LBNAME('SLPW')/=0
IF(SLIPVL) THEN
LBSLU=0 ; LBSLV=0 ;LBSLW=0
IF(SOLVE(3).AND.SOLVE(4)) LBSLU=LBNAME('SLPU')
IF(SOLVE(5).AND.SOLVE(6)) LBSLV=LBNAME('SLPV')
IF(SOLVE(7).AND.SOLVE(8)) LBSLW=LBNAME('SLPW')
ENDIF
IF(NEZ(CPIP)) CALL GXINTP
IF(NEZ(CLIFT)) CALL GXLIFT
ENDIF
IF(BLIN) CALL GXBLIN ! Boundary-layer inlet profile facility
IF(FDFSOL) CALL GXFDFS ! Fully-developed floq
IF(HOCS) CALL GXHOCS
IF(NOX1) CALL GXNOX
IF(NOX2) CALL GXNOX
IF(CHMKIN) CALL GXCHKI(0,0.0)
!!! IF(LSG87) CALL GXTCALC(0,0.0)
IF(LBARRT/=0) CALL GXARRT
IF((GRNDNO(9,TMP1).OR.GRNDNO(10,TMP1)).AND..NOT.LSG87) THEN
ESCRS=.TRUE.
IF(GRNDNO(9,TMP1)) THEN
C Define operating pressure in atmospheres via CHSOC for CHEMKIN
CHSOC=PRESS0/1.01325E5
CALL GXCHKI(0,0.0)
ENDIF
CALL GXSCRI
ELSE
ESCRS=.FALSE.
ENDIF
EXPLOS=.FALSE.
C***************************************************************
CALL UCASZZ(NAMGRD,48)
SPCLGR= YPLS.OR.WALPRN.OR.S2SR
1 .OR.NAMGRD=='FURN'.OR.NAMGRD=='ESTR'.OR.NAMGRD=='FLAR'
1 .OR.NAMGRD=='HTBX'.OR.NAMGRD=='TACT'.OR.NAMGRD=='MICA'
1 .OR.NAMGRD=='CONV'.OR.NAMGRD=='CVD '.OR.NAMGRD=='F1'
1 .OR.NAMGRD=='WATS'.OR.NAMGRD=='WAVE'
LL=LENGZZ(NAMGRD)
IF(SPCLGR.AND.NAMGRD/='NONE') THEN
IF(.NOT.NULLPR) CALL WRIT40
1 ('Special GROUND '//NAMGRD(1:LL)//' is in use ')
ELSEIF(NAMGRD/='NONE') THEN
CALL WRITST
CALL WRITBL
CALL WRIT40
1 ('Stop because '//NAMGRD(1:LL)//' is not known to GREX3 ')
CALL WRYT40
1 ('Stop because '//NAMGRD(1:LL)//' is not known to GREX3 ')
CALL SET_ERR(213,'Error. See result file.',1)
ENDIF
IF(MYID==0) CALL LATEST_DUMP(1,0,0,' ',' ',' ')
GO TO 25
C
C***********************************************************************
C
C * -----------GROUP 1 SECTION 3 ---------------------------
c
C---- Use Section 3 to create storage via MAKE or GXMAKE which it is
C not necessary to dump to the PHI(da) file for restarts.
C
1003 CONTINUE
C
C.... Provision of special EARTH arrays by means of the
C storage-setting subroutine MAKE, used thus:
C CALL MAKE(variable name)
C
C---------------------------------------------------------BFC=T
IF(BFC) THEN
CALL GXBFC
C.... Segment EASP4 is used for buoyancy sources when BFC=T
C (see GXBUOY)
CALL MAKE(EASP4)
C
ENDIF
IF(SURFT) CALL GXSURFT
C
C-------------------------------------------------EASP1, EASP2, EASP3
C
C.... F-array segments with zero-locations EASP1, EASP2 & EASP3 are
C always available, for use as local temporary stores
C (See for example their uses in GXWALL & GXPORA)
CALL MAKE(EASP1)
CALL MAKE(EASP2)
CALL MAKE(EASP3)
CALL MAKE(EASP5) ! to store SKINT in gxwall.for
IF(UCONV.AND.NZ>1) CALL MAKE(EASP20)
C
C-------------------------------------------------------EASP7
C.... Store for k-eps source-term linearisation
IF(BFC.OR.KELIN==1) CALL MAKE(EASP7)
C
C---------------------------------------------------Letter Mask
MXL=MXXLET
CINDEX=' '
DO 111 IPAT=1,NUMPAT
KK=INDEX(NAMPAT(IPAT)(:2),'%')
IF(KK/=0) THEN
WREF=NAMPAT(IPAT)(KK+1:KK+2)
C* WREF(:1) = $ LETTER MASK IN XY PLANE
C* = & LETTER MASK IN ZY PLANE
C* = # LETTER MASK IN XZ PLANE
IF(WREF(:1)=='$' .OR. WREF(:1)=='&'
1 .OR. WREF(:1)=='#') THEN
K=INDEX(CINDEX,WREF(2:2))
IF(K==0) CALL LETMSK(CINDEX,WREF(2:2),L0C,MXL)
ENDIF
ENDIF
111 CONTINUE
C---------------------------------------------------------LGEN1
C.... Provision of storage for square-of-velocity-gradient
C expressions for viscous-dissipation or turbulence-
C energy source.
C
IF(GENK .OR. DUDX.OR.DVDX.OR.DWDX.OR.DUDY.OR.DVDY.OR.DWDY.OR.
1 DUDZ.OR.DVDZ.OR.DWDZ) THEN
IF(BFC) THEN
LIMIT1=1
LIMIT2=18
ELSE
LIMIT1=5
LIMIT2=9
ENDIF
DO 1000 I= LIMIT1,LIMIT2
CALL MAKE(LXYEA(I))
1000 CONTINUE
ENDIF
IF((SOLVE(H1).AND.(GRNDNO(5,TMP1).OR.GRNDNO(6,TMP1)))
1 .OR.(SOLVE(H2).AND.(GRNDNO(5,TMP2).OR.GRNDNO(6,TMP2)))) NCRT=1
IF(IENUTA.EQ.14) CALL GX_REALISABLE_KE
IF(IENUTA.EQ.10.OR.IENUTA.EQ.11) CALL GXKW_WILCOX ! not needed.
IF(IENUTA.EQ.15.OR.IENUTA.EQ.16) CALL GXKW_WILCOX
IF(IENUTA.GE.17.AND.IENUTA.LE.20) CALL GXKW_MENTER
GO TO 25
C
C***********************************************************************
C * -----------GROUP 1 SECTION 2 Further initializations which
C require that storage shall aready have been allocated
c
1002 CONTINUE
c call fcheck(10,namsub)
IF(KEBUOY) CALL GXKEGB
c call fcheck(20,namsub)
namsub='grex3'
IF(FSQSOR) CALL GXFSQS
C.... Wall-function initializations
IF(TURB .AND. NMWALL>0) CALL GXWFUN
IF(HOCS) CALL GXHOCS
IF(.NOT.ONEPHS) LBWAVE=LBNAME('WAVE')
C.... Possibly set BFC grid by way of coding in GXBFGR
IF(SETBFC.OR.MOVBFC) CALL GXBFGR
IF(LGXIO) CALL GXIO
IF(LSUN) CALL GXSUN
IF(BLIN) CALL GXBLIN ! Boundary-layer inlet profile facility
CALL GXOUTLET
CALL WRTPAT
GO TO 25
c
C*****************************************************************
C
C
C--- GROUP 2. Transience; time-step specification
C
2 CONTINUE
C * Set DT if TLAST has been made <=GRND in satellite
C... Set timestep from Courant Number
IF(TLAST-TFIRST==GRND1) CALL GXTIMESTEP
GO TO 25
C*****************************************************************
C
C
C--- GROUP 3. X-direction grid specification
C
3 CONTINUE
C * Set XRAT if AZXU has been made <=GRND in satellite
GO TO 25
C*****************************************************************
C
C
C--- GROUP 4. Y-direction grid specification
C
4 CONTINUE
C * Set YRAT if AZYV has been made <=GRND in satellite
C
GO TO 25
C*****************************************************************
C
C
C--- GROUP 5. Z-direction grid specification
C
5 CONTINUE
C * Set DZ if AZDZ has been made <=GRND in satellite
IF(GRNDNO(1,AZDZ)) THEN ! azdz = grnd1 or propx
DZ=DZW1*XULAST ! i.e. proportional to xulast
ELSEIF(GRNDNO(2,AZDZ)) THEN ! azdz = grnd2 or propy
DZ=DZW1*YVLAST ! i.e. proportional to yvlast
ENDIF ! dzw1 = proportionality factor
c * Insert DZ in the appropriate slab-wise array
CALL FN1(LXYDZ,DZ)
GO TO 25
C*****************************************************************
C
C
C--- GROUP 6. Body-fitted coordinates or grid distortion
C This group is visited when UGEOM is set .TRUE.. The visits
C occur at the start of each z slab after its geometry has been
C computed. Hence, this is the right place to access AEAST, AHIGH,..
C etc.
6 CONTINUE
C
GO TO 25
C*****************************************************************
C * Make changes for this group only in group 19.
C
C--- GROUP 7. Variables stored, solved & named
C*****************************************************************
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,816)
1 ISC
C * Add velocities relative to grid for phase 1 and/or phase 2
C * -----------GROUP 8 SECTION 1 ---------------------------
C--- for U1AD<=GRND--- phase 1 additional velocity (VELAD).
81 CONTINUE
C--- for U1AD<=GRND--- phase 1 additional velocity (VELAD).
c.... The option U1AD=GRND1 adds to U1 a component equal -W1*TAN(ALF)
c where ALF is the east-face grid angle. It is needed for the
c parabolic solver for grids for which constant-x lines are
c not precisely at right angles to constant-z lines. The
c correction is also available for U2, V1 and V2 (see below).
c
IF(GRNDNO(1,U1AD)) CALL GXPVAD(XULAST,XRAT,XU2D,W1,'X')
GO TO 25
C * -----------GROUP 8 SECTION 2 ---------------------------
C--- for U2AD<=GRND--- phase 2 additional velocity (VELAD).
82 CONTINUE
IF(GRNDNO(1,U2AD)) CALL GXPVAD(XULAST,XRAT,XU2D,W2,'X')
GO TO 25
C * -----------GROUP 8 SECTION 3 ---------------------------
83 CONTINUE
C--- for V1AD<=GRND--- phase 1 additional velocity (VELAD).
IF(GRNDNO(1,V1AD)) CALL GXPVAD(YVLAST,YRAT,YV2D,W1,'Y')
GO TO 25
C * -----------GROUP 8 SECTION 4 ---------------------------
C--- for V2AD<=GRND--- phase 2 additional velocity (VELAD).
84 CONTINUE
IF(GRNDNO(1,V2AD)) CALL GXPVAD(YVLAST,YRAT,YV2D,W2,'Y')
GO TO 25
C * -----------GROUP 8 SECTION 5 ---------------------------
C--- for W1AD<=GRND--- phase 1 additional velocity (VELAD).
85 CONTINUE
GO TO 25
C * ----------GROUP 8 SECTION 6 ---------------------------
C--- for W2AD<=GRND--- phase 2 additional velocity (VELAD).
86 CONTINUE
GO TO 25
C * -----------GROUP 8 SECTION 7 ---- VOLUMETRIC SOURCE FOR GALA
C---- Entered when GALA =.TRUE.; block-location index is LSU or LSU2
87 CONTINUE
GO TO 25
C * -----------GROUP 8 SECTION 8 --- CONVECTION COEFFICIENTS
C--- Entered when UCONV =.TRUE.
88 CONTINUE
c call fcheck(881,'grex3 ')
C.... Correction applicable to velocities in compressible flow
IF(COMPRS) THEN ! COMPRS is a PIL variable
IF(INDVAR>=U1.AND.INDVAR<=W2) THEN ! do this only for velocities
IF((INDVAR==U1.OR.INDVAR==U2).AND.NDIREC==3) THEN
CALL GXCMPR(LD11,EAST(P1),P1,1.4)
CALL GXCMPR(LD12,P1,EAST(P1),1.4)
ELSEIF((INDVAR==V1.OR.INDVAR==V2).AND.NDIREC==1) THEN
CALL GXCMPR(LD11,NORTH(P1),P1,1.4)
CALL GXCMPR(LD12,P1,NORTH(P1),1.4)
ELSEIF(INDVAR==W1.OR.INDVAR==W2) THEN
IF(NDIREC==5) THEN
CALL GXCMPR(EASP20,HIGH(P1),P1,1.4)
ELSEIF(NDIREC==6) THEN
CALL GXCMPR(EASP20,P1,HIGH(P1),1.4)
ENDIF
ENDIF
ENDIF
ENDIF
c call fcheck(882,'grex3 ')
C
GO TO 25
C * -----------GROUP 8 SECTION 9 --- DIFFUSION COEFFICIENTS
C.... Entered when UDIFF =.TRUE.; block-location indices are:
C LAE for east, LAN for north, and LD11 for high.
C The user should provide INDVAR and NDIREC IF's as appropriate.
C EARTH will apply the DIFCUT and GP12 modifications after the user
C has made his settings.
C
89 CONTINUE
IF(CHMKIN) CALL GXCHKI(0,0.0)
GO TO 25
C * -----------GROUP 8 SECTION 10 --- CONVECTION NEIGHBOURS
C---- Entered when UCONNE =.TRUE.; block-location indices are LD7
810 CONTINUE
GO TO 25
C * -----------GROUP 8 SECTION 11 --- DIFFUSION NEIGHBOURS
C---- Entered when UDIFNE =.TRUE.; block-location indices are LD7
C for south, west, high or low, and LD8 for north or east.
C User should provide INDVAR and NDIREC IF's as above.
811 CONTINUE
IF(CHMKIN) CALL GXCHKI(0,0.0)
GO TO 25
C * -----------GROUP 8 SECTION 12 --- LINEARISED SOURCES
C---- Entered when USOURC =.TRUE.
812 CONTINUE
IF(FDFSOL) CALL GXFDFS
IF(CHMKIN) CALL GXCHKI(0,0.0)
GO TO 25
C * -----------GROUP 8 SECTION 13 --- CORRECTION COEFFICIENTS
C---- Entered when UCORCO =.TRUE.
813 CONTINUE
GO TO 25
c
C * -----------GROUP 8 SECTION 14 --- USER'S SOLVER
C---- Entered when USOLVE =.TRUE.;
814 CONTINUE
c........................ SLVR and CSG3 are EQUIVALENCEd in GRDLOC
C------------------------- Call Gauss-Seidel solver
IF(SLVR=='GAUS'.AND..NOT.SLBSOL) CALL GXGAUS(OVRRLX,
1 LITER(INDVAR),IXMON,IYMON,IZMON,ENDIT(INDVAR),XCYCLE,LUPR1,LUPR3)
C
USER=.FALSE.
C
CALLSL = INDVAR==1 .OR. INDVAR==3 .OR. INDVAR==5 .OR.
1 INDVAR==7 .OR.(INDVAR>=ICNGRB .AND. INDVAR<=ICNGRC)
IF(CALLSL) THEN
NZZ=NZ
IF( SLBSOL ) NZZ=1
IF( SLVR=='CGGR') THEN
CALL GXCGGR(NX,NY,NZZ)
ELSEIF( SLVR=='CRGR') THEN
CALL GXCRGR(NX,NY,NZZ)
ENDIF
C.... entry to GXCHKI for N-R solvers
ELSEIF(CHMKIN) then
CALL GXCHKI(0,0.0)
ENDIF
C
GO TO 25
C * -----------GROUP 8 SECTION 15 --- Change solution result
C---- Entered when UCORR =.TRUE.; block-location indices are as above.
815 CONTINUE
C.... entry to GXCHKI for N-R solvers
C
IF(CHMKIN) CALL GXCHKI(0,0.0)
GO TO 25
816 CONTINUE
C * ------------------- SECTION 16 --- Change DVEL/DPs if UDVDP = .TRUE.
GO TO 25
C * Make all other group-8 changes in group 19.
C***********************************************************************
C
C
C--- GROUP 9. Properties of the medium (or media)
C
C... EARTH no longer calls Group 9 and 10 of GREX for properties. Instead
C it directly calls subroutine SLBPRP (in EARTH) which itself calls the
c the appropriate one of the following:
c SLBDEN (in GXDENS.HTM) for DENST1
c SLBDRP CMPRS1
c SLBTHE THRME1
c SLBVST VISTRB
c SLBVSL VISCLM
c SLBPRL PRNLAM
c SLBTMP TEMPR1
c SLBSPH SPEHT1
c SLBLEN MIXLN1
c SLBDEN DENST2
c SLBDRP CMPRS2
c SLBTHE THRME2
c SLBTMP TEMPR2
c SLBSPH SPEHT2
c SLBLEN MIXLN2
c SLBIVL INTVL1
c SLBFRC INTFCO
c SLBIMS INTMTR
c SLBIDF INTRC1
c SLBVRM IVRMCO
c
C property flag is set to one of the GRND values. However, if a property
C flag is set to GRND, SLBPRP calls the GROUND subroutine, expecting to
C find there user-provided coding for the property in question.
C
9 GO TO 25
C***********************************************************************
C
C
C--- GROUP 10. Inter-phase-transfer processes and properties
C
C... See comment to Group 9.
10 GO TO 25
C***********************************************************************
C
C
C--- GROUP 11. Initialization of variable or porosity fields
C
11 IF(NPATCH(1:4)=='IBFC') THEN ! Initial field of uniform
CALL GXBFC ! flow in a curvilinear grid
ELSEIF(NPATCH(1:6)=='INIPOL') THEN
CALL GXIPOL
ELSEIF(NPATCH(1:6)=='ICHEMK') THEN
IF(CHMKIN) CALL GXCHKI(0,0.0)
ELSEIF(NPATCH(1:4)=='BLIN'.OR.NPATCH(1:4)=='GXBL') THEN ! Boundary-layer profile
CALL GXBLIN ! Boundary-layer inlet profile facility
ENDIF
GO TO 25
C*****************************************************************
C
C
C--- GROUP 12. Convection and diffusion adjustments
C
12 CONTINUE
GO TO 25
C*****************************************************************
C
C
C--- GROUP 13. Boundary conditions and special sources
C
C * XCYCLE may be changed in group 19.
C
C Sections 1 to 11 for CO
C Sections 12 to 22 for VAL
13 NPAT= NPATCH(1:4)
c.... "Un-comment" call to accutm('grp 13') in three places below, in
c order to compute and print cumulative time spent in group 13.
c Such calls can be placed anywhere; but, since accutm itself takes
c a not-insignificant amount of computer time, they are best
c "commented" when not needed
c
c call accutm('grp 13',1)
c
IF(ISVROB) THEN ! Momentum sources due to swirling fan object
IF(NAMGRD/='HTBX'.AND.NAMGRD/='FLAR') CALL SWIRFAN
ENDIF
IF(NPATCH(1:4)=='BLIN'.OR.NPATCH(1:4)=='GXBL') THEN ! Boundary-layer inlet profile
CALL GXBLIN
ELSEIF(NPATCH(1:3)=='FIR') THEN
CALL GXFIRE ! FIRE
ELSEIF(NPATCH(1:6)=='DENDIF') THEN
IF(INDVAR>2.AND.INDVAR<9) CALL GXDENDIF ! density difference
ELSEIF(NPAT=='KESO') THEN ! Sources of turbulence energy and
!dissipation rate GXKESO
c.... "Un-comment" call to accutm('ke sor') in three places below, in
c order to compute and print cumulative time spent in GXKESO
c call accutm('ke sor',1)
CALL GXKESO(VIST,LEN1,VISL,FIXVAL)
IF(.NOT.LSG4.AND.ISG60==1.AND.ISG62==0)
1 CALL PB2_KESO(VIST,LEN1,VISL,FIXVAL)
c call accutm('ke sor',2)
c if(isweep==lsweep-1.and.istep==lstep)
c 1 call accutm('ke sor',1)
ELSEIF(NPAT=='KEBU') THEN
CALL GXKEGB
C
ELSEIF(NPAT=='RNGM') THEN ! Extra source of EP in RNG k-e model
IF(INDVAR==EP) CALL GXRNGM(VIST,LEN1,FIXFLU)
ELSEIF(NPAT=='KECH') THEN ! Extra source of EP in Chen-Kim k-e model
IF(INDVAR==EP.AND.ISC==16) THEN
L0G=L0F(GEN1); L0K=L0F(KE); L0VAL=L0F(VAL); L0VIST=L0F(VIST)
L0CO=L0F(CO)
DOIT=IENUTA==3.OR.IENUTA==4
DO I=1,NXNY
IF(SLD(I).OR.F(KAP+I)>FIXVAL) THEN
F(L0VAL+I)=0.0; F(L0CO+I)=0.0
ELSE
F(L0VAL+I)=0.25*F(L0G+I)*F(L0G+I)*F(L0VIST+I)*F(L0VIST+I)
1 /(F(L0K+I)+TINY)
IF(DOIT) F(L0VAL+I)=F(L0VAL+I)*F(L0FONE+I)
ENDIF
ENDDO
ENDIF
ELSEIF(NPAT=='REKE') THEN
IF(IENUTA==14) CALL GX_REALISABLE_KE
ELSEIF(NPAT=='KWSO') THEN
IF(IENUTA==10.OR.IENUTA==11) CALL GXKW_WILCOX
IF(IENUTA==15.OR.IENUTA==16) CALL GXKW_WILCOX
IF(IENUTA.GE.17.AND.IENUTA.LE.20) CALL GXKW_MENTER
ELSEIF(NPAT=='TSKE') THEN
IF(IENUTA==7) CALL GXTSKE(VIST,LEN1,FIXFLU)
ELSEIF(NPAT=='EYAP') THEN
IF(ILTLS>0) CALL GXEYAP
ELSEIF(NPAT=='BUOY') THEN ! Sources of momentum caused by buoyancy GXBUOY
CALL GXBUOY(BFC,CARTES,PARAB,DEN1,DEN2)
ELSEIF(NPAT=='STEN') THEN ! Sources of momentum caused by surface tension GXSURFT
CALL GXSURFT
ELSEIF(NPAT=='CLDA') THEN ! Conservative low-dispersion algorithm GXCLDA
CALL GXCLDA
ELSEIF(NPAT=='HOCS') THEN ! Higher-order schemes for convection GXHOCS
CALL GXHOCS
C.... Exit boundary condition for Supersonic flow Z-direction.
ELSEIF(NPAT=='OUTL') THEN
IF(INDVAR==P1) THEN
IF(ISC==13) THEN
INDW=ITWO(W1,LBNAME('WCRT'),.NOT.BFC)
CALL FN21(VAL,LOW(INDW),LOW(DEN1),0.0,-0.995)
ENDIF
ENDIF
C
C.... The following call to GXNEPA
C allows the "values" used in a COVAL command to be chosen as
C the "neighbour-values" of the cell in specified space, time or
C "variable-space" directions. The condition for the call is that
C the PATCH name should begin with the characters NE.
ELSEIF(NPAT(1:2)=='NE') THEN
CALL GXNEPA(NPAT)
C
C.... The following call to GXPOLR
C fixes the u- and v-velocities at the circumferential boundary of a
C cylindrical-polar domain for uniform flow at speed POLRA
ELSEIF(NPAT(2:4)=='POL') THEN
IF(.NOT.CARTES) CALL GXPOLR(NPAT(1:2))
C
C.... Sources used to set profiles GXPROF
ELSEIF(NPAT=='PROF') THEN
CALL GXPROF(DEN1,RHO1,FIINIT(H1))
C
C.... Calculate radiation sources . GXRADI
ELSEIF(NPAT=='RADI') THEN
CALL GXRADI(CARTES,TEMP1,DEN1,RHO1)
C
C.... Centrifugal & Coriolis forces due to rotation of
C coordinate system about axis of cylindrical-polar
C system. GXROTA
ELSEIF(NPAT=='ROTA' .AND. .NOT.GCV) THEN
IF(.NOT.CARTES.OR.BFC) CALL GXROTA(BFC,NX)
C
C.... Velocity resolutes at a curved inlet boundary GXBFC
ELSEIF(NPAT(1:3)=='BFC') THEN
CALL GXBFC
C
C.... Time-varying sources GXTIM
ELSEIF(NPAT(1:3)=='TIM') THEN
IF(NEZ(TIMA)) CALL GXTIM
C
C.... Sources provided for single-slab solution of fully-
C developed flow GXFDFS
ELSEIF(NPAT(1:3)=='FDF') THEN
CALL GXFDFS
C
C.... Chemical-reaction source GXCHSO
ELSEIF(NPAT=='CHSO') THEN
IF(.NOT.EXPLOS) THEN
CALL GXCHSO(TEMP1)
ENDIF
C
C.... Scalar fluctuations source GXFSQS
ELSEIF(NPAT=='FSQS') THEN
CALL GXFSQS
C.... Chemical-reaction sources, inlet properties and
C something else - from CHEMKIN GXCHKI
ELSEIF(NPATCH(1:5)=='CHEMK'.OR.NPATCH(1:4)=='NOCP'.OR.
1 NPATCH(1:6)=='CHEMID'.OR.NPATCH(1:5)=='CKWAL') THEN
IF(CHMKIN) CALL GXCHKI(0,0.0)
C.... Finite-rate chemical source terms for extended SCRS model
ELSEIF(ESCRS.AND.NPAT=='SCRS') THEN
CALL GXSCRI
!!! ELSEIF(LSG87) THEN
!!! CALL GXTCALC(0,0.0)
C.... Two-phase-flow patches
ELSEIF(.NOT.ONEPHS) THEN
C.... Momentum source caused by lateral gravity
C in 2-phase flow GXLATG
IF(NPAT=='LATG') THEN
CALL GXLATG
C.... Tentative length-scale source for two-fluid turbulence model
ELSEIF(NPAT=='LESO') THEN
CALL GXLESO('SOUTH')
C.... Shear source of lateral velocity for the 2-fluid
C turbulence model GXSHSO
ELSEIF(NPAT=='SHSO') THEN
CALL GXSHSO(INTFRC,LEN1)
C.... Sources providing for upwinding of the interphase
C transport terms GXIPST
ELSEIF(NPAT(1:4)=='IPST') THEN
CALL GXIPST(INTFRC,CINT(INDVAR),NPATCH)
C.... Sources of momentum for interfacial pressure GXINTP
ELSEIF(NPAT=='INTP') THEN
IF(NEZ(CPIP).AND.ISC==16) CALL GXINTP
C.... Sources of momentum for interfacial lift GXLIFT
ELSEIF(NPAT=='LIFT') THEN
IF(NEZ(CLIFT).AND.ISC==16) CALL GXLIFT
C.... Turbulence modulation due to presence of particles GXDISP
ELSEIF(NPAT=='KEDI') THEN
CALL GXDISP
ENDIF
ENDIF
IF(ISC==13.AND.INDVAR>=U1.AND.INDVAR<=W2) THEN
C.... Deduced inflow velocity at fixed-pressure outlet GXOUTLET
CALL GXOUTLET
ENDIF
C.... Angled inlet and outlet GXIO
IF(LGXIO) CALL GXIO
IF(LSUN) CALL GXSUN
C
IF(LSG62) CALL GXNOX
IF(LSG63) CALL GXNOX
c call accutm('grp 13',2)
c if(isweep==lsweep-1.and.istep==lstep) call accutm('grp 13',3)
GO TO 25
C***************************************************************
C
C
C--- GROUP 14. Downstream pressure for PARAB=.TRUE.
C
14 CONTINUE
GO TO 25
C***************************************************************
C * Make changes for these groups 15, 16, 17 & 18 only in group 19.
C***************************************************************
C
C
C--- GROUP 19. Special calls to GROUND from EARTH
C
19 GO TO (191,192,193,194,195,196,197,198,199,1910,1911),ISC
C
C * ----------GROUP 19 SECTION 1 ---- START OF TIME STEP.
191 CONTINUE
C.... Print or compute new grid and geometry
IF(MOVBFC) CALL GXBFGR
C.... Call to GXMOBS to modify blockages:
C * either to introduce complex initial field for steady cases;
C * or moving obstacles/blockages for unsteady cases; or both.
IF(IPORIB==1.OR.IPORIB==2) CALL GXMOBS
C.... Multiply DYG2D by EPOR, for thin-film lubrication
IF(ADJEPR) CALL FN26(DYG2D,EPOR)
C.... Call GXPIST to expand and contract the grid in accordance
C with the kinematics of crankshaft-connecting-rod mechanisms.
IF(NEZ(W1AD)) THEN
IF(NEZ(AZW1))
1 CALL GXPIST(ISTEP,NZ,TIM,DT,ISWEEP,LSWEEP,NPRINT,NTPRIN)
ENDIF
IF(RESET.AND..NOT.EXPLOS) THEN
CALL GXRSET
ENDIF
IF(.NOT.STEADY.AND.LSUN) CALL GXSUN
IF(.NOT.STEADY.AND.BLIN) CALL GXBLIN
C
GO TO 25
c * ----------group 19 section 2 ---- start of sweep.
192 CONTINUE
IF(SURFT) CALL GXSURFT
IF(MOVBFC) CALL GXBFGR
IF(CHMKIN) CALL GXCHKI(0,0.0)
IF(.NOT.DONE192.AND.IDISPA/=0.AND.ISWEEP==1) THEN ! dumping on start of sweep 1
I1=ITWO(3,2,PNAM(1:2)=='SW')
IF(PNAM(I1:I1+1)=='IN') THEN ! PNAM ends in IN so dump initial field
IF(STEADY.OR.(.NOT.STEADY.AND.(I1==3.OR.(I1==2.AND.
1 ISTEP==1)))) THEN
ISTP0=FSTEP ; ISWP0=FSWEEP ; FSTEP=0 ; FSWEEP=0
ISTP1=ISTEP ; ISWP1=ISWEEP ; ISTEP=0 ; ISWEEP=0
TIM0=TIM; TIM=0.0
CALL DUMPS(PNAM(1:1),GNAM(1:1),0,1,0,1)
FSTEP = ISTP0 ; FSWEEP=ISWP0 ; DONE192=.TRUE.
ISTEP = ISTP1 ; ISWEEP=ISWP1 ; TIM=TIM0
ENDIF
ENDIF
ENDIF
GO TO 25
C * ----------GROUP 19 SECTION 3 ---- START OF IZ SLAB.
C.... Modify porosities
193 IF(IPORIA/=0) THEN
CALL GXPORA(SETPOR)
IF(BFC) THEN
CALL BGEOM(1)
CALL BGEOM(2)
ENDIF
ENDIF
C
IF(LBMARK/=0) L0MARK= L0F(LBMARK)
GO TO 25
C
C * ------------------- SECTION 4 ---- Start of iterations over slab.
C
194 CONTINUE
IF(SURFT) CALL GXSURFT
IF(KEBUOY) CALL GXKEGB
C.... Set DOSKIN to true for WALL patches and correct G-function
IF(TURB) THEN
IF(NMWALL>0) CALL GXWFUN
ENDIF
C.... GXBLIN compute PASQUILL-F profiles for buoyancy forces
IF(PASQBUOY) CALL GXBLIN
IF(IENUTA==10.OR.IENUTA==11) CALL GXKW_WILCOX
IF(IENUTA==15.OR.IENUTA==16) CALL GXKW_WILCOX
IF(IENUTA.GE.17.AND.IENUTA.LE.20) CALL GXKW_MENTER
Call realisable k-e model after upgrading of G-function in E2conn for walls
IF(IENUTA.EQ.14) CALL GX_REALISABLE_KE
C.... Calculate phase-average z-direction velocities
C Copy coding if velocities for other directions are required
IF(.NOT.ONEPHS) THEN
IF(LBWAVE>0) THEN
L0WAVE=L0F(LBWAVE)
L0R1=L0F(9); L0R2=L0F(10)
L0W1=L0F(7); L0W2=L0F(8)
RH1=1.0/RHO1A; RH2=1.0/RHO2A
DO 1940 I=1,NXNY
F(L0WAVE+I)=(F(L0W1+I)*F(L0R1+I)*RH1 +
1 F(L0W2+I)*F(L0R2+I)*RH2) /
1 (F(L0R1+I)*RH1 + F(L0R2+I)*RH2)
1940 CONTINUE
if(dbgrnd) call prn('wave',lbwave)
ENDIF
C
IF(STORE(P2)) THEN
IF(STORE(VPOR)) THEN
C.... Extra second-phase pressure proportional to volume fraction
C divided by porosity.
CALL FN15(P2,R2,VPOR,0.0,PHS2PA)
ELSE
IF(PHS2P) THEN
C.... Extra second-phase pressure proportional to:
C a*(r2-r2crit)**b if r2 gt r2crit.
L0R2=L0F(R2)
L0P2=L0F(P2)
RLX=ABS(DTFALS(P2))
DO 1941 I=1,NXNY
IF(F(L0R2+I)>PHS2PA) THEN
GP2=PHS2PB*(F(L0R2+I)-PHS2PA)**PHS2PC
ELSE
GP2=0.0
ENDIF
F(L0P2+I)=RLX*GP2+(1.0-RLX)*F(L0P2+I)
1941 CONTINUE
ELSE
C.... Extra second-phase pressure proportional to volume fraction
CALL FN2(P2,R2,0.0,PHS2PA)
ENDIF
ENDIF
ENDIF
ENDIF
IF(LSG62) CALL GXNOX
IF(LSG63) CALL GXNOX
C.... Calculate derivatives of phase 1 temperature
IF(ITEM1>0.OR.TEMP1>0) THEN
LBSCAL= ITWO(ITEM1,TEMP1,ITEM1>0); L0SCAL=L0F(LBSCAL)
IF(ISDTDX.AND.NX>1) THEN
L0DTDX= L0F(LBDTDX)
CALL FN1(-L0DTDX,0.0)
CALL FNDFDX(L0SCAL,L0DTDX,IZSTEP,.TRUE.)
ENDIF
IF(ISDTDY.AND.NY>1) THEN
L0DTDY= L0F(LBDTDY)
CALL FN1(-L0DTDY,0.0)
CALL FNDFDY(L0SCAL,L0DTDY,IZSTEP,.TRUE.)
ENDIF
IF(ISDTDZ.AND.NZ/=1.AND..NOT.PARAB) THEN
L0DTDZ= L0F(LBDTDZ)
CALL FN1(-L0DTDZ,0.0)
CALL FNDFDZ(L0SCAL,L0DTDZ,IZSTEP,.TRUE.)
ENDIF
ENDIF
C.... Calculate derivatives of phase 2 temperature
IF(ITEM2>0.OR.TEMP2>0) THEN
LBSCAL= ITWO(ITEM2,TEMP2,ITEM2>0); L0SCAL=L0F(LBSCAL)
IF(ISDTX2.AND.NX>1) THEN
L0DTDX= L0F(LBDTX2)
CALL FN1(L0DTDX,0.0)
CALL FNDFDX(L0SCAL,L0DTDX,IZSTEP,.TRUE.)
ENDIF
IF(ISDTY2.AND.NY>1) THEN
L0DTDY= L0F(LBDTY2)
CALL FN1(L0DTDY,0.0)
CALL FNDFDY(L0SCAL,L0DTDY,IZSTEP,.TRUE.)
ENDIF
IF(ISDTZ2.AND.NZ/=1.AND..NOT.PARAB) THEN
L0DTDZ= L0F(LBDTZ2)
CALL FN1(L0DTDZ,0.0)
CALL FNDFDZ(L0SCAL,L0DTDZ,IZSTEP,.TRUE.)
ENDIF
ENDIF
GO TO 25
1911 CONTINUE
C * ------------------- SECTION 11---- After calculation of convection
C fluxes for scalars, and of volume
C fractions, but before calculation of
C scalars or velocities
IF(.NOT.ONEPHS) THEN
IF(NEZ(CPIP)) CALL GXINTP
IF(NEZ(CLIFT)) CALL GXLIFT
ENDIF
GO TO 25
199 CONTINUE
C * ------------------- SECTION 9 ---- Start of solution sequence for
C a variable
IF(FSQSOR.AND.INDVAR==LBNAME('FSQ')) CALL GXFSQS
GO TO 25
1910 CONTINUE
C * ------------------- SECTION 10---- Finish of solution sequence for
C a variable
IF(CHMKIN) CALL GXCHKI(0,0.0)
GO TO 25
C
C * ------------------- SECTION 5 ---- Finish of iterations over slab.
195 IF(POTVEL) THEN ! Compute velocities from potential differences.
IF(ISWEEP/=1) THEN
CALL GXPOTV(EPOR,NPOR,HPOR,NZ,DZG,BFC)
IF(POTCMP.AND.ISWEEP>2)
1 CALL GXPOTC(EPOR,NPOR,HPOR,PORIA,PORIB,NZ) ! allow for
ENDIF ! compressibility flow.
ENDIF
IF(PSIVEL) CALL GXPSIV(BFC) ! deduce velocities from stream function psi
IF(FDFSOL) CALL GXFDFS ! fully-developed flow
IF(CHMKIN) CALL GXCHKI(0,0.0) ! chemkin
C.... Turbulence-energy generation rates
IF(LBGENK/=0) CALL FN21(LBGENK,GEN1,VIST,0.0,1.0) ! phase 1
IF(LBGNK2/=0) CALL FN21(LBGNK2,GEN2,VIST,0.0,1.0) ! phase 2
IF(LBVLSQ/=0) CALL FNVLSQ(LBVLSQ,1) ! sum of squares of phase1 velocity
IF(LBVABS/=0) THEN ! absolute velocity, i.e. SQRT(VLSQ)
IF(LBVLSQ/=0) THEN
CALL FN30A(LBVABS,LBVLSQ) ! VABS is SQRT(VLSQ)
ELSE
CALL FNVLSQ(LBVABS,1) ! put VLSQ in VABS, then take square root
CALL FN30A(LBVABS,LBVABS) ! VABS is SQRT(VLSQ)
ENDIF
IF(ISCREY) THEN
L0VOL =L0F(VOL) ; L0VIST=L0F(VIST) ; L0VISL=L0F(VISL)
L0CREY=L0F(LBCREY) ; L0VABS=L0F(LBVABS)
DO I=1,NXNY
IF(.NOT.SLD(I)) F(L0CREY+I)=F(L0VABS+I)*F(L0VOL+I)**0.3333/
1 (F(L0VIST+I)+F(L0VISL+I))
ENDDO
ENDIF
ENDIF
IF(LBVAB2/=0) THEN ! absolute velocity, phase 2
CALL FNVLSQ(LBVAB2,2) ! put VLSQ in VABS, then take square root
CALL FN30A(LBVAB2,LBVAB2) ! VABS is SQRT(VLSQ)
ENDIF
IF(ISTREY) THEN
L0VIST=L0F(VIST); L0VISL=L0F(VISL); L0TREY=L0F(LBTREY)
DO I=1,NXNY
IF(.NOT.SLD(I)) F(L0TREY+I) = F(L0VIST+I)/F(L0VISL+I)
ENDDO
ENDIF
IF(ESCRS) CALL GXSCRI
C....Mach number calculation for Compressible Flow; it must be
c done here as parab=t field printout precedes gr 19 sc 6
IF(PARAB) THEN
IF(LBMACH/=0.OR.LBPTOT/=0)
1 CALL GXMACH(LBMACH,LBPTOT,LBVABS,0.,1)
IF(.NOT.ONEPHS) THEN
IF(LBMAC2/=0.OR.LBPTO2/=0)
1 CALL GXMACH(LBMAC2,LBPTO2,LBVAB2,0.,2)
IF(SLIPVL) THEN ! conditionally compute the slip velocities
IF(LBSLU/=0) CALL FN10(LBSLU,U1,U2,0.0,1.0,-1.0)
IF(LBSLV/=0) CALL FN10(LBSLV,V1,V2,0.0,1.0,-1.0)
IF(LBSLW/=0) CALL FN10(LBSLW,W1,W2,0.0,1.0,-1.0)
ENDIF
ENDIF
ENDIF
GO TO 25
C
C GROUP 19 SECTION 6 Finish of iz slab
196 IF(PARAB) THEN
IF(IDISPA>0) CALL GXPARA
ELSEIF(ARRIVT) THEN
IF(ISWEEP==LSWEEP) CALL GXARRT
ENDIF
C....Mach number calculation for Compressible Flow
IF(.NOT.PARAB) THEN
IF(LBMACH/=0.OR.LBPTOT/=0)
1 CALL GXMACH(LBMACH,LBPTOT,LBVABS,0.,1)
IF(.NOT.ONEPHS) THEN
IF(LBMAC2/=0.OR.LBPTO2/=0)
1 CALL GXMACH(LBMAC2,LBPTO2,LBVAB2,0.,2)
C.... For two-phase flow, conditionally compute the slip velocities
IF(SLIPVL) THEN
IF(LBSLU/=0) CALL FN10(LBSLU,U1,U2,0.0,1.0,-1.0)
IF(LBSLV/=0) CALL FN10(LBSLV,V1,V2,0.0,1.0,-1.0)
IF(LBSLW/=0) CALL FN10(LBSLW,W1,W2,0.0,1.0,-1.0)
ENDIF
ENDIF
ENDIF
IF(LSG62) CALL GXNOX
IF(LSG63) CALL GXNOX
IF(BLIN) CALL GXBLIN ! Boundary-layer inlet profile facility
IF(IZ0) CALL DMPDSP
ENDIF
GO TO 25
C
C GROUP 19 SECTION 7 FINISH OF SWEEP.
197 CONTINUE
C
C.... Sequence permitting intervention to terminate execution
C when a file named STOPJOB exists, for example created by
C editing in a second window
C Note that the name stopjob will not suffice on case-sensitive
C systems such as UNIX
C Note that this may not work on all computer systems, in which
C case the coding from here to STOPJOB END below should be
C deactivated.
IERR=-1
ISTOP=5
IF(MOD(ISWEEP,ISTOP)==0) THEN
IF(MIMD.AND.NPROC>1) THEN
IF(MYID==0) THEN
IF(LUNIT(26)/=0) CALL OPENZZ(26,IERR)
ENDIF
c... pass first proc IERR value to others
CALL RTFSII(IERR)
ELSE
IF(LUNIT(26)/=0) CALL OPENZZ(26,IERR)
ENDIF
ENDIF
IF(IERR==0) THEN
ENUFSW=.TRUE.
LSTEP=ISTEP
LSWEEP=ISWEEP+1
LUNIT(26)=0
CALL WRIT40('file stopjob found; execution stopped ')
CALL WRYT40('file stopjob found; execution will be ')
CALL WRYT40('stopped after completion of print-out ')
ENDIF
c.... stopjob end
IF(CHMKIN) CALL GXCHKI(0,0.0)
C.... Print-out of Y+, stress and Stanton number at WALL-type patches
IF(TURB) THEN
IF(NMWALL>0) THEN
IF(.NOT.NULLPR.AND.(YPLS .OR. WALPRN)) CALL GXYPLS
ENDIF
ENDIF
c
C.... "Time-out" sequence activated by setting the PIL variable ISG20
C equal to the maximum number of seconds which the run may take.
IF(ISG20>0) THEN
CALL SECONDS(NUMSEC)
IF(NUMSEC>=ISG20) THEN
CALL WRITST
CALL WRIT40('Execution "timed out" at numsec seconds')
CALL WRIT1I(' numsec',ISG20)
CALL WRITST
ENUFSW=.TRUE. ; LSTEP=ISTEP
IF(ISWEEP/=LSWEEP) LSWEEP=ISWEEP+1
ENDIF
ENDIF
c
c.... print-out elicited by SPEDAT('PRINT', ... statements in Q1
IF(ISWEEP==FSWEEP) THEN
IF(ISTEP==1) THEN
c.... get number of special-print commands
NUMSPPR=0
CALL GETSDI('PRINT','NUMBER',NUMSPPR)
ENDIF
ENDIF
IF(TIMPRN.AND.SWPRIN.AND..NOT.NULLPR) THEN
C
IF(LBVOXY>0) L0VOXY=M0F1(LBVOXY)
IF(LBVOYZ>0) L0VOYZ=M0F1(LBVOYZ)
IF(LBVOZX>0) L0VOZX=M0F1(LBVOZX)
C
IF(LBVOXY>0.OR.LBVOYZ>0.OR.LBVOZX>0) THEN
DO IZZ = 1,NZ
DO IX=1,NX
DO IY=1,NY
I3D=IY+(IX-1)*NY+(IZZ-1)*NFM
IF(LBVOXY/=0) F(L0VOXY+I3D)=VORTXY(1,IX,IY,IZZ)
IF(LBVOYZ/=0) F(L0VOYZ+I3D)=VORTYZ(1,IX,IY,IZZ)
IF(LBVOZX/=0) F(L0VOZX+I3D)=VORTZX(1,IX,IY,IZZ)
ENDDO
ENDDO
ENDDO
ENDIF
ENDIF
IF(NUMSPPR>0) THEN
IF(TIMPRN.AND.SWPRIN.AND..NOT.NULLPR) THEN
call writ40('Special print-out via spedat commands')
c.... how many commands?
DO I=1,NUMSPPR
c.... get the command
CALL GETSDC('PRINT','COMMAND'//DIGIT(I),COMMAND)
c.... carry out the commands
IF(COMMAND(1:3)=='POW') THEN
c.... do this before the istep test so that result is included in phi
CALL GETSDR('PRINT','POWER',POWER)
IF(MOD(ISTEP,NTPRIN)==0) THEN
CALL WRIT40(COMMAND(10:13)//
1 'holds "power-lawed" values of '
1 //COMMAND(5:8))
CALL WRIT1R('power ',POWER)
ENDIF
CALL PUT_POWER( COMMAND(5:8) , COMMAND(10:13),POWER )
CALL WRITBL
ELSEIF(MOD(ISTEP,NTPRIN)==0) THEN
CALL WRIT40('Spedat print command is '//command)
IF(COMMAND(1:3)=='AVE') THEN
CALL AVERAGE(AVE,COMMAND(5:8))
CALL WRIT40
1 ('mass-average value of '//COMMAND(5:8))
CALL WRIT1R('average ',AVE)
CALL WRITBL
ELSEIF(COMMAND(1:5)=='TOTAL') THEN
CALL SUMALL(SUM,COMMAND(7:10))
CALL WRIT40('total whole-domain amount of '//
1 COMMAND(7:10))
CALL WRIT1R('TOTAL ',SUM)
CALL WRITBL
ELSEIF(COMMAND(1:6)=='MINMAX') THEN
CALL WRIT40('hilo3d called for '//COMMAND(8:11))
CALL HILO3D(LBNAME(COMMAND(8:11)))
CALL WRIT2R('MinValue',RLOW3D,'MaxValue',HI3D)
CALL WRITBL
ELSEIF(COMMAND(1:5)=='VALUE') THEN
CALL GETSDI('PRINT','IXLOC',IXLOC)
CALL GETSDI('PRINT','IYLOC',IYLOC)
CALL GETSDI('PRINT','IZLOC',IZLOC)
CALL GETVAL(COMMAND(7:10),IXLOC,IYLOC,IZLOC)
ENDIF
ENDIF
ENDDO
ENDIF
ENDIF
C.... Dump fields after each IDISPA sweeps, eg for PHOTON viewing
IF(IDISPA/=0) THEN
IF(PNAM(1:2)=='SW') THEN
IF(STEADY) THEN
CALL DUMPS(PNAM(1:1),GNAM(1:1),ISWEEP,IDISPA,IDISPB,IDISPC)
ELSE
STEADY=.TRUE.
CALL DUMPS(PNAM(1:1),GNAM(1:1),ISWEEP,IDISPA,1,LSWEEP)
STEADY=.FALSE.
DONE192=.FALSE.
ENDIF
IF(CALFOR) THEN
IDSPB=ITWO(1,IDISPB,IDISPB==0)
IDSPC=ITWO(LSWEEP,IDISPC,IDISPC==0)
IDSPA=ITWO(1,IDISPA,IDISPA==0)
IF(ISWEEP>=IDSPB.AND.ISWEEP<=IDSPC.AND.
1 MOD(ISWEEP,IDSPA)==0) THEN
CALL FORCES(1) ! Integrate
CALL FORCES(3) ! Output to table file
ENDIF
ENDIF
ENDIF
ENDIF
C.... Angled inlet and outlet GXIO
IF(LGXIO) CALL GXIO
IF(BLIN) CALL GXBLIN ! Boundary-layer inlet profile facility
GO TO 25
C
C GROUP 19 SECTION 8 FINISH OF TIME STEP.
198 CONTINUE
C.... Save fields and cell-corner coordinates to a series of files,
C for examination via VR-Viewer or PHOTON
C.... PNAM and GNAM are the names of the field and (if BFC)
C grid files, and IDISPA is the frequency of dumping.
C They are EQUIVALENCEd in the included Common file GRDLOC to the
C PIL variables CSG1 and CSG2, and can therefore be set in the
C satellite, eq by the line:
C CSG1=PHI;CSG2=XYZ;IDISPA=1
C which will create files p1da or p1, p2da or p2, etc, & x1da or x1,
C etc according to whether phida and xyzda are true (the default)
C or false in the PREFIX file
C
!!! IF(LSG87) CALL GXTCALC(0,0.0)
IF(PNAM(1:2)/='SW'.AND.IDISPA/=0) THEN
LGTEMP=BFC ; BFC=BFC.OR.NOTBFC
CALL DUMPS(PNAM(1:1),GNAM(1:1),ISTEP,IDISPA,IDISPB,IDISPC)
BFC=LGTEMP
IF(CALFOR) THEN
IDSPB=ITWO(1,IDISPB,IDISPB==0)
IDSPC=ITWO(LSTEP,IDISPC,IDISPC==0)
IDSPA=ITWO(1,IDISPA,IDISPA==0)
IF(ISTEP>=IDSPB.AND.ISTEP<=IDSPC.AND.MOD(ISTEP,IDSPA)==0) THEN
CALL FORCES(1) ! Integrate
CALL FORCES(3) ! Output to table file
ENDIF
ENDIF
C... Save monitor plot image for each time step
GXMON=(TSTSWP/=0)
CALL GETSDL('GXMONI','TRANSIENT',GXMON)
IF(GXMON) THEN
IFST=ITWO(FSTEP,IDISPB,IDISPB==0) ! first step to dump
ILST=ITWO(LSTEP,IDISPC,IDISPC==0) ! last step to dump
GXMON=ISTEP>=IFST.AND.ISTEP<=ILST.AND.
1 MOD(ISTEP,IDISPA)==0
ENDIF
IF(GXMON) THEN
IF(PLTCLSX) THEN
COMMAND='gxmoni'
IF (ISTEP>=100000) THEN
WRITE (COMMAND(7:12),FMT='(BN,I6)') ISTEP
ELSEIF(ISTEP>=10000) THEN
WRITE (COMMAND(7:11),FMT='(BN,I5)') ISTEP
ELSEIF(ISTEP>=1000) THEN
WRITE (COMMAND(7:10),FMT='(BN,I4)') ISTEP
ELSEIF(ISTEP>=100) THEN
WRITE (COMMAND(7:9),FMT='(BN,I3)') ISTEP
ELSEIF(ISTEP>=10) THEN
WRITE (COMMAND(7:8),FMT='(BN,I2)') ISTEP
ELSE
WRITE (COMMAND(7:7),FMT='(BN,I1)') ISTEP
ENDIF
LL=LENGZZ(COMMAND)
CALL DUMPZZ(COMMAND(1:LL),LL)
ELSE
!DEC$ IF .NOT.DEFINED(LINUX) .AND. DEFINED(_X64_)
!DEC$ ATTRIBUTES C, ALIAS : 'AllStepsONC' :: SDASCON
!cs interface code
!! IF(ISTEP==IFST.OR.ISTEP==ILST) THEN
CALL SDASCON() ! toggle dumping on
!------------------------------------------------------
!this block is necessary for the super fast cases
!it pauses earth immediately after sending signal to take screenshots
!otherwise the data in the gui may skip a few sweeps before it starts the screenshot r
!corresponding unpause call in abstractions.screenShotAll(special==1)
!DEC$ ATTRIBUTES C, ALIAS : 'earTogglePause' :: ETP
CALL ETP()
!------------------------------------------------------
!! ENDIF
!DEC$ ENDIF
ENDIF
ELSE
!DEC$ IF .NOT.DEFINED(LINUX) .AND. DEFINED(_X64_)
IF(.NOT.PLTCLSX) THEN
!! IF(ISTEP==IFST.OR.ISTEP==ILST) THEN
!DEC$ ATTRIBUTES C, ALIAS : 'dumpAllStepsOFFC' :: SDASCOFF
CALL SDASCOFF() ! toggle dumping off
!! ENDIF
ENDIF
!DEC$ ENDIF
ENDIF
ENDIF
C
IF(MOVBFC) THEN
C.... Print-out versus time by use of GRAPH
TIME(ISTEP)=TIM
C.... set ordinate values
OR(ISTEP,1)=VARONE(LBNAME('VOLU'),1,1)
OR(ISTEP,2)=VARONE(P1,1,1)
IF(ISTEP/=LSTEP) GO TO 25
C.... call GRAPH
C.... set ordinate labels
LABOR(1)='VOLU' ; LABOR(2)='PRES'
IOSC(1)=1 ; IOSC(2)=2 ; IOSC(3)=3 ; IOSC(4)=4
NABDIM=100 ; NORDIM=10 ; NOR=2
ABSIZ=0.6 ; ORSIZ=0.4
C
CALL GRAPH(TIME,NABDIM,ISTEP,'TIME',OR,NORDIM,NOR,
1 LABOR,IOSC,ABSIZ,ORSIZ,3,LUPR1)
C.... end of GRAPH sequence
ENDIF
C
IF(.NOT.STEADY) THEN
IF(.NOT.PARAB.AND.IDISPA>0.AND.PNAM==' ') CALL GXPARA
IF(IDISPD>0) CALL DMPDSP
ENDIF
GO TO 25
C***************************************************************
C
C
C--- GROUP 20. Preliminary print-out
C
20 IF(ECHO) THEN ! echo all 25 data groups to RESULT file
CALL DATPRN(Y,Y,Y,Y, Y,Y,Y,Y, Y,Y,Y,Y, Y,Y,Y,Y,
1 Y,Y,Y,Y, Y,Y,Y,Y, Y)
ELSE ! echo only data group 1 to RESULT
CALL DATPRN(Y,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N)
ENDIF
GO TO 25
C***************************************************************
C
C--- GROUP 21. Special print-out to screen
21 CONTINUE
GO TO 25
C***************************************************************
C * Make changes for these groups 22 and 23 only in group 19.
C***************************************************************
C
C
C--- GROUP 24. Dumps for restarts
C
24 CONTINUE
C * Insert CALL DUMP where required in group 19.
C
25 CONTINUE
if(test) call fcheck(25,'grex3 ')
C***************************************************************
C
IF(SPCLGR) THEN
C.......................................... special GROUNDS
C.... S2SR, for surface-to-surface radiation
IF(S2SR) THEN
RADCVD = .FALSE.
CALL GETSPD('CVD','RADCVD',3,RDUM,IDUM,RADCVD,' ',IERR)
C.... Only call GXS2SR GR.9, SEC 7 when INDVAR is temperature
IF(INDVAR==ITEM1.OR..NOT.(IGR==9.AND.ISC==7)) THEN
IF(RADCVD) THEN
CALL GXS2SR_CVD
ELSE
CALL WRIT40('Please set:')
CALL WRIT40('SPEDAT(SET,CVD,RADCVD,L,T)')
CALL WRIT40('in Q1')
CALL SET_ERR(587,'Old S2SR not supported',1)
ENDIF
ENDIF
ENDIF
IF(NAMGRD=='NONE') THEN
C....NONE
ELSEIF(NAMGRD=='FURN') THEN
C....FURN
CALL FURNGR
C.... ESTER, the aluminium-smelter module
ELSEIF(NAMGRD=='ESTR') THEN
CALL ESTRGR
C....HOTBOX
ELSEIF(NAMGRD=='HTBX') THEN
CALL HTBXGR
C....TACT
ELSEIF(NAMGRD=='TACT') THEN
CALL TACTGR
C....COMPRESSIBLE
ELSEIF(NAMGRD=='CONV') THEN
CALL GXCONV
C....CVD
ELSEIF(NAMGRD=='CVD ') THEN
CALL GXCVD
C....MICA special GROUND's
ELSEIF(NAMGRD=='MICA') THEN
CALL MICAGR
C... F1 - calculates drag and lift for formula-1 car
ELSEIF(NAMGRD=='F1') THEN
CALL F1GRD
C... Water cooled steam condenser WATSTCON
ELSEIF(NAMGRD=='WATS') THEN
CALL GXCOND
ELSEIF(NAMGRD=='WAVE') THEN
CALL GXWAVE
ENDIF
ENDIF
if(test) call fcheck(26,'grex3 ')
NAMSUB='grex3'
if(dbgloc) then
call banner(2,namsub,120121)
endif
END
c