cGxs2sr.for c c

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