c
c
C To activate the coding set S2SR=T in the Q1 file.
C
C Limitations: Note that the coding assumes that the radiative
C exchange coefficients provided by the user take into account
C the emissivities, the Stefan-Boltzmann constant, but not the
C zone area.
C
C Click here
c for full information
C-----------------------------------------------------------------------
SUBROUTINE GXS2SR
INCLUDE 'patcmn'
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'grdbfc'
INCLUDE 'parear'
PARAMETER (NLG=100, NIG=100, NRG=200, NCG=100)
C.... Set dimension of patch-name array here. WARNING: the array
C NAMPAT in the MAIN program of the satellite, and EARTH must have
C the same dimension.
PARAMETER (NPNAM=200)
C
C.... To be able to access patchwise coefficients and values for qdot bc.
COMMON/ICOVL/M04,I0PHI
C
COMMON/NAMFN/NAMFUN,NAMSUB
COMMON/LGRND/LG(NLG)/IGRND/IG(NIG)/RGRND/RG(NRG)/CGRND/CG(NCG)
COMMON/LUNITS/LUNIT(60)/INFXY/INFO(20)
COMMON/GENL/LGFL1(43),RESTRT,LGFL45(16)
COMMON/GENI/NXNY,IGFL1(44),LOOPZ,IGFL47(14)
LOGICAL LGFL1,RESTRT,LGFL45
C
C.... Common LUSF for zero location of F-array
COMMON/LUSF/L0RPNO,L0AVRF,L0AVTS,L0AVSH,L0RC,L0CC,L0CCF,
1 L0QFLX,L0RFAC,L0COND,L0WPNO,L0CP1,L0HTC,L0AE,L0AN,L0AS,L0AH,
1 L0EXRD,L0EXLN,L0QEXT,L0XCOF,L0KDM,L0NTS
C
LOGICAL LG,LFOUND,LTURB,LINIL,GTZ,LTZ,EQZ,QEQ,NEZ,LEXR
LOGICAL LDUM, NOTEXISTZONE, EXISTSPATCH
CHARACTER*4 CG
CHARACTER*6 NAMFUN,NAMSUB
C
SAVE NTZ,JTEM1,JPRPS,IVISIT,LURFAC,LUPHI,IGRP97,IGP914,
1 LTURB,LINIL,LEXR
C
NAMSUB='GXS2SR'
C-----------------------------------------------------------------------
C.... Group 1 Section 1, for opening data files |
C-----------------------------------------------------------------------
IF(IGR.EQ.1.AND.ISC.EQ.1) THEN
IF(.NOT.NULLPR) THEN
IF(MIMD.AND.NPROC.GT.1)THEN
IF(MYID.EQ.0)THEN
CALL WRYT40('GROUND file is GXS2SR.FOR of: 150416 ')
CALL WRIT40('GROUND file is GXS2SR.FOR of: 150416 ')
ENDIF
ELSE
CALL WRYT40('GROUND file is GXS2SR.FOR of: 150416 ')
CALL WRIT40('GROUND file is GXS2SR.FOR of: 150416 ')
ENDIF
ENDIF
C.... Initialise IVISIT
IVISIT=0
C.... LINIL determines whether the users part of the f-array needs
C to be initialised
LINIL=.TRUE.
C.... Assign units for file openings
LUPHI=LUNIT(23)
IF(LUPHI.EQ.0) LUPHI=LUNIT(22)
LURFAC=33
C.... Open data file containing internal radiative exchange coefs
C INQUIRE if RADI.DAT exists in the local directory else CALL
C GX2DVF
IF(MIMD.AND.NPROC.GT.1)THEN
IF(MYID.EQ.0)INQUIRE(FILE='RADI.DAT',EXIST=LEXR)
CALL RTFALI(LEXR,1)
IF(LEXR)THEN
IF(MYID.EQ.0) OPEN(LURFAC,FILE='RADI.DAT',STATUS='OLD',
1 ACCESS='SEQUENTIAL',FORM='FORMATTED')
ELSE
CALL SET_ERR(560,
* 'Error. File RADI.DAT is absent at parallel run.',1)
RETURN
ENDIF
ELSE
INQUIRE(FILE='RADI.DAT',EXIST=LEXR)
IF(LEXR) OPEN(LURFAC,FILE='RADI.DAT',STATUS='OLD',
1 ACCESS='SEQUENTIAL',FORM='FORMATTED')
ENDIF
C
C-----------------------------------------------------------------------
C.... Group 1 Section 2, for PRELIMINARIES such as define zero |
C location of user part of F-array, Reading in the exchange |
C coefficients, storing the PHOENICS PATCH numbers of each |
C THERMAL zone PATCH. |
C-----------------------------------------------------------------------
C.... Perform preliminaries only on the first visit
IF(IVISIT.EQ.0) THEN
IVISIT=1
C.... Determine if flow is turbulent
LTURB=.FALSE.
IF(NEZ(ENUT)) LTURB=.TRUE.
C.... Read in header and number of thermal zones
IF(LEXR) THEN
IF(MIMD.AND.NPROC.GT.1)THEN
IF(MYID.EQ.0)THEN
READ(LURFAC,*)
READ(LURFAC,*) NTZ
ENDIF
CALL RTFAII(NTZ,1)
ELSE
READ(LURFAC,*)
READ(LURFAC,*) NTZ
ENDIF
ELSE
C.... Loop through all patches to determine patch numbers for radiative
C patches. At present assumes that one patch for each thermal zone
C and that each patch commences with '@RI'
NTZ=0
DO 10 IPN=1,NUMPAT
IF(NAMPAT(IPN)(1:3).EQ.'@RI') NTZ=NTZ + 1
10 CONTINUE
ENDIF
C
C.... Define zero locations for the relevant 1-d stores
C L0RPNO for PHOENICS PATCH numbers of all radiative patches
C L0AVRF for mean radiative flux W/m2
C L0AVTS for mean surface radiation temperatures K
C L0AVSH for total source of heat flux, conductive and convective W/m2
C L0RC for radiative coefficients W/m2/K
C L0CC for conductive coefficients W/m2/K
C L0CCF for heat transfer coefficients for fluid side W/m2/K
C L0QFLX for value of prescribed heat flux zone conditions W/m2
C L0RFAC for radiative exchange coefficient
C L0COND for laminar thermal conductivity W/m2
C L0EXRD for external radiative exchange factor and temperature
C L0EXLN for external linear heat transfer coefficient and temperature
C L0QEXT for sum of external linear and radiative heat flux
C L0KDM for conductivity divided by metal thickness for thin plates
C L0NTS for node number of adjacent thermal node of thin plates
C L0XCOF for sum of external linear and radiative coefficient
C L0WPNO for PHOENICS WALL PATCH numbers corrsponding to the
C radiative patches
C L0CP1 for specific heat at constant pressure J/kg/K
C L0HTC for heat transfer coefficient when turbulent flow
C
CALL GXMAKE(L0RPNO,NTZ,'RPNO')
CALL GXMAKE(L0AVRF,NTZ,'AVRF')
CALL GXMAKE(L0AVTS,NTZ,'AVTS')
CALL GXMAKE(L0AVSH,NTZ,'AVSH')
CALL GXMAKE(L0RC ,NTZ,'RC ')
CALL GXMAKE(L0CC ,NTZ,'CC ')
CALL GXMAKE(L0CCF ,NTZ,'CCF ')
CALL GXMAKE(L0QFLX,NTZ,'QFLX')
CALL GXMAKE(L0RFAC,NTZ*NTZ,'RFAC')
CALL GXMAKE(L0COND,NXNY*NZ,'COND')
CALL GXMAKE(L0EXRD,2*NTZ,'EXRD')
CALL GXMAKE(L0EXLN,2*NTZ,'EXLN')
CALL GXMAKE(L0QEXT,NTZ,'QEXT')
CALL GXMAKE(L0XCOF,NTZ,'XCOF')
CALL GXMAKE(L0KDM ,NTZ,'KDM ')
CALL GXMAKE(L0NTS ,NTZ,'NTS ')
IF(LTURB) THEN
CALL GXMAKE(L0WPNO,2*NTZ,'WPNO')
CALL GXMAKE(L0CP1 ,NXNY*NZ,'CP1 ')
CALL GXMAKE(L0HTC ,NXNY*NZ,'HTC ')
CALL ZERNUM(L0RPNO,L0HTC+NXNY*NZ - L0RPNO)
ELSE
CALL ZERNUM(L0RPNO,L0NTS+NTZ - L0RPNO)
ENDIF
C
C.... Loop through all patches to obtain patch numbers for radiative
C patches. At present it is assumed that there is one patch for
C each thermal zone .
C
IR=0
DO 11 IPN=1,NUMPAT
C.... It is assumed that internal radiative thermal zones have
C PATCH name commencing with '@RI' for internal radiation.
IF(NAMPAT(IPN)(1:3).EQ.'@RI') THEN
IR = IR+1
F(L0RPNO+IR) = FLOAT(IPN)
ENDIF
11 CONTINUE
C
C.... Check that number of patches found for thermal zones are equal
C to the number of thermal zones read in from files. Note that
C this is not a thorough check when LEXR is .FALSE..
C
IF(IR.NE.NTZ) THEN
IF(MIMD.AND.NPROC.GT.1)THEN
IF(MYID.EQ.0)THEN
CALL WRITST
CALL WRITBL
WRITE(LUPR1,*) ' user set number of radiative zone = ',
1 NTZ
WRITE(LUPR1,*) 'Number of PATCHES found in Q1 for '
1 //'radiative zone'
WRITE(LUPR1,*) 'is insufficient. Number of PATCHES'
1 //' found = ',IR
WRITE(LUPR1,*) 'Specify remaining PATCHES '
ENDIF
ELSE
CALL WRITST
CALL WRITBL
WRITE(LUPR1,*) ' user set number of radiative zone = ',NTZ
WRITE(LUPR1,*) 'Number of PATCHES found in Q1 for'
1 //' radiative zone'
WRITE(LUPR1,*) 'is insufficient. Number of PATCHES'
1 //' found = ',IR
WRITE(LUPR1,*) 'Specify remaining PATCHES '
ENDIF
CALL SET_ERR(298, 'Error. See result file.',1)
ENDIF
ENDIF
ELSEIF(IGR.EQ.1.AND.ISC.EQ.2) THEN
C
C.... Read in, or calculate, Radiative exchange coefficient
C
IF(LEXR) THEN
IF(.NOT.NULLPR) THEN
IF(MIMD.AND.NPROC.GT.1)THEN
IF(MYID.EQ.0)THEN
CALL WRITBL
CALL WRITST
CALL WRITBL
WRITE(LUPR1,*) 'Viewfactors read from RADI.DAT '
CALL WRYT40('Viewfactors read from RADI.DAT ')
CALL WRITBL
CALL WRITST
ENDIF
ELSE
CALL WRITBL
CALL WRITST
CALL WRITBL
WRITE(LUPR1,*) 'Viewfactors read from RADI.DAT '
CALL WRYT40('Viewfactors read from RADI.DAT ')
CALL WRITBL
CALL WRITST
ENDIF
ENDIF
IF(MIMD.AND.NPROC.GT.1)THEN
IF(MYID.EQ.0)THEN
DO ITZ=1,NTZ
READ(LURFAC,30) (F(L0RFAC+ITZ+NTZ*(JTZ-1)),JTZ=1,NTZ)
ENDDO
CLOSE(LURFAC)
ENDIF
CALL FArrayCast(L0RFAC+1,NTZ*NTZ,0)
ELSE
DO ITZ=1,NTZ
READ(LURFAC,30) (F(L0RFAC+ITZ+NTZ*(JTZ-1)),JTZ=1,NTZ)
ENDDO
CLOSE(LURFAC)
ENDIF
30 FORMAT (5(1PE13.6))
ELSEIF (RESTRT) THEN
IF(.NOT.NULLPR) THEN
IF(MIMD.AND.NPROC.GT.1)THEN
IF(MYID.EQ.0)THEN
CALL WRITBL
CALL WRIT40('View-factors are taken from PHI(DA) file')
CALL WRITBL
ENDIF
ELSE
CALL WRITBL
CALL WRIT40('View-factors are taken from PHI(DA) file')
CALL WRITBL
ENDIF
ENDIF
ELSE
IF(.NOT.NULLPR) THEN
CALL WRITBL
CALL WRITST
CALL WRITBL
CALL WRYT40('viewfactor calculation begins ')
CALL WRIT40('viewfactor calculation begins ')
CALL WRITBL
CALL WRITST
ENDIF
C.... Calculate viewfactor
CALL GX2DVF
IF(.NOT.NULLPR) THEN
CALL WRITST
CALL WRITBL
CALL WRYT40('viewfactor calculation ends ')
CALL WRIT40('viewfactor calculation ends ')
CALL WRITBL
CALL WRITST
ENDIF
ENDIF
IF(F(L0RFAC+1).LT.0.0.AND..NOT.NULLPR)
1 CALL WRYT40('GEBHARD exchange matrix formulation used')
C.... INDVAR
JTEM1 = LBNAME('TEM1')
JPRPS = LBNAME('PRPS')
C.... Zero location for LXDX, LYDY, LZDZ
L0AE = L0F(LXYAE)
L0AN = L0F(LXYAN)
L0AS = L0AN-1 ! ok for all except iy=1. See use below
L0AH = L0F(LXYAH)
ELSEIF(IGR.EQ.9.AND.ISC.EQ.7) THEN
C-----------------------------------------------------------------------
C.... Group 9 section 7: This part of the coding set the conductivity |
C into 3-d user store. |
C-----------------------------------------------------------------------
C.... Check that it is the first visit
IF(IGRP97.EQ.0) THEN
IGRP97=1
L0KLAM=L0F(LAMPR)
C.... Store conductivity into user f-array if first visit
ICZ=(IZ-1)*NXNY
DO 97 J=1,NY*NX
ICC=J + ICZ
F(L0COND+ICC)=-F(L0KLAM+J)
C.... Commented out initialisation of HTC store because this store
C needs to be initialised as a patch wise store.
c IF(LTURB) F(L0HTC+ICC)=F(L0COND+ICC)
97 CONTINUE
ENDIF
C-----------------------------------------------------------------------
C.... Group 9 section 14: This part of the coding set the specific |
C heat into 3-d user store. |
C-----------------------------------------------------------------------
ELSEIF(IGR.EQ.9.AND.ISC.EQ.14) THEN
IF(LTURB.AND.(IGP914.EQ.0.OR.IZ.EQ.2)) THEN
IGP914=1
L0CP=L0F(-14)
C.... Store conductivity in user F-array if first visit
ICZ=(IZ-1)*NXNY
DO 914 J=1,NY*NX
ICC=J + ICZ
F(L0CP1+ICC)=F(L0CP+J)
914 CONTINUE
ENDIF
ELSEIF(IGR.EQ.13) THEN
C-----------------------------------------------------------------------
C.... GROUP 13 |
C-----------------------------------------------------------------------
IF(INDVAR.EQ.JTEM1) THEN
C.... Store computed stanton*urel if turbulent flow and wallty patch
IF(MIMD.AND.NPROC.GT.1)THEN
LDUM = EXISTSPATCH(NPATCH)
ELSE
LDUM = .TRUE.
ENDIF
IF(LTURB.AND.LDUM) THEN
IPNO=IPNAME(NPATCH)
L0PTCH=I0PAT+10*(IPNO-1)
ITYPE=F(L0PTCH+2)
IF(ITYPE.GT.16.AND.ITYPE.LT.23) THEN
IF(ISC.GT.1.AND.ISC.LT.5) THEN
C-----------------------------------------------------------------------
C.... SECTION 2, or 3, or 4: Retrieve wall coefficient as computed by |
C FNCOEFF ST*UREL |
C-----------------------------------------------------------------------
L0CO=L0F(CO)
C.... Obtain zero location of patchwise store
ITC=-1
L0HTCC=L0HTC
DO 135 ITZ=1,NTZ
IF(MIMD.AND.NPROC.GT.1)THEN
IF(NOTEXISTZONE(ITZ))GOTO 135
ENDIF
ITC=ITC+2
IF(INT(F(L0WPNO+ITC)).EQ.IPNO) THEN
IPN=F(L0RPNO+ITZ)+0.1
L0PTCH=I0PAT+10*(IPN-1)
ITYPE=F(L0PTCH+2)+0.1
JXF=F(L0PTCH+3)+0.1
JXL=F(L0PTCH+4)+0.1
JYF=F(L0PTCH+5)+0.1
JYL=F(L0PTCH+6)+0.1
JZF=F(L0PTCH+7)+0.1
JZL=F(L0PTCH+8)+0.1
C.... Check if patch is defined in solid
L0PRPS=ABS(ANYZ(JPRPS,JZF))
ICELCK=JYF+(JXF-1)*NY
IF(F(L0PRPS+ICELCK).GT.99.) THEN
C.... If radiative zone defined in solid set required wall type
C and limits
IF(ITYPE.EQ.2) THEN
JXF=JXF+1
JXL=JXL+1
ELSEIF(ITYPE.EQ.3) THEN
JXF=JXF-1
JXL=JXL-1
ELSEIF(ITYPE.EQ.4) THEN
JYF=JYF+1
JYL=JYL+1
ELSEIF(ITYPE.EQ.5) THEN
JYF=JYF-1
JYL=JYL-1
ELSEIF(ITYPE.EQ.6) THEN
JZF=JZF+1
JZL=JZL+1
ELSEIF(ITYPE.EQ.7) THEN
JZF=JZF-1
JZL=JZL-1
ENDIF
ENDIF
IF(IZ.GE.JZF.AND.IZ.LE.JZL) THEN
CALL GETCOV(NAMPAT(IPNO),JTEM1,GCOF,GVAL)
L0HTCC=F(L0WPNO+ITC+1)+0.1
C.... If the patch number does not match any of those stored in the
C user array, RETURN because the wall patch does not coincide with
C a radiative thermal zone
C
C.... multiply arrays by first-phase density
L0RHO=L0F(DEN1)
ICZ=(IZ-1)*NXNY
ICZPV=(IZ-JZF)
DO 133 JX=JXF,JXL
ICX=(JX-1)*NY
DO 133 JY=JYF,JYL
IC=ICX+JY
ICC=IC+ICZ
ICZPV=ICZPV+1
GRHO=F(L0RHO+IC)
F(L0HTCC+ICZPV)=F(L0CP1+ICC)*GRHO*F(L0CO+IC)
C.... Zero the wall coeff when wall b.c at edge of domain and is not
C a fixed-temperature condition, because appropriate source will
C then be set by the radiative patch.
IF(EQZ(GVAL)) F(L0CO+IC)=0.0
133 CONTINUE
ENDIF
ENDIF
135 CONTINUE
ENDIF
ENDIF
ENDIF
IF(NPATCH(1:3).EQ.'@RI'.AND.LDUM) THEN
C-----------------------------------------------------------------------
C.... Group 13 section 2 and 12: This part of the coding sets the CO |
C and VAL. |
C-----------------------------------------------------------------------
C.... NB
C The method of retrieving the CC and TS is not a good one because
C it assumes that the Nth zone is numbered N. This may not always
C be the case, for example having defined the THERMAL ZONES the user
C may decide to skip the first two in a subsequent run. Therefore a
C means of obtaining the correct F-array index incremeant needs to be
C formulated. Eg loop over all zones until the PATCH number stored
C in store F(L0RPNO+......) matches the PATCH number IPNUM defined
C above, but this would be costly speed wise.
C
C.... Zone number
READ(NPATCH(1:6),'(3X,I3)',ERR=134) IZONE
IF(ISC.EQ.2) THEN
C-----------------------------------------------------------------------
C.... SECTION 2: Set coefficient CC in CO |
C-----------------------------------------------------------------------
C.... Retrieve conductive coefficient
GCC=F(L0CC +IZONE)
GCF=F(L0CCF+IZONE)
c write(14,*) 'l0cc,l0ccf,izone,gcc,gcf',
c 1 l0cc,l0ccf,izone,gcc,gcf
C.... Retrieve surface temperature TS
GTSURF=F(L0AVTS+IZONE)
C.... Zero location of SU, AP, TEM1, PRPS
L0SU=L0F(LSU)
L0AP=L0F(LAP)
L0TEM1=L0F(JTEM1)
L0TEMZ=L0TEM1
L0PRPS=L0F(JPRPS)
IF(.NOT.CARTES) L0YV=L0F(YV2D)
C.... Obtain Patch type
IPN=F(L0RPNO+IZONE)+0.1
L0PTCH=I0PAT+10*(IPN-1)
ITYPE=F(L0PTCH+2)
c.... Check that Patch is within a solid
ICELCK=IYF+(IXF-1)*NY
IF(F(L0PRPS+ICELCK).GT.99.0) THEN
C.... If patch in solid then solid/fluid interface, hence set heat
C source for fluid as well. ICINC is the increment or decrement
C needed to obtain the correct F-array index for the fluid cell.
C
IF(ITYPE.EQ.2.AND.IXF.NE.NX) THEN
ICINC= NY
ELSEIF(ITYPE.EQ.3.AND.IXF.NE.1) THEN
ICINC=-NY
ELSEIF(ITYPE.EQ.4.AND.IYF.NE.NY) THEN
ICINC= 1
ELSEIF(ITYPE.EQ.5.AND.IYF.NE.1) THEN
ICINC=-1
ELSEIF(ITYPE.EQ.6.AND.IZ.NE.NZ) THEN
ICINC= NXNY
L0TEMZ=L0F(ANYZ(JTEM1,IZ+1))
ELSEIF(ITYPE.EQ.7.AND.IZ.NE.1) THEN
ICINC=-NXNY
L0TEMZ=L0F(ANYZ(JTEM1,IZ-1))
ENDIF
DO 130 JX=IXF,IXL
ICX=(JX-1)*NY
DO 130 JY=IYF,IYL
IC=JY+ICX
C.... Get the unblocked areas to match the Patch type
IF(BFC) THEN
GAREA=F(L0B(ITYPE+3)+IC+NY*NX*(IZ-1))
ELSEIF(ITYPE.EQ.4) THEN
GAREA=F(L0AN+IC)
ELSEIF(ITYPE.EQ.5) THEN
IF(JY.NE.1) THEN ! for all except south domain boundary
GAREA=F(L0AS+IC)
ELSE
GAREA=F(L0AS+IC+1) ! at south domain boundary, use AN
ENDIF
ELSEIF(ITYPE.EQ.2 .OR. ITYPE.EQ.3) THEN
GAREA=F(L0AE+IC)
ELSEIF(ITYPE.EQ.6 .OR. ITYPE.EQ.7) THEN
GAREA=F(L0AH+IC)
ENDIF
C.... Obtain correct cell index for the fluid cell
IC=IC+ICINC
C.... ICZ is used for slab-wise stores ie for Temperature
ICZ=IC
C.... When PATCH is of type LOW or HIGH the zero location for temperature
C is obtained for the low or high slab and so ICZ does not need to be
C altered.
IF(ITYPE.EQ.6.OR.ITYPE.EQ.7) ICZ=ICZ-ICINC
F(L0SU+IC)=F(L0SU+IC)+GCF*GAREA*(GTSURF-F(L0TEMZ+ICZ))
F(L0AP+IC)=F(L0AP+IC)+GCF*GAREA
130 CONTINUE
ENDIF
DO 131 JX=IXF,IXL
ICX=(JX-1)*NY
DO 131 JY=IYF,IYL
IC=JY+ICX
C.... Get the unblocked areas to match the Patch type
IF(BFC) THEN
GAREA=F(L0B(ITYPE+3)+IC+NY*NX*(IZ-1))
ELSEIF(ITYPE.EQ.4) THEN
GAREA=F(L0AN+IC)
ELSEIF(ITYPE.EQ.5) THEN
IF(IC.NE.1) THEN ! for all except south domain boundary
GAREA=F(L0AS+IC)
ELSE
GAREA=F(L0AS+IC+1) ! at south domain boundary, use AN
ENDIF
ELSEIF(ITYPE.EQ.2 .OR. ITYPE.EQ.3) THEN
GAREA=F(L0AE+IC)
ELSEIF(ITYPE.EQ.6 .OR. ITYPE.EQ.7) THEN
GAREA=F(L0AH+IC)
ENDIF
GCAREA=GCC*GAREA
F(L0SU+IC)=F(L0SU+IC)+GCAREA*(GTSURF-F(L0TEM1+IC))
F(L0AP+IC)=F(L0AP+IC)+GCAREA
131 CONTINUE
CALL FN1(CO,0.0)
ELSEIF(ISC.EQ.13) THEN
C-----------------------------------------------------------------------
C.... Set value |
C-----------------------------------------------------------------------
CALL FN0(VAL,JTEM1)
ENDIF
ELSEIF(.NOT.LDUM)THEN
IF(ISC.EQ.2)THEN
CALL FN1(CO,0.0)
ELSEIF(ISC.EQ.13)THEN
CALL FN0(VAL,JTEM1)
ENDIF
ENDIF
ENDIF
C-----------------------------------------------------------------------
C.... GROUP 19 |
c-----------------------------------------------------------------------
ELSEIF(IGR.EQ.19.AND.ISC.EQ.2) THEN ! Start of sweep
C-----------------------------------------------------------------------
C.... SECTION 2: Initialise user F-array for TS, RF, SH, RC, and CC on |
C first sweep and update TS |
C-----------------------------------------------------------------------
C.... Find wall patches if any for the thermal patches
IF(ISWEEP.EQ.FSWEEP) THEN
C.... Check that the PATCH type of radiative zones is area
DO 1921 ITZ=1,NTZ
C.... Obtain Patch type and limits
IPN=F(L0RPNO+ITZ)+0.1
L0PTCH=I0PAT+10*(IPN-1)
ITYPE=F(L0PTCH+2)
C.... Error message if wrong patch type
IF(ITYPE.LT.2.OR.ITYPE.GT.7) GO TO 1924
1921 CONTINUE
C.... For turbulent flows obtain patch number of wall type patches for
C the radiative patches.
IF(LTURB) THEN
C.... Initialise patch limits of pervious wall patch used in allocating
C patch wise index
CALL SUB3(JPXF,0,JPXL,0,JPYF,0)
CALL SUB3(JPYL,0,JPZF,0,JPZL,0)
C.... Initialise index pointer
L0HTCC=L0HTC
ICTZ=-1
DO 40 ITZ=1,NTZ
IF(MIMD.AND.NPROC.GT.1)THEN
IF(NOTEXISTZONE(ITZ))GOTO 40
ENDIF
C.... Obtain Patch type and limits
ICTZ=ICTZ+1
IPN=F(L0RPNO+ITZ)+0.1
L0PTCH=I0PAT+10*(IPN-1)
ITYPE=F(L0PTCH+2)
JXF=F(L0PTCH+3)+0.1
JXL=F(L0PTCH+4)+0.1
JYF=F(L0PTCH+5)+0.1
JYL=F(L0PTCH+6)+0.1
JZF=F(L0PTCH+7)+0.1
JZL=F(L0PTCH+8)+0.1
C.... Check if patch is defined in solid
L0PRPS=ABS(ANYZ(JPRPS,JZF))
ICELCK=JYF+(JXF-1)*NY
IF(F(L0PRPS+ICELCK).GT.99.0) THEN
C.... If radiative zone defined in solid set required wall type
C and limits
IF(ITYPE.EQ.2) THEN
CALL SUB3(JWXF,JXF+1,JWXL,JXL+1,JWYF,JYF)
CALL SUB3(JWYL,JYL,JWZF,JZF,JWZL,JZL)
IWLTYP=18
ELSEIF(ITYPE.EQ.3) THEN
CALL SUB3(JWXF,JXF-1,JWXL,JXL-1,JWYF,JYF)
CALL SUB3(JWYL,JYL,JWZF,JZF,JWZL,JZL)
IWLTYP=17
ELSEIF(ITYPE.EQ.4) THEN
CALL SUB3(JWXF,JXF,JWXL,JXL,JWYF,JYF+1)
CALL SUB3(JWYL,JYL+1,JWZF,JZF,JWZL,JZL)
IWLTYP=20
ELSEIF(ITYPE.EQ.5) THEN
CALL SUB3(JWXF,JXF,JWXL,JXL,JWYF,JYF-1)
CALL SUB3(JWYL,JYL-1,JWZF,JZF,JWZL,JZL)
IWLTYP=19
ELSEIF(ITYPE.EQ.6) THEN
CALL SUB3(JWXF,JXF,JWXL,JXL,JWYF,JYF)
CALL SUB3(JWYL,JYL,JWZF,JZF+1,JWZL,JZL+1)
IWLTYP=22
ELSEIF(ITYPE.EQ.7) THEN
CALL SUB3(JWXF,JXF,JWXL,JXL,JWYF,JYF)
CALL SUB3(JWYL,JYL,JWZF,JZF-1,JWZL,JZL-1)
IWLTYP=21
ENDIF
ELSE
C.... Otherwise limits for the wall type are the
C same as for the radiative patch.
CALL SUB3(JWXF,JXF,JWXL,JXL,JWYF,JYF)
CALL SUB3(JWYL,JYL,JWZF,JZF,JWZL,JZL)
IWLTYP=ITYPE+15
ENDIF
F(L0WPNO+ITZ+ICTZ)=0.0
F(L0WPNO+ITZ+ICTZ+1)=L0HTC
C.... loop over all patches to find required wall patch
LFOUND=.FALSE.
DO 50 IPN0=1,NUMPAT
L0PTCH=I0PAT+10*(IPN0-1)
ITYPE=F(L0PTCH+2)
IF(ITYPE.EQ.IWLTYP) THEN
JXF=F(L0PTCH+3)+0.1
JXL=F(L0PTCH+4)+0.1
JYF=F(L0PTCH+5)+0.1
JYL=F(L0PTCH+6)+0.1
JZF=F(L0PTCH+7)+0.1
JZL=F(L0PTCH+8)+0.1
C.... The following checks allows a wall patch to coincide with several
C thermal zones, but it does not allow for a single thermal zone to
C coinicide with several wall patches.
IF(ITYPE.EQ.17.OR.ITYPE.EQ.18) THEN
C.... EAST OR WEST type so jxf=jxl=jwxf=jwxl
IF(JXF.EQ.JWXF.AND.JXL.EQ.JWXL.AND.JYF.LE.JWYF.
1 AND.JYL.GE.JWYL.AND.JZF.LE.JWZF.AND.JZL.GE.JWZL)
1 LFOUND=.TRUE.
ELSEIF(ITYPE.EQ.19.OR.ITYPE.EQ.20) THEN
C.... NORTH OR SOUTH type so jyf=jyl=jwyf=jwyl
IF(JYF.EQ.JWYF.AND.JYL.EQ.JWYL.AND.JXF.LE.JWXF.
1 AND.JXL.GE.JWXL.AND.JZF.LE.JWZF.AND.JZL.GE.JWZL)
1 LFOUND=.TRUE.
ELSEIF(ITYPE.EQ.21.OR.ITYPE.EQ.22) THEN
C.... HIGH OR LOW type so jzf=jzl=jwzf=jwzl
IF(JZF.EQ.JWZF.AND.JZL.EQ.JWZL.AND.JYF.LE.JWYF.
1 AND.JYL.GE.JWYL.AND.JXF.LE.JWXF.AND.JXL.GE.JWXL)
1 LFOUND=.TRUE.
ENDIF
ENDIF
IF(LFOUND) THEN
LFOUND=.FALSE.
F(L0WPNO+ITZ+ICTZ)=FLOAT(IPN0)
GO TO 60
ENDIF
50 CONTINUE
GO TO 40
C.... Set index
60 L0HTCC=L0HTCC+(JPXL-JPXF+1)*(JPYL-JPYF+1)*(JPZL-JPZF+1)
F(L0WPNO+ITZ+ICTZ+1)=FLOAT(L0HTCC)
CALL SUB3(JPXF,JWXF,JPXL,JWXL,JPYF,JWYF)
CALL SUB3(JPYL,JWYL,JPZF,JWZF,JPZL,JWZL)
40 CONTINUE
ENDIF
C.... Initialise user F-array
IF(LINIL) THEN
LINIL=.FALSE.
IF(QEQ(FIINIT(JTEM1),READFI)) THEN
C.... When restart then for fixed flux radiative zones need to reset
C the VALue to GRND1
DO 1910 ITZ=1,NTZ
IPN=F(L0RPNO+ITZ)+0.1
CALL GETCOV(NAMPAT(IPN),JTEM1,GCOF,GVAL)
IF(.NOT.GRNDNO(1,GVAL).AND.GRNDNO(1,GCOF)) THEN
F(L0QFLX+ITZ)=GVAL
F(I0PHI+3)=GRND1
ELSE
F(L0QFLX+ITZ)=0.0
ENDIF
C.... Store k/d for thin plates and node number of adjacent node
IF(GTZ(GVAL).AND.GTZ(GCOF)) THEN
F(L0KDM+ITZ)=GCOF
F(L0NTS+ITZ)=GVAL
F(I0PHI+3)=GRND1
F(I0PHI+2)=GRND1
ELSE
F(L0KDM+ITZ)=0.0
F(L0NTS+ITZ)=0.0
ENDIF
1910 CONTINUE
ELSE
ICNT=-1
DO 191 ITZ=1,NTZ
ICNT=ICNT+2
IPN=F(L0RPNO+ITZ)+0.1
CALL GETCOV(NAMPAT(IPN),JTEM1,GCOF,GVAL)
IF(GRNDNO(1,GVAL).OR.GRNDNO(1,GCOF)) THEN
F(L0AVTS+ITZ)=FIINIT(JTEM1)
ELSE
C.... assign fixed wall temp 140792
F(L0AVTS+ITZ)=GVAL
IF(LTZ(GVAL)) F(L0AVTS+ITZ)=FIINIT(JTEM1)
ENDIF
C.... Store heat flux value if applicable and reset VALue to GRND1
IF(.NOT.GRNDNO(1,GVAL).AND.GRNDNO(1,GCOF)) THEN
F(L0QFLX+ITZ)=GVAL
F(I0PHI+3)=GRND1
ELSE
F(L0QFLX+ITZ)=0.0
ENDIF
C.... Store k/d for thin plates and node number of adjacent node
IF(GTZ(GVAL).AND.GTZ(GCOF)) THEN
F(L0KDM+ITZ)=GCOF
F(L0NTS+ITZ)=GVAL
F(I0PHI+3)=GRND1
F(I0PHI+2)=GRND1
ELSE
F(L0KDM+ITZ)=0.0
F(L0NTS+ITZ)=0.0
ENDIF
CALL SUB3R(F(L0AVRF+ITZ),F(L0QFLX+ITZ),F(L0AVSH+ITZ),
1 0.0,F(L0RC+ITZ),0.0)
CALL SUB3R(F(L0CC+ITZ),0.0,F(L0CCF+ITZ),0.0,
1 F(L0QEXT+ITZ),0.0)
CALL SUB3R(F(L0XCOF+ITZ),0.0,F(L0EXRD+ICNT),0.0,
1 F(L0EXRD+ICNT+1),0.0)
F(L0EXLN+ICNT)=0.0
F(L0EXLN+ICNT+1)=0.0
191 CONTINUE
ENDIF
C.... If there are external heat transfer PATCHES then obtain and store
C the coefficient and external temperature.
DO 1922 IEXPN=1,NUMPAT
IF(NAMPAT(IEXPN)(1:3).EQ.'@ER' .OR.
1 NAMPAT(IEXPN)(1:3).EQ.'@EL') THEN
IF(MIMD.AND.NPROC.GT.1) THEN
IF(.NOT.EXISTSPATCH(NAMPAT(IEXPN))) GOTO 1922
ENDIF
L0PTCH=I0PAT+10*(IEXPN-1)
ITYPE=F(L0PTCH+2)
C.... Error message if wrong patch type
IF(ITYPE.LT.2.OR.ITYPE.GT.7) GO TO 1925
JXF=F(L0PTCH+3)+0.1
JXL=F(L0PTCH+4)+0.1
JYF=F(L0PTCH+5)+0.1
JYL=F(L0PTCH+6)+0.1
JZF=F(L0PTCH+7)+0.1
JZL=F(L0PTCH+8)+0.1
ICNT=-1
DO 1923 ITZ=1,NTZ
IF(MIMD.AND.NPROC.GT.1)THEN
IF(NOTEXISTZONE(ITZ))GOTO 1923
ENDIF
ICNT=ICNT+2
C.... Obtain Patch type and limits
IPN=F(L0RPNO+ITZ)+0.1
L0PTCH=I0PAT+10*(IPN-1)
IRTYPE=F(L0PTCH+2)
F(L0EXRD+ICNT )=0.0
F(L0EXRD+ICNT+1)=0.0
F(L0EXLN+ICNT )=0.0
F(L0EXLN+ICNT+1)=0.0
IF(IRTYPE.EQ.ITYPE) THEN
JRXF=F(L0PTCH+3)+0.1
JRXL=F(L0PTCH+4)+0.1
JRYF=F(L0PTCH+5)+0.1
JRYL=F(L0PTCH+6)+0.1
JRZF=F(L0PTCH+7)+0.1
JRZL=F(L0PTCH+8)+0.1
IF(JXF.EQ.JRXF.AND.JXL.EQ.JRXL.AND.JYF.EQ.JRYF.AND.
1 JYL.EQ.JRYL.AND.JZF.EQ.JRZF.AND.JZL.EQ.JRZL) THEN
C.... The internal radiative patch is found so store the coefficient
C and the external temperature.
CALL GETCOV(NAMPAT(IEXPN),JTEM1,GCOF,GVAL)
IF(NAMPAT(IEXPN)(1:3).EQ.'@ER') THEN
F(L0EXRD+ICNT )=GCOF
F(L0EXRD+ICNT+1)=GVAL
C.... Zero the COefficient to avoid duplication of source by EARTH
F(I0PHI+2)=0.0
ELSEIF(NAMPAT(IEXPN)(1:3).EQ.'@EL') THEN
F(L0EXLN+ICNT )=GCOF
F(L0EXLN+ICNT+1)=GVAL
C.... Zero the COefficient to avoid duplication of source by EARTH
F(I0PHI+2)=0.0
ENDIF
ENDIF
ENDIF
1923 CONTINUE
ENDIF
1922 CONTINUE
ENDIF
ENDIF
C.... Update surface temperatures
IF(ISWEEP.NE.FSWEEP.AND.LOOPZ.EQ.1) THEN
DO 192 ITZ=1,NTZ
GQMET=0.0
IF(GTZ(F(L0NTS+ITZ))) THEN
INTZ=F(L0NTS+ITZ)+0.1
GQMET=F(L0KDM+ITZ)*(F(L0AVTS+ITZ)-F(L0AVTS+INTZ))
ENDIF
F(L0AVTS+ITZ)=F(L0AVTS+ITZ)+
1 (F(L0AVRF+ITZ)-F(L0AVSH+ITZ)+
1 F(L0QFLX+ITZ)+F(L0QEXT+ITZ)-GQMET)/
1 (TINY+F(L0CC +ITZ)+F(L0KDM +ITZ)-
1 F(L0RC +ITZ)+F(L0CCF +ITZ)-F(L0XCOF+ITZ))
192 CONTINUE
ENDIF
C
C.... Update aperture zone temperature
C
DO 198 ITZ=1,NTZ
IPN=F(L0RPNO+ITZ)+0.1
CALL GETCOV(NAMPAT(IPN),JTEM1,GCOF,GVAL)
IF(GRNDNO(1,GVAL) .AND. EQZ(GCOF)) THEN
IF(MIMD.AND.NPROC.GT.1)THEN
IF(NOTEXISTZONE(ITZ))GOTO 200
ENDIF
C.... obtain patch type and limits
L0PTCH=I0PAT+10*(IPN-1)
ITYPE=F(L0PTCH+2)
JXF=F(L0PTCH+3)+0.1
JXL=F(L0PTCH+4)+0.1
JYF=F(L0PTCH+5)+0.1
JYL=F(L0PTCH+6)+0.1
JZF=F(L0PTCH+7)+0.1
JZL=F(L0PTCH+8)+0.1
CALL SUB2R(GSUMBT,0.0,GSUMVL,0.0)
DO 199 JZ=JZF,JZL
L0TEMP=L0F(ANYZ(JTEM1,JZ))
L0VOL=L0F(ANYZ(VOL,JZ))
IF(BFC) L0VOL=L0B(VOLUME)
ICZ = (JZ-1)*NXNY
DO 199 JX=JXF,JXL
ICX=(JX-1)*NY
DO 199 JY=JYF,JYL
IC=JY + ICX
IF(BFC) THEN
GSUMVL=GSUMVL+F(L0VOL+IC+ICZ)
GSUMBT=GSUMBT+F(L0TEMP+IC)*F(L0VOL+IC+ICZ)
ELSE
GSUMVL=GSUMVL+F(L0VOL+IC)
GSUMBT=GSUMBT+F(L0TEMP+IC)*F(L0VOL+IC)
ENDIF
199 CONTINUE
200 IF(MIMD.AND.NPROC.GT.1)CALL GLSUM2(GSUMBT,GSUMVL)
F(L0AVTS+ITZ)=GSUMBT/(GSUMVL+TINY)
ENDIF
198 CONTINUE
ELSEIF(IGR.EQ.19.AND.ISC.EQ.3) THEN
C.... Group 9 section 7 visit counter
IGRP97=0
C.... Group 9 section 14 visit counter
IGP914=0
ELSEIF(IGR.EQ.19.AND.ISC.EQ.7) THEN
C-----------------------------------------------------------------------
C.... SECTION 7: Calculate RF, RC, SH and CC, the last can be a |
C function of temperature if conductivity is a f(tem1) |
C-----------------------------------------------------------------------
ITC=-1
DO 193 ITZ=1,NTZ
ITC=ITC+2
CALL SUB3R(GSUMSH,0.0,GSUMAR,0.0,GSUMCC,0.0)
CALL SUB3R(GSUMRC,0.0,GSUMRF,0.0,GTSI,F(L0AVTS+ITZ))
c write(14,*)'gtsi=',gtsi
c write(14,*)'l0avts, itz ',l0avts, itz
CALL SUB3R(GSUMSS,0.0,GSUMSF,0.0,GSUMCF,0.0)
c call fcheck(100,'gxs2sr')
C
C.... Compute RF and RC. Modified formalation for net radiative flux,
C based on specification of GEBHARD exchange matrix, in which
C diagonal elements will be -ve.
C
IF(F(L0RFAC+1).LT.0.0) THEN
DO 197 JTZ=1,NTZ
GTSJ=F(L0AVTS+JTZ)
GRDFIJ=F(L0RFAC+ITZ+NTZ*(JTZ-1))
IF(JTZ.NE.ITZ) GSUMRC=GSUMRC+GRDFIJ
GSUMRF=GSUMRF+GRDFIJ*(GTSJ+TMP1A)**4
197 CONTINUE
call fcheck(101,'gxs2sr')
ELSE
DO 194 JTZ=1,NTZ
GTSJ=F(L0AVTS+JTZ)
GRDFIJ=F(L0RFAC+ITZ+NTZ*(JTZ-1))
GSUMRC=GSUMRC+GRDFIJ
GSUMRF=GSUMRF+GRDFIJ*((GTSJ+TMP1A)**4-(GTSI+TMP1A)**4)
194 CONTINUE
ENDIF
call fcheck(102,'gxs2sr')
F(L0RC+ITZ)=-4.*(GTSI+TMP1A)**3*GSUMRC
F(L0AVRF+ITZ)=GSUMRF
IPN=F(L0RPNO+ITZ)+0.1
CALL GETCOV(NAMPAT(IPN),JTEM1,GCOF,GVAL)
IF(GRNDNO(1,GVAL)) THEN
IF(MIMD.AND.NPROC.GT.1)THEN
IF(NOTEXISTZONE(ITZ))GOTO 201
ENDIF
C.... Obtain patch type and limits
L0PTCH=I0PAT+10*(IPN-1)
ITYPE=F(L0PTCH+2)+0.1
JXF = F(L0PTCH+3)+0.1
JXL = F(L0PTCH+4)+0.1
JYF = F(L0PTCH+5)+0.1
JYL = F(L0PTCH+6)+0.1
JZF = F(L0PTCH+7)+0.1
JZL = F(L0PTCH+8)+0.1
C
C.... ICINC is either an increment or a decrement required
C to locate the fluid cells
ICINC=0
IF(ITYPE.EQ.2.OR.ITYPE.EQ.3) THEN
IF(BFC) THEN
L0DIS=L0B(DHXPE)
ELSE
L0DIS=L0F(DXU2D)
ENDIF
IF(ITYPE.EQ.2.AND.JXF.NE.NX) ICINC= NY
IF(ITYPE.EQ.3.AND.JXF.NE. 1) ICINC=-NY
ELSEIF(ITYPE.EQ.6.OR.ITYPE.EQ.7) THEN
IF(BFC) THEN
L0DIS=L0B(DHZPH)
ELSE
L0DIS=L0F(DZWNZ)
ENDIF
ELSEIF(ITYPE.EQ.4.OR.ITYPE.EQ.5) THEN
IF(BFC) THEN
L0DIS=L0B(DHYPN)
ELSE
L0DIS=L0F(DYV2D)
c write(14,*)'dyv2d,l0dis',dyv2d,l0dis
ENDIF
IF(ITYPE.EQ.4.AND.JYF.NE.NY) ICINC= 1
IF(ITYPE.EQ.5.AND.JYF.NE. 1) ICINC=-1
ENDIF
c call fcheck(102,'gxs2sr')
IF(.NOT.CARTES) L0YV=L0F(YV2D)
C.... Calculate SH and CC
DO 195 JZ=JZF,JZL
C.... Zero locations
L0TEMP=L0F(ANYZ(JTEM1,JZ))
L0PRPS=L0F(ANYZ(JPRPS,JZ))
L0TEMZ=L0TEMP
IF(ITYPE.EQ.6.AND.JZF.NE.NZ) THEN
ICINC=NXNY
JZZ=JZ+1
L0TEMZ=L0F(ANYZ(JTEM1,JZZ))
ELSEIF(ITYPE.EQ.7.AND.JZF.NE.1) THEN
ICINC=-NXNY
JZZ=JZ-1
L0TEMZ=L0F(ANYZ(JTEM1,JZZ))
ENDIF
ICZ = (JZ-1)*NXNY
C.... L0HTCC and IHC are index and increment for heat transfer coefficient
C which is stored as a patch wise variable in the users f-array.
IF(LTURB) THEN
L0HTCC=F(L0WPNO+ITC+1)+0.1
IHC=(JZ-JZF)
ENDIF
DO 195 JX=JXF,JXL
ICX = (JX-1)*NY
DO 195 JY=JYF,JYL
IC =JY + ICX
IC1 =IC + ICINC
IF(BFC) THEN
GAREA=F(L0B(ITYPE+3)+IC+NY*NX*(IZ-1))
c GAREA=F(L0B(ITYPE+3)+IC+NY*NZ*(IZ-1))
ELSEIF(ITYPE.EQ.4) THEN
GAREA=F(L0AN+IC)
ELSEIF(ITYPE.EQ.5) THEN
GAREA=F(L0AS+IC)
IF(JY.NE.1) THEN ! for all except south domain boundary
GAREA=F(L0AS+IC)
ELSE
GAREA=F(L0AS+IC+1) ! at south domain boundary, use AN
ENDIF
ELSEIF(ITYPE.EQ.2 .OR. ITYPE.EQ.3) THEN
GAREA=F(L0AE+IC)
ELSEIF(ITYPE.EQ.6 .OR. ITYPE.EQ.7) THEN
GAREA=F(L0AH+IC)
ENDIF
IF(LTURB) IHC=IHC+1
GRCD=0.5*F(L0DIS+IC)
c write(14,*)'l0dis, ic, grcd 1',l0dis, ic, grcd
IF(ITYPE.EQ.6.OR.ITYPE.EQ.7) GRCD=0.5*F(L0DIS+JZ)
c write(14,*)'grcd 2',grcd
IF(BFC) GRCD=0.5*F(L0DIS+IC+ICZ)
GRCD =1.0/(GRCD+TINY)
c write(14,*)'grcd 3',grcd
GCONZ=F(L0COND+IC+ICZ)
IF(F(L0PRPS+IC).LE.99.0) THEN
IF(LTURB.AND.GTZ(F(L0WPNO+ITC))) THEN
GCONZ=F(L0HTCC+IHC)
GRCD=1.0
ENDIF
ENDIF
GSUMSH=((GTSI-F(L0TEMP+IC))*GCONZ*GRCD*GAREA)+GSUMSH
c write(14,*)' F(L0TEMP+IC)', F(L0TEMP+IC)
c write(14,*)' GCONZ ', GCONZ
c write(14,*)' GRCD ', GRCD
c write(14,*)' GAREA ', GAREA
c write(14,*) 'ic, gsumsh',ic,gsumsh
GSUMCC=GSUMCC+GCONZ*GRCD*GAREA
C.... When at solid/fluid interface compute fluid side heat source due
C to radiation
IF(F(L0PRPS+IC).GT.99.0) THEN
GCONZ=F(L0COND+IC1+ICZ)
IF(LTURB.AND.GTZ(F(L0WPNO+ITC))) THEN
GRCD=1.0
ELSE
IF(BFC) THEN
GRCD=F(L0DIS+IC1+ICZ)
ELSEIF(ITYPE.EQ.6) THEN
GRCD=F(L0DIS+JZ+1)
ELSEIF(ITYPE.EQ.7) THEN
GRCD=F(L0DIS+JZ-1)
ELSE
GRCD =F(L0DIS+IC1)
ENDIF
GRCD=2.0/(GRCD + TINY)
ENDIF
IF(LTURB) GCONZ=F(L0HTCC+IHC)
IF(ITYPE.EQ.6.OR.ITYPE.EQ.7) IC1=IC1-ICINC
GCRARE=GCONZ*GRCD*GAREA
GSUMSF=GSUMSF + (GTSI-F(L0TEMZ+IC1)) * GCRARE
GSUMCF=GSUMCF + GCRARE
ENDIF
GSUMAR = GSUMAR + GAREA
c write(14,*)'gsumar',gsumar
195 CONTINUE
c call fcheck(104,'gxs2sr')
201 IF(MIMD.AND.NPROC.GT.1)THEN
CALL GLSUM5(GSUMSH,GSUMAR,GSUMCC,GSUMSF,GSUMCF)
ENDIF
RECDEN=1.0/(GSUMAR+TINY)
c call writ3r('recden ',recden,'gsumsh ',gsumsh,
c 1 'gsumsf ',gsumsf)
F(L0AVSH+ITZ) = GSUMSH*RECDEN + GSUMSF*RECDEN
c call writ1r('avsh ',F(L0AVSH+ITZ))
F(L0CC +ITZ) = GSUMCC*RECDEN
F(L0CCF +ITZ) = GSUMCF*RECDEN
c call fcheck(105,'gxs2sr')
C.... Compute external heat transfer
C.... Radiative
GQEXRD=F(L0EXRD+ITC)*((F(L0EXRD+ITC+1)+TMP1A)**4-
1 (GTSI+TMP1A)**4)
C.... Linear
GQEXLN=F(L0EXLN+ITC)*(F(L0EXLN+ITC+1)-GTSI)
F(L0QEXT+ITZ)=GQEXRD+GQEXLN
C.... Update external coefficient for the heat balance equation
F(L0XCOF+ITZ)=-4.0*F(L0EXRD+ITC)*(GTSI+TMP1A)**3-
1 F(L0EXLN+ITC)
ELSE
C.... When fixed temp set SH to RF
F(L0AVSH+ITZ)=F(L0AVRF+ITZ)
ENDIF
193 CONTINUE
C
C.... Print out for each zone the number, temperature, radiative flux
C and conductive flux
C
IF(ISWEEP.EQ.LSWEEP.AND..NOT.NULLPR) THEN
IF(MIMD.AND.NPROC.GT.1)THEN
IF(MYID.EQ.0)THEN
WRITE(LUPR1,*) ' zone printout '
WRITE(LUPR1,*) ' zone tsur radf'
1 //' htrf'
DO ITZ=1,NTZ
WRITE(LUPR1,30) FLOAT(ITZ),F(L0AVTS+ITZ),F(L0AVRF+ITZ),
1 F(L0AVSH+ITZ)
ENDDO
ENDIF
ELSE
WRITE(LUPR1,*) ' zone printout '
WRITE(LUPR1,*) ' zone tsur radf htrf'
DO ITZ=1,NTZ
WRITE(LUPR1,30) FLOAT(ITZ),F(L0AVTS+ITZ),F(L0AVRF+ITZ),
1 F(L0AVSH+ITZ)
ENDDO
ENDIF
ENDIF
ENDIF
NAMSUB='gxs2sr'
RETURN
C
C.... Write out error message if PATCH name is not of
C the format @RI###, where ### are digits.
C
134 CONTINUE
IF(MIMD.AND.NPROC.GT.1)THEN
IF(MYID.EQ.0)THEN
CALL WRITST
CALL WRITBL
WRITE(LUPR1,*) ' THERMAL ZONE PATCH NAME ',NPATCH
WRITE(LUPR1,*) ' IS NOT OF THE REQUIRED FORM. '
WRITE(LUPR1,*) ' IT SHOULD START WITH @RI FOLLOWED BY
1 THREE DIGITS.'
WRITE(LUPR1,*) ' EXAMPLES OF ACCEPTED FORMS ARE: '
WRITE(LUPR1,*) ' @RI001 or @RI100.'
CALL WRITST
CALL WRITBL
ENDIF
ELSE
CALL WRITST
CALL WRITBL
WRITE(LUPR1,*) ' THERMAL ZONE PATCH NAME ',NPATCH
WRITE(LUPR1,*) ' IS NOT OF THE REQUIRED FORM. '
WRITE(LUPR1,*) ' IT SHOULD START WITH @RI FOLLOWED BY
1 THREE DIGITS.'
WRITE(LUPR1,*) ' EXAMPLES OF ACCEPTED FORMS ARE: '
WRITE(LUPR1,*) ' @RI001 or @RI100.'
CALL WRITST
CALL WRITBL
ENDIF
CALL SET_ERR(299, 'Error. See result file.',1)
RETURN
1924 CONTINUE
IF(MIMD.AND.NPROC.GT.1)THEN
IF(MYID.EQ.0)THEN
CALL WRITST
CALL WRITBL
WRITE(LUPR1,*) ' Group ',IGR,'Section ',ISC
WRITE(LUPR1,*) ' Incorrect PATCH TYPE for the',
1 ' radiating surface'
WRITE(LUPR1,*) 'PHOENICS PATCH NUMBER =',IPN
WRITE(LUPR1,*) 'PATCH has to be of an area type '
ENDIF
ELSE
CALL WRITST
CALL WRITBL
WRITE(LUPR1,*) ' Group ',IGR,'Section ',ISC
WRITE(LUPR1,*) ' Incorrect PATCH TYPE for the',
1 ' radiating surface'
WRITE(LUPR1,*) 'PHOENICS PATCH NUMBER =',IPN
WRITE(LUPR1,*) 'PATCH has to be of an area type '
ENDIF
CALL SET_ERR(300, 'Error. See result file.',1)
RETURN
1925 CONTINUE
IF(MIMD.AND.NPROC.GT.1)THEN
IF(MYID.EQ.0)THEN
CALL WRITST
CALL WRITBL
WRITE(LUPR1,*) ' Group ',IGR,'Section ',ISC
WRITE(LUPR1,*) ' Incorrect PATCH TYPE for the',
1 ' external heat transfer surface'
WRITE(LUPR1,*) 'PHOENICS PATCH NUMBER =',IEXPN
WRITE(LUPR1,*) 'PATCH has to be of an area type '
ENDIF
ELSE
CALL WRITST
CALL WRITBL
WRITE(LUPR1,*) ' Group ',IGR,'Section ',ISC
WRITE(LUPR1,*) ' Incorrect PATCH TYPE for the',
1 ' external heat transfer surface'
WRITE(LUPR1,*) 'PHOENICS PATCH NUMBER =',IEXPN
WRITE(LUPR1,*) 'PATCH has to be of an area type '
ENDIF
CALL SET_ERR(301, 'Error. See result file.',1)
END
c