cgxbfgr.for c c

c
c
c
C     How the subroutine is activated: GX2DVF is called from
C     GXS2SR in g1 s2 if the user has NOT created a file called
C     'RADI.DAT' in the local directory containing the required
C     radiation exchange matrix.
C
C     Limitations: 1-d and 2-d cartesian/bfc coordinates only. Note
C     that the coding assumes that the radiative exchange
C     coefficients provided by the user must include emissivities and
C     the Stefan-Boltzmann constant.
C     Click here
c     for full information
C-----------------------------------------------------------------
      SUBROUTINE GX2DVF
      INCLUDE 'patcmn'
      INCLUDE 'farray'
      INCLUDE 'grdloc'
      INCLUDE 'satgrd'
      INCLUDE 'grdear'
      INCLUDE 'satear'
      INCLUDE 'parear'
      PARAMETER (NPNAM=200)
      PARAMETER (NLG=100, NIG=100, NRG=200, NCG=100)
      PARAMETER (JRTMAX=5)
C.... maximal number of viewfactor matrices =
C     maximal vnumber of spectral bands for str. media
      PARAMETER (JLMX=5)
C.... maximal number of radiating patches
      PARAMETER (JPTHMX=50)
C.... number of spectral intervals
      PARAMETER (JSPMX=60)
C.... number of spectral intervals
      PARAMETER (JBNDMX=7)
C.... maximal number of different materials in this
C     calculation (to limit the storage for optical data)
      PARAMETER (JMTLIM=10)
C.... range of material numbers in data file
      PARAMETER (JMTIND=500)
C
      COMMON/LGRND/LG(NLG)/IGRND/IG(NIG)/RGRND/RG(NRG)/CGRND/CG(NCG)
      LOGICAL LG
      CHARACTER*4 CG
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
      CHARACTER*6 NAMSUB
      LOGICAL STPRPS,STEMMS,STEPOR,STNPOR,STHPOR,QEQ,QNE,GBHMTX
C
      INTEGER JINT(JSPMX),JSTRNT(JSPMX),JMTCTR(JMTIND)
      REAL GRELP(JPTHMX,JBNDMX),GSURF(JPTHMX,JBNDMX)
      REAL OPT(JMTLIM,JSPMX,9),GE(JPTHMX,JSPMX),GR(JPTHMX,JSPMX)
C.... New array providing the spectral weight
C     for every surface must be introduced
      REAL GW(JPTHMX,JSPMX)
      LOGICAL JSPFLG,SURFL
C
      REAL GCOR1(100,100),GCOR2(100,100)
      REAL A(JPTHMX,JPTHMX),B(JPTHMX,JPTHMX),R(JPTHMX,JPTHMX)
C...  Local arrays, integers etc.
      REAL GVIEW(JPTHMX,JPTHMX,JLMX),GPATA(JPTHMX)
      REAL GLOJ1R(JRTMAX),GLOJ2R(JRTMAX)
      REAL GHIJ1R(JRTMAX),GHIJ2R(JRTMAX)
      REAL GLOI1R(JRTMAX),GLOI2R(JRTMAX)
      REAL GHII1R(JRTMAX),GHII2R(JRTMAX)
      REAL GCNTBF(JLMX,JRTMAX,JRTMAX),GTEST(JLMX),GDIR(2)
      LOGICAL BEFORE,ABSORP
      INTEGER J2DDIR(7)
      INTEGER JSHFT1(2,2),JSHFT2(2,2)
      INTEGER J2DLO1(2,2),J2DLO2(2,2),J2DHI1(2,2),J2DHI2(2,2)
      INTEGER J2DA1(2,2,2),J2DA2(2,2,2),J2DB1(2,2,2,3),J2DB2(2,2,2,3)
      INTEGER J2DSH1(2,2,2,3),J2DSH2(2,2,2,3),J2DTYP(2,3,7,3)
      DATA JSHFT1 / -1, 0, 1, 0 /
      DATA JSHFT2 /  0,-1, 0, 1 /
      DATA J2DSH1 /1, 0,-1, 0, 1, 0,-1, 0, 0,-1, 0,-1, 0, 1, 0, 1,
     1             0, 1, 0, 1, 0,-1, 0,-1 /
      DATA J2DSH2 /0, 1, 0,-1, 0, 1, 0,-1,-1, 0,-1, 0, 1, 0, 1, 0,
     1             1, 0, 1, 0,-1, 0,-1, 0 /
      DATA J2DLO1 /0,0,1,0/
      DATA J2DLO2 /0,0,0,1/
      DATA J2DHI1 /0,1,1,1/
      DATA J2DHI2 /1,0,1,1/
      DATA J2DA1 / 1,0,0,0,1,1,0,1 /
      DATA J2DA2 / 0,1,0,0,1,1,1,0 /
      DATA J2DB1 / 1,1,0,1,1,0,0,0,0,0,1,0,0,1,1,1,0,1,1,1,0,0,1,0 /
      DATA J2DB2 / 1,1,1,0,0,1,0,0,0,0,0,1,1,0,1,1,1,0,1,1,0,0,0,1 /
      DATA J2DTYP
     1 / 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,4,6,7,7,6,5,5,6,7,
     1   7,6,6,6,4,5,5,4,7,7,4,5,5,4,
     1   0,0,0,0,0,0,2,2,6,7,7,6,3,3,6,7,7,6,0,0,0,0,0,0,0,0,0,0,
     1   0,0,6,6,2,3,3,2,7,7,2,3,3,2,
     1   0,0,0,0,0,0,2,2,4,5,5,4,3,3,4,5,5,4,4,4,2,3,3,2,5,5,2,3,
     1   3,2,0,0,0,0,0,0,0,0,0,0,0,0 /
C
C.... number of testrays on every test surface, defines accuracy
C     of search for obstructions to the second approxmation
      JRT = 2
C
      NAMSUB = 'VIEW2D'
      OPEN(3,FILE='rad.dat',STATUS='UNKNOWN')
      IZ = 1
      GVIEW=0.0; R=0.0
      STPRPS =.FALSE.
      STEMMS =.FALSE.
      STEPOR =.FALSE.
      STNPOR =.FALSE.
      STHPOR =.FALSE.
      JPRPS  = LBNAME('PRPS')
      JEMMS  = LBNAME('EMMS')
      JEPOR  = LBNAME('EPOR')
      JNPOR  = LBNAME('NPOR')
      JHPOR  = LBNAME('HPOR')
      IF(STORE(JPRPS)) STPRPS =.TRUE.
      IF(STORE(JEMMS)) STEMMS =.TRUE.
      IF(STORE(JEPOR)) STEPOR =.TRUE.
      IF(STORE(JNPOR)) STNPOR =.TRUE.
      IF(STORE(JHPOR)) STHPOR =.TRUE.
      IF(.NOT.NULLPR) THEN
        WRITE(LUPR1,*) 'Print-out from sub-routine VIEW2D'
        CALL WRITST
        WRITE(LUPR1,*) 'Store PRPS=',STPRPS
        WRITE(LUPR1,*) 'Store EMMS=',STEMMS
        WRITE(LUPR1,*) 'Store EPOR=',STEPOR
        WRITE(LUPR1,*) 'Store NPOR=',STNPOR
        WRITE(LUPR1,*) 'Store HPOR=',STHPOR
      ENDIF
      LURAD=3
      LURAD=LUPR1
      JSPFLG    = .false.
      GBHMTX    = .false.
      IF(.NOT.NULLPR) THEN
        WRITE(LUPR1,*) 'GBHMTX=',GBHMTX
        WRITE(LUPR1,*) 'JSPECF=',JSPFLG
      ENDIF
C
C.... Loop through all patches to obtain patch numbers for radiative
C     patches. At present assumes that one patch for each thermal zone.
C     This calculation is repeated in GXS2SR.
C
      NTZ=0
      DO 100 IPAT=1,NUMPAT
C
C.... It is assumed that internal radiative thermal zones have
C     PATCH name commencing with '@RI' for internal radiation.
C
        IF(NAMPAT(IPAT)(1:3).EQ.'@RI') THEN
          NTZ = NTZ + 1
          F(L0RPNO+NTZ) = FLOAT(IPAT)
        ENDIF
 100  CONTINUE
C
C.... Coherence length of the radiation. Default value:
      GCOHL=1.E-03
C
      IF(CARTES) THEN
        IF(NX.EQ.1.OR.NY.EQ.1.OR.NZ.EQ.1) THEN
C-----------------------------------------------------------------------
C.... 1D or 2D cartesian viewfactors, begin                            |
C-----------------------------------------------------------------------
C.... GCOR1(I1=1,N1+1) and GCOR2(I2=1,N2+1) contain the cordinates of
C     faces in the '1' and '2' directions.
          IF(NZ.EQ.1) THEN
C.... If domain is x or y (1-d), or if domain is x-y (2-d)
            CALL SUB4(J2DDIR(2),1,J2DDIR(3),1,J2DDIR(4),2,J2DDIR(5),2)
            CALL SUB3(JN1,NX,JN2,NY,JM6,3)
            IF(.NOT.BFC) THEN
              DO 10101 I1=0,JN1
                DO 10101 I2=0,JN2
                  GCOR1(I1+1,I2+1) = F(L0F(LXYXU)+1+NY*(I1-1))
                  GCOR2(I1+1,I2+1) = F(L0F(LXYYV)+I2)
10101         CONTINUE
            ELSE
              DO 10111 I1=1,JN1+1
                DO 10111 I2=1,JN2+1
                  GCOR1(I1,I2) = XC(I1,I2,1)
                  GCOR2(I1,I2) = YC(I1,I2,1)
10111         CONTINUE
            ENDIF
          ELSEIF(NY.EQ.1) THEN
C.... if domain is x or z (1-d), or if domain is x-z (2-d)
            CALL SUB4(J2DDIR(2),1,J2DDIR(3),1,J2DDIR(6),2,J2DDIR(7),2)
            CALL SUB3(JN1,NX,JN2,NZ,JM6,2)
            IF(.NOT.BFC) THEN
              DO 10102 I1=0,JN1
                DO 10102 I2=0,JN2
                  GCOR1(I1+1,I2+1) = F(L0F(LXYXU)+I1)
                  GCOR2(I1+1,I2+1) = F(L0F(LZZW) +I2)
10102         CONTINUE
            ELSE
              DO 10112 I1=1,JN1+1
                DO 10112 I2=1,JN2+1
                  GCOR1(I1,I2) = XC(I1,1,I2)
                  GCOR2(I1,I2) = ZC(I1,1,I2)
10112         CONTINUE
            ENDIF
          ELSEIF(NX.EQ.1) THEN
C.... If domain is y or z (1-d), or if domain is y-z (2-d)
            CALL SUB4(J2DDIR(4),1,J2DDIR(5),1,J2DDIR(6),2,J2DDIR(7),2)
            CALL SUB3(JN1,NY,JN2,NZ,JM6,1)
            IF(.NOT.BFC) THEN
              DO 10103 I1=0,JN1
                DO 10103 I2=0,JN2
                  GCOR1(I1+1,I2+1) = F(L0F(LXYYV)+I1)
                  GCOR2(I1+1,I2+1) = F(L0F(LZZW) +I2)
10103         CONTINUE
            ELSE
              DO 10113 I1=1,JN1+1
                DO 10113 I2=1,JN2+1
                  GCOR1(I1,I2) = YC(1,I1,I2)
                  GCOR2(I1,I2) = ZC(1,I1,I2)
10113         CONTINUE
            ENDIF
          ENDIF
          IF(.NOT.BFC) THEN
            DO 10104 I1=1,JN1+1
10104         GCOR2(I1,1) = 0.0
            DO 10105 I2=1,JN2+1
10105         GCOR1(1,I2) = 0.0
          ENDIF
C
C....  1-d - length in 'second' dimension needs to be set to a
C      very large number for purpose of viewfactor calculation
          IF(JN1.EQ.1) THEN
            DO 10106 I2=1,JN2+1
10106       GCOR1(2,I2) = 1.E5
          ENDIF
          IF(JN2.EQ.1) THEN
            DO 10107 I1=1,JN1+1
10107       GCOR2(I1,2) = 1.E5
          ENDIF
C
          IF(STEMMS) THEN
            JMTCNT = 0
C.... Loop over patches to find surface information begins
            DO 1201 JTZ=1,NTZ
              IF(NZ.EQ.1) THEN
                CALL GETPTC(NAMPAT(NINT(F(L0RPNO+JTZ))),TYPE,
     1                      JJ1F,JJ1L,JJ2F,JJ2L,JDMY,JDMY,JDMY,JDMY)
              ELSEIF(NY.EQ.1) THEN
                CALL GETPTC(NAMPAT(NINT(F(L0RPNO+JTZ))),TYPE,
     1                      JJ1F,JJ1L,JDMY,JDMY,JJ2F,JJ2L,JDMY,JDMY)
              ELSEIF(NX.EQ.1) THEN
                CALL GETPTC(NAMPAT(NINT(F(L0RPNO+JTZ))),TYPE,
     1                      JDMY,JDMY,JJ1F,JJ1L,JJ2F,JJ2L,JDMY,JDMY)
              ENDIF
              JJTYPE = MOD(INT(TYPE),15)
C
C.... The patch NAMPAT(F(L0RPNO+JTZ)) contains the following CO and
C     VA for the variable EMMS1:
C
C     CO = material label of surface
C
C         | radiation temperature (different from surf.temp), IF GRDTMP > 0.0
C     VA =| 0.0, indicates that rad. temp. equals surf.temp , IF GRDTMP = 0.0
C         | thickness of the surface layer in micron        , IF GRDTMP < 0.0
C
              CALL GETCOV(NAMPAT(NINT(F(L0RPNO+JTZ))),
     1                    LBNAME('EMMS'),GMAT,GRDTMP)
              IF(GMAT.LT.100.) PRINT*,'SURFACE',JTZ,' CANNOT BE A GAS!'
              JMAT = INT(GMAT)
C
C.... This variable counts the materials present in the calculation
C     as a function of the material label
              IF(JMTCTR(JMAT).EQ.0) THEN
                JMTCNT = JMTCNT+1
                JMTCTR(JMAT) = JMTCNT
              ENDIF
C
C.... Store a surface layer thickness eventually present
C     (in unused part of array)
              SURFL =.FALSE.
              IF(GRDTMP.LT.0.0) THEN
                GSURF(JTZ,3) = ABS(GRDTMP)
                SURFL =.TRUE.
              ENDIF
C.... Store material label = GMAT (in unused part of array)
              GSURF(JTZ,2) = GMAT
C.... Store radiation temperature MAX(GRDTMP,F(L0FAVTS+JTZ))
C     in unused part of array
              IF(GRDTMP.GE.0.0) GSURF(JTZ,5) = GRDTMP
              GSURF(JTZ,5) = MAX(GRDTMP,300.)
C
C.... Calculation of min and max surface temperature TMPMIN and TMPMAX
C     of the present calculation, initialization:
              IF(JTZ.EQ.1) THEN
                TMPMIN = GSURF(JTZ,5)
                TMPMAX = GSURF(JTZ,5)
              ENDIF
C.... Update of TMPMIN and TMPMAX
              TMPMIN = MIN(TMPMIN,GSURF(JTZ,5))
              TMPMAX = MAX(TMPMAX,GSURF(JTZ,5))
C
C.... Next comes the calculation of absorption length inside the str.
C     medium for this calculation, move patch inside the str. medium
C     (except when the patch is already inside conductive material)
C
              IF(STPRPS) THEN
                L0PRPS = L0F(JPRPS)
                ICELL  = JJ2F+NY*(JJ1F-1)
                IF(NZ.NE.1) THEN
                  L0PRPS = L0F(ANYZ(JPRPS,JJ2F))
                  ICELL  = JJ1F
                ENDIF
                M1J = J2DDIR(JJTYPE)
                M2J = MOD(JJTYPE-1,2)+1
                JJTYPE = JJTYPE-2*MOD(JJTYPE,2)+1
                JJ1F = JJ1F+JSHFT1(M1J,M2J)
                JJ1L = JJ1L+JSHFT1(M1J,M2J)
                JJ2F = JJ2F+JSHFT2(M1J,M2J)
                JJ2L = JJ2L+JSHFT2(M1J,M2J)
                ICELL  = JJ2F+NY*(JJ1F-1)
                IF(NZ.NE.1) THEN
                  L0PRPS = L0F(ANYZ(JPRPS,JJ2F))
                  ICELL  = JJ1F
                ENDIF
              ENDIF
C
C.... In case of SURFL, find out whether it is LAYER ON SOLID or a THIN WALL
C
              IF(SURFL) THEN
                IMATC = 1000
                IF(JJ1F.EQ.0.OR.JJ1F.EQ.JN1+1
     1             .OR.JJ1F.EQ.0.OR.JJ1F.EQ.JN1+1) THEN
                  IMATC = 0
                ELSE
                  IMATC = F(L0PRPS+ICELL)
                ENDIF
                IF(IMATC.LT.100) THEN
C.... If Thinwall, store 1.0 in otherwise unused array
                  GRELP(JTZ,1) = 1.
C.... The material behind the surface is NOT SOLID. Hence the surface
C     defines a THIN WALL the absorption length must not be calculated,
C     but the patch must be transformed back
                  M1J = J2DDIR(JJTYPE)
                  M2J = MOD(JJTYPE-1,2)+1
                  JJTYPE = JJTYPE-2*MOD(JJTYPE,2)+1
                  JJ1F = JJ1F+JSHFT1(M1J,M2J)
                  JJ1L = JJ1L+JSHFT1(M1J,M2J)
                  JJ2F = JJ2F+JSHFT2(M1J,M2J)
                  JJ2L = JJ2L+JSHFT2(M1J,M2J)
                  ICELL  = JJ2F+NY*(JJ1F-1)
                  IF(NZ.NE.1) THEN
                    L0PRPS = L0F(ANYZ(JPRPS,JJ2F))
                    ICELL  = JJ1F
                  ENDIF
C.... The thickness of the thin wall gives the absorption length
                  IF(JMAT.GE.400) THEN
                    GLNGTH = ABS(GRDTMP)
                  ELSE
                    GLNGTH = 1.E+06
                  ENDIF
                  GO TO 1218
                ENDIF
              ENDIF
C
C.... Surface and volume material can be different. The volume material
C     is indicated by the value of PRPS in the cell. The calculation of
C     the absorption length GLNGTH is only done, when the surface
C     material is semitransparent and the volume material is
C     semitransparent. The calculation of the absorption length is done,
C     until the volume material changes or ends.
C
              IF(JMAT.LT.400) THEN
C.... Surface material is opaque, the absorption length is large compared
C     with the thickness.
                GLNGTH = 1.e+06
              ELSE
C.... Surface material is semitransparent
                IMATO = F(L0PRPS+ICELL)
C.... Volume material label is stored
                GSURF(JTZ,4) = IMATO
                IF(IMATO.NE.JMAT) THEN
                  IF(IMATO.LT.400) THEN
C.... Surface is semitransparent, but solid opaque.
C     No calculation of absorption length
                    GLNGTH = 1.e+06
                    GO TO 1218
                  ELSE
C.... Solid is semitransparent and of different material.
C     Check, whether a thickness of the surface layer was given.
                    IF(.NOT.SURFL) CALL SET_ERR(445,
     *                    'Error. THICKNESS OF LAYER NOT GIVEN.',1)
C                    IF(.NOT.SURFL) STOP 'THICKNESS OF LAYER NOT GIVEN '
                  ENDIF
                ENDIF
C.... Go from (jj1,jj2) inside the str. medium
                JJ1 = JJ1F+(JJ1L-JJ1F)/2
                JJ2 = JJ2F+(JJ2L-JJ2F)/2
                M1J = J2DDIR(JJTYPE)
                M2J = MOD(JJTYPE-1,2)+1
C.... The initial coordinates are
                B1 = GCOR1(JJ1+J2DB1(M1J,M2J,1,2),
     1                     JJ2+J2DB2(M1J,M2J,1,2))
                B2 = GCOR2(JJ1+J2DB1(M1J,M2J,1,2),
     1                     JJ2+J2DB2(M1J,M2J,1,2))
                M3 = (3+INT(SIGN(1.,GDIR(3-M1J))))/2
1219            J1OLD = JJ1
                J2OLD = JJ2
C.... Move to next cell
                JJ1 = JJ1+J2DSH1(M1J,M2J,M3,1)
                JJ2 = JJ2+J2DSH2(M1J,M2J,M3,1)
                IF(JJ1.NE.0.AND.JJ1.NE.JN1+1
     1             .AND.JJ2.NE.0.AND.JJ2.NE.JN2+1) THEN
                  ICELL  = JJ2+NY*(JJ1-1)
                  IF(NZ.NE.1) THEN
                    L0PRPS = L0F(ANYZ(JPRPS,JJ2))
                    ICELL  = JJ1
                  ENDIF
                  IMATC=F(L0PRPS+ICELL)
                  IF(IMATC.EQ.IMATO) GO TO 1219
                ENDIF
C.... End of str. material reached, calculate length of path
                A1 = GCOR1(J1OLD+J2DA1(M1J,M2J,1),
     1                     J2OLD+J2DA2(M1J,M2J,1))
                A2 = GCOR2(J1OLD+J2DA1(M1J,M2J,1),
     1                     J2OLD+J2DA2(M1J,M2J,1))
                GLNGTH = SQRT((A1-B1)**2+(A2-B2)**2)
              ENDIF
C
C.... Store the patch no. in EMMS in the following way to make accessible
C     the correct patchnumber avoiding the problem of double counting.
C
1218          DO 1215 J1=JJ1F,JJ1L
                IF(J1.EQ.0.OR.J1.EQ.JN1+1) GO TO 1215
                DO 1216 J2=JJ2F,JJ2L
                  IF(J2.EQ.0.OR.J2.EQ.JN2+1) GO TO 1216
C.... Direct f_array addressing of EMMS
                  L0EMMS = L0F(JEMMS)
                  ICELL  = J2+NY*(J1-1)
                  IF(NZ.NE.1) THEN
                    L0EMMS = L0F(ANYZ(JEMMS,J2))
                    ICELL  = J1
                  ENDIF
                  IF(QNE(F(L0EMMS+ICELL),0.0)) THEN
                    GRELP(INT(F(L0EMMS+ICELL)),JJTYPE) = JTZ
                    GRELP(         JTZ        ,JJTYPE) = JTZ
                  ELSE
                    GRELP(JTZ,JJTYPE) = JTZ
                    IF(.NOT.SURFL)
     1                 GRELP(JTZ,(JJTYPE-2*MOD(JJTYPE,2)+1))=-1.
C
C.... This indicates that a ray enters a region with EMMS(jj1,jj2)=0, but
C     the previous EMMS(j1old,j2old) referenced to a patchnumber -1.0,
C     the ray DID NOT leave a solid region, but is still in it.
C
                    F(L0EMMS+ICELL) = FLOAT(JTZ)
                  ENDIF
                  GSURF(INT(GRELP(INT(F(L0EMMS+ICELL)),JJTYPE)),1)=
     1                                                      GLNGTH
C
C.... The patchnumber belonging to a cell (j1,j2) and a face JJTYPE
C     of the cell (normal in the inside direction) can now be accessed
C     in the following way: look for the patchnumber in EMMS(J1,J2).
C     Now look in the array GRELP(EMMS(J1,J2),JJTYPE). This contains
C     the patchnumber of the proper patch, double counting is avoided.
C     We call this patchnumber 'THE RELATIVE PATCHNUMBER JJTYPE'.
C     The absorption length of patch JTZ is stored in GSURF(JTZ,1),
C     the radiation temperature was stored in GSURF(JTZ,5)
1216            CONTINUE
1215          CONTINUE
1201        CONTINUE
C.... loop over patches to find surface information ends
            IF(.NOT.NULLPR) WRITE(LURAD,*) 'EMMS'
            DO 1205 I=1,JN1
              IF(.NOT.NULLPR) WRITE(LURAD,*)
     1              (F(L0EMMS+J+NY*(I-1)),J=1,JN2)
1205        CONTINUE
C
C.... JMTLBL is the material label, extending the PHOENICS convension,
C     fluid:          0
      SUBROUTINE GSPECT(JL,JBAND,JINT,JSTRNT,JMTCTR,
     1              JSPFLG,JMTINC,TMPMIN,TMPMAX,OPT)
C
C.... Maximal number of viewfactor matrices = maximal
C             number of spectral bands for str. media
      PARAMETER (JLMX=5)
C.... Number of spectral intervals
      PARAMETER (JSPMX=60)
C.... Maximal number of different materials in this calculation
C     (to limit the storage for optical data)
      PARAMETER (JMTLIM=10)
C.... Range of material numbers in data file
      PARAMETER (JMTIND=500)
      INTEGER JINT(JSPMX),JSTRNT(JSPMX)
      INTEGER JMTCTR(JMTIND)
      REAL OPT(JMTLIM,JSPMX,9)
C.... a new array providing the spectral weight for every
C     surface must be introduced
      LOGICAL RDFLAG,JSPFLG
      CHARACTER*9 GLINE
      CHARACTER*80 BUFF
      COMMON /LDATA/ LDAT1(31),NULLPR,LDAT2(30),DISTIL,LDAT3(21)
      LOGICAL LDAT1,NULLPR,LDAT2,DISTIL,LDAT3
C----------------------------------------------------------------------
C     DATA IN THE DATA FILE 'optical.data'
C
C     Data file contains the complex refractive index:
C        n(v,T) + i k(v,T)
C     Spectral dependence is represented in JSPMX=60 intervals:
C        the spectrum in the range 2=log10( 10**2 *cm ) to
C                                  5=log10( 10**5 *cm ) is equally divided
C        into 60 spectral intervals
C     Temperature dependence is contained in the power series representation:
C        n(i,T) = ni(0) + ni(1)*((T-300)/1000) +
C                         ni(2)*((T-300)/1000)**2 +
C                         ni(3)*((T-300)/1000)**3
C        k(i,T) = ki(0) + ki(1)*((T-300)/1000) +
C                         ki(2)*((T-300)/1000)**2 +
C                         ki(3)*((T-300)/1000)**3
C
C     1. line:  # comment
C     2. line:  materiallabel, const. EMMSsivity, const. reflectivity,
C        label of reflection law, param. 1, param. 2, param. 3, param. 4
C        number of band (default), n(0), n(1), n(2), n(3),
C                                  k(0), k(1), k(2), k(3),
C                               interval number (not used)
C     3. line:                  data for  1. spectral interval
C         ...                                            ...
C    62. line:                  data for 60. spectral interval
C
C                 NEXT MATERIAL
C    63. line:  # comment
C    64. line:  materiallabel, label of reflection law, parameter 1,
C                             parameter 2, parameter 3, parameter 4
C    65. line:                  data for  1. spectral interval
C         ...                                            ...
C   124. line:                  data for 60. spectral interval
C                 NEXT MATERIAL
C                  ......
C-----------------------------------------------------------------------
C     Read only the data of materials, which are present in the actual
C     problem, into the array OPT( counter,spectralinterval,property )
      OPEN(2,FILE='optical.data')
C
      lurad=14
      JMTINC = 0
      DO 10 J=1,JMTIND
        RDFLAG =.FALSE.
        READ(2,'(a)',END=11) GLINE
        IF(GLINE.EQ.'endofdata') GO TO 11
        READ(2,*) JMTLBL,GECNST,GRCNST
C
      DO 20 JMAT=1,JMTIND
        IF(JMTCTR(JMAT).NE.0.AND.JMAT.EQ.JMTLBL) RDFLAG =.TRUE.
20    CONTINUE
C
      IF(RDFLAG) THEN
        IF(JSPFLG) THEN
          JMTINC = JMTINC + 1
          IF(JMTINC.GT.JMTLIM) CALL SET_ERR(449,
     *          'Error. Too many materials.',1)
C          IF(JMTINC.GT.JMTLIM) STOP 'too many materials!'
          DO 30 JIN=1,JSPMX
            READ(2,*) (OPT(JMTCTR(JMTLBL),JIN,JPROP),JPROP=1,9)
            JDMY = INT(OPT(JMTCTR(JMTLBL),JIN,1))
            IF(JMTLBL.GE.400) THEN
              JINT(JIN)   = JINT(JIN)  +JDMY
              JSTRNT(JIN) = JSTRNT(JIN)+JDMY
            ELSE
              JINT(JIN)   = JINT(JIN)  +JDMY
            ENDIF
30        CONTINUE
        ELSE
          OPT(JMTCTR(JMTLBL),1,1) = GECNST
          OPT(JMTCTR(JMTLBL),1,2) = GRCNST
          DO 40 JIN=1,JSPMX
            READ(2,*)
40        CONTINUE
        ENDIF
      ELSE
        DO 41 JIN=1,JSPMX
          READ(2,*)
41      CONTINUE
      ENDIF
10    CONTINUE
11    CLOSE(2)
C
C.... Calculate the lowest spectral interval from the lowest temperature
C     in the problem, and the highest spectral interval from the highest
C     temperature in the problem
      JDMY = INT(20*(LOG10(TMPMIN/0.3668)-2.9))
      IF(JDMY.LT.0) THEN
        WRITE(BUFF,'(A,1X,1PE10.3)') 'MINIMAL TEMPERATURE=',TMPMIN
        CALL PUT_LINE(BUFF,.false.)
        CALL PUT_LINE('WARNING, MINIMAL TEMPERATURE BELOW LIMIT',.true.)
      ENDIF
      INTMIN = MAX(0,JDMY)+1
      JDMY = INT(20*(LOG10(TMPMAX/0.3668)-1.4))
      IF(JDMY.GT.59) THEN
        WRITE(BUFF,'(A,1X,1PE10.3)') 'MAXIMAL TEMPERATURE=',TMPMAX
        CALL PUT_LINE(BUFF,.false.)
        CALL PUT_LINE('WARNING, MAXIMAL TEMPERATURE ABOVE LIMIT',.true.)
      ENDIF
      INTMAX = MAX(0,JDMY)+1
C
C.... Calculate the bandstructure from the default values from
C     'optical.data' the bands for str. materials and for the whole
C     calculation must be determined.
C.... JINT(I) contains the band number for the spectral interval I
C     JSTRNT(I) contains the band number for the spectral interval
C     I of a str. medium
C
      DO 1301 I=1,JSPMX-1
        JINT(I)   = MIN(1,JINT(I+1)-JINT(I))*
     1              MAX(0,SIGN(1,I-INTMIN))*MAX(0,SIGN(1,INTMAX-I-1))
        JSTRNT(I) = MIN(1,JSTRNT(I+1)-JSTRNT(I))*
     1              MAX(0,SIGN(1,I-INTMIN))*MAX(0,SIGN(1,INTMAX-I-1))
1301  CONTINUE
      JINT(JSPMX)   = 0
      JSTRNT(JSPMX) = 0
      DO 1302 I=INTMAX,INTMIN+1,-1
        JINT(I)   = JINT(I-1)
        JSTRNT(I) = JSTRNT(I-1)
1302  CONTINUE
      JINT(INTMIN)   = 1
      JSTRNT(INTMIN) = 1
      DO 1303 I=INTMIN+1,INTMAX
        JINT(I)   = JINT(I-1)+JINT(I)
        JSTRNT(I) = JSTRNT(I-1)+JSTRNT(I)
1303  CONTINUE
C.... Simplification, if there is no spectral dependence required
      IF(.NOT.JSPFLG) THEN
      DO 1304 I=1,JSPMX
        JINT(I)   = 1
        JSTRNT(I) = 1
1304  CONTINUE
      JL = 1
      JBAND = 1
      ENDIF
C.... Number of spectral bands for str. medium and hence number of
C     viewfactormatrices
      JL = JSTRNT(INTMAX)
C.... Number of spectral bands for the surfaces and hence for the
C     radiation exchange
      JBAND = JINT(INTMAX)
C
      IF(JL.GT.JLMX) CALL SET_ERR(450,
     *     'Error. Increase JLMX to allow more view factor matrices.',1)
C      IF(JL.GT.JLMX)STOP'increase JLMX to allow more viewfactormatrices'
C
      IF(.NOT.NULLPR) THEN
        WRITE(LURAD,*)
     1 '&&&&&& output from the spectral calculations &&&&&&'
     1,' SPECTRAL DEPENDENCE=',JSPFLG
        WRITE(LURAD,*) JMTINC,' different materials in this calcul.n'
        WRITE(LURAD,*) JL,' bands for str. materials'
        WRITE(LURAD,*) JBAND,' bands for all materials'
        WRITE(LURAD,*) 'the band-structure for str. materials is:'
        WRITE(LURAD,*) (JSTRNT(J),J=1,JSPMX)
        WRITE(LURAD,*) 'the band-structure for all materials is:'
        WRITE(LURAD,*) (JINT(J),J=1,JSPMX)
      ENDIF
C
      END
C----------------------------------------------------------------------
c
      SUBROUTINE GBAND(NTZ,JL,JBAND,JINT,JSTRNT,JMTCTR,GCOHL,
     1                 JSPFLG,SURFL,GRELP,GSURF,GR,GE,GW,OPT)
      COMMON /LDATA/ LDAT1(31),NULLPR,LDAT2(30),DISTIL,LDAT3(21)
      LOGICAL LDAT1,NULLPR,LDAT2,DISTIL,LDAT3
C
C.... Maximal number of radiating patches
      PARAMETER (JPTHMX=50)
C.... Number of spectral intervals
      PARAMETER (JSPMX=60)
C.... Number of spectral intervals
      PARAMETER (JBNDMX=7)
C.... Maximal number of different materials in this
C     calculation (to limit the storage for optical data)
      PARAMETER (JMTLIM=10)
C.... Range of material numbers in data file
      PARAMETER (JMTIND=500)
C
      INTEGER JINT(JSPMX),JSTRNT(JSPMX)
      INTEGER JMTCTR(JMTIND)
      REAL OPT( JMTLIM , JSPMX , 9 )
      REAL GRELP(JPTHMX,JBNDMX),GSURF(JPTHMX,JBNDMX)
      REAL GE(JPTHMX,JSPMX),GR(JPTHMX,JSPMX)
C.... a new array providing the spectral weight for every surface
C     must be introduced
      REAL GW(JPTHMX,JSPMX)
      LOGICAL JSPFLG,SURFL,QNE,QEQ,QGT
C.... unit: cm K
      DATA CC2 / 1.438769 /
C-----------------------------------------------------------------------
C.... Patchwise arrays GRELP(NTZ,JBNDMX) and GSURF(NTZ,JBNDMX)
C     contain the following information:
C
C.... The last entries of the array GE:
C     GRELP(JTZ,7) = relative patch# from direction 7, JTZ is value of EMMS
C     GRELP(JTZ,6) = relative patch# from direction 6, JTZ is value of EMMS
C     GRELP(JTZ,5) = relative patch# from direction 5, JTZ is value of EMMS
C     GRELP(JTZ,4) = relative patch# from direction 4, JTZ is value of EMMS
C     GRELP(JTZ,3) = relative patch# from direction 3, JTZ is value of EMMS
C     GRELP(JTZ,2) = relative patch# from direction 2, JTZ is value of EMMS
C     GRELP(JTZ,1) = 1. for thinwall, 0. otherwise
C
C.... The first entries of the array GE:
C     GE(JTZ,1)     = the EMMSsivity of the 1. spectral band
C     GE(JTZ,2)     = the EMMSsivity of the 2. spectral band
C       ...
C     GE(JTZ,JBAND) = the EMMSsivity of the last spectral band
C
C.... The last entries of the array GR:
C     GSURF(JTZ,4) = the material label of material below layer
C     GSURF(JTZ,3) = eventually the thickness of a surface layer
C     GSURF(JTZ,2) = the material label
C     GSURF(JTZ,1) = the absorption length of the radiation in the patch
C     GSURF(JTZ,5)   = the radiation temperature of the patch
C
C.... The first entries of the array GE:
C     GR(JTZ,1)     = the reflectivity of the 1. spectral band
C     GR(JTZ,2)     = the reflectivity of the 2. spectral band
C       ...
C     GR(JTZ,JBAND) = the reflectivity of the last spectral band
C
C.... Since the number of spectral bands will be much smaller than
C     the number of spectral intervals, there will be no overwriting.
C-----------------------------------------------------------------------
C
      lurad=14
      IF(.NOT.NULLPR) WRITE(LURAD,*)
     1 '&&& information contained in the patchwise arrays &&&'
C
C.... Loop over patches to determine surface properties of the bands starts
      DO 1202 JTZ=1,NTZ
C     GE(JTZ,JBAND) contains the EMMSsivity of patch JTZ in the band JBAND,
C     GR(JTZ,JBAND) contains the reflectivity of patch JTZ in the band JBAND
C     1-GR-GE is then the transmissivity of patch JTZ in the band JBAND
C
        IF(.NOT.NULLPR)
     1    WRITE(LURAD,*)'&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& patch ',
     1    jtz,' &&'
C
        IF(JSPFLG) THEN
C.... Spectral dependence activated
C
          JBANDO = 1
          DO 1212 JS=1,JSPMX
            IF(JINT(JS).EQ.0) GO TO 1212
C
C.... Establish ref. from the label of the spectr. band to the label of the
C     view-factor matrix
            JINTJS = JINT(JS)
            JSTRNT(JINTJS) = JSTRNT(JS)
C
            IF(JINTJS.NE.JBANDO) THEN
C.... finish interval averaging of the previous band by normalization with
C     energy fraction in band
              GE(JTZ,JBANDO) = GE(JTZ,JBANDO)/GW(JTZ,JBANDO)
              GR(JTZ,JBANDO) = GR(JTZ,JBANDO)/GW(JTZ,JBANDO)
              GW(JTZ,JINTJS) = 0.0
              GE(JTZ,JINTJS) = 0.0
              GR(JTZ,JINTJS) = 0.0
            ENDIF
C
C.... GWAVNO is the wavenumber of the upper end of the interval
C     mult. of GWAVNO with 10**(-1.0/20.)=0.891250938
C     gives the waven. of the lower end of the interval
C     mult. of GWAVNO with 10**(-0.5/20.)=0.944060876
C     gives the waven. of the middle of the interval
C     PLINT is the fraction of radiation in the spectral interval
            GWAVNO = 10**(2.+FLOAT(JS)/20.)
            GT = GSURF(JTZ,5)
            PLINT = FRACT(CC2*GWAVNO*0.891250938/GT)-
     1              FRACT(CC2*GWAVNO/GT)
            GWAVNO = GWAVNO*0.944060876
C
C.... extract n + i k of bulk from the coefficients in OPT(JM,JS,JPROP)
            JM = JMTCTR(INT(GSURF(JTZ,2)))
            GN = OPT(JM,JS,2)+OPT(JM,JS,3)*1.E-03*(GT-300.)
     1        +OPT(JM,JS,4)*1.E-06*(GT-300.)**2
     1        +OPT(JM,JS,5)*1.E-09*(GT-300.)**3
            GK = OPT(JM,JS,6)+OPT(JM,JS,7)*1.E-03*(GT-300.)
     1        +OPT(JM,JS,8)*1.E-06*(GT-300.)**2
     1        +OPT(JM,JS,9)*1.E-09*(GT-300.)**3
            GLNGTH=GSURF(JTZ,1)
            IF(INT(GSURF(JTZ,2)).LT.400) GLNGTH=1.E+06
            GLAYER = 0.0
            IF(SURFL.AND.QNE(GRELP(JTZ,1),1.0)) THEN
C.... (when GRELP(JTZ,1) = 1. the surface is not a SURFL but a THINWALL)
C     extract ns + i ks of layer from the coefficients in OPT(JM,JS,JPROP)
              JMS = JMTCTR(INT(GSURF(JTZ,4)))
              GNS = OPT(JMS,JS,2)+OPT(JMS,JS,3)*1.E-03*(GT-300.)
     1         +OPT(JMS,JS,4)*1.E-06*(GT-300.)**2
     1         +OPT(JMS,JS,5)*1.E-09*(GT-300.)**3
              GKS = OPT(JMS,JS,6)+OPT(JMS,JS,7)*1.E-03*(GT-300.)
     1         +OPT(JMS,JS,8)*1.E-06*(GT-300.)**2
     1         +OPT(JMS,JS,9)*1.E-09*(GT-300.)**3
              GLAYER=GSURF(JTZ,3)
            ENDIF
C
C.... Calculation of TRANSMITTANCE, REFLECTANCE and EMITTANCE when the
C     coherence lenght of the radiation COHL (specified in the beginning
C     of GX2DVF, the parameter should be 0.1 mm be default) is larger
C     than the wavelength of radiation, interference effects in the layer
C     (if one is present) are included. Normal incidence of radiation is
C     assumed here.
            CALL COATNG(0.0,GWAVNO,GCOHL,GLAYER,GLNGTH,
     1                      GN,GK,GNS,GKS,GEM,GRES)
C
            IF(.NOT.NULLPR) THEN
              WRITE(LURAD,*)'** band=',jintjs,', mat=',GSURF(JTZ,2)
     1        ,', T=',gt,', v[1/cm]=',gwavno,', n=',gn,', k=',gk
            ENDIF
C
C.... Sum the properties of the interval to the band properties
            GW(JTZ,JINTJS) = GW(JTZ,JINTJS)+PLINT
            GE(JTZ,JINTJS) = GE(JTZ,JINTJS)+GEM*PLINT
            GR(JTZ,JINTJS) = GR(JTZ,JINTJS)+GRES*PLINT
C.... Keep the bandnumber for the final normalization
            JBANDO = JINTJS
1212      CONTINUE
C.... Normalization of the last band must not be forgotten
          GE(JTZ,JBANDO) = GE(JTZ,JBANDO)/GW(JTZ,JBANDO)
          GR(JTZ,JBANDO) = GR(JTZ,JBANDO)/GW(JTZ,JBANDO)
        ELSE
C.... No spectral dependence activated
          GW(JTZ,1) = 1.0
          GE(JTZ,1) = OPT(JMTCTR(INT(GSURF(JTZ,2))),1,1)
          GR(JTZ,1) = OPT(JMTCTR(INT(GSURF(JTZ,2))),1,2)
        ENDIF
        IF(.NOT.NULLPR) THEN
          IF(INT(GSURF(JTZ,2)).GE.100) THEN
            WRITE(LURAD,*) '****** WALL:'
            WRITE(LURAD,*) 'the absorption length=',GSURF(JTZ,1)
            WRITE(LURAD,*) 'relative patch# 7=',GRELP(JTZ,7)
            WRITE(LURAD,*) 'relative patch# 6=',GRELP(JTZ,6)
            WRITE(LURAD,*) 'relative patch# 5=',GRELP(JTZ,5)
            WRITE(LURAD,*) 'relative patch# 4=',GRELP(JTZ,4)
            WRITE(LURAD,*) 'relative patch# 3=',GRELP(JTZ,3)
            WRITE(LURAD,*) 'relative patch# 2=',GRELP(JTZ,2)
            IF(QGT(GSURF(JTZ,3),0.0)) WRITE(LURAD,*) 'SURFL'
            IF(QEQ(GRELP(JTZ,1),1.0)) WRITE(LURAD,*) 'THINWALL'
          ELSE
            WRITE(LURAD,*) 'SOLID WALL:'
          ENDIF
          WRITE(LURAD,*) '****** band dependent properties'
          WRITE(LURAD,*) 'PL(average)   ',(GW(JTZ,JS),JS=1,JBAND)
          WRITE(LURAD,*) 'EM(average)   ',(GE(JTZ,JS),JS=1,JBAND)
          WRITE(LURAD,*) 'RES(average)  ',(GR(JTZ,JS),JS=1,JBAND)
          WRITE(LURAD,*) 'TRS(average)  ',
     1                      (1.-GR(JTZ,JS)-GE(JTZ,JS),JS=1,JBAND)
        ENDIF
C
1202  CONTINUE
C
      IF(.NOT.NULLPR) WRITE(LURAD,*)
     1 '&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&'
C
      END
C----------------------------------------------------------------------
      FUNCTION FRACT(ARG)
C----------------------------------------------------------------------
C.... fract is the energy of black body radiation
C     from 0 to V = CC2*GWAVNO/GTEMP.
C----------------------------------------------------------------------
      REAL(8) PI,V,ACT,SM,F1
      DATA PI/3.141592654D0/
      V=DBLE(ARG)
      IF(V.GE.2.0D0) THEN
        ACT=0.
        DO 5 M=1,100
          SM=DBLE(M)
          F1=(DEXP(-SM*V)/SM**4)*(((SM*V+3.0D0)*SM*V+6.0D0)*SM*V+6.0D0)
          ACT=ACT+F1
          IF(DABS(F1).LE.1.D-07) GO TO 10
5       CONTINUE
10      ACT=ACT*(15.D0/PI**4)
      ELSE
        ACT = 1.0D0-(15.D0/PI**4)*V**3*(1.D0/3.D0-V/8.D0+V**2/60.D0-
     1                 V**4/5040.D0+V**6/272160.D0-V**8/13305600.D0)
      ENDIF
      FRACT=SNGL(ACT)
      END
C-----------------------------------------------------------------------
c
      SUBROUTINE COATNG(GSNTH,GWAVNO,GCOHL,GLAYER,GLNGTH,
     1                            GN,GK,GNS,GKS,GEM,GRES)
C
      COMPLEX I,AMPLTD,D2,D3,RS,RP,TS,TP
      COMPLEX RS1,RS2,RS3,TS1,TS2,TS3,RP1,RP2,RP3,TP1,TP2,TP3
      COMPLEX CD2,CD3,SALP,SBET,SRHO,SSIG,SEPS,PALP,PBET,PRHO,PSIG,PEPS
      COMPLEX G1CMPX,G2CMPX,G3CMPX
      COMPLEX GSNTH1,GSNTH2,GSNTH3,GCSTH1,GCSTH2,GCSTH3
      LOGICAL QEQ
      DATA PI /3.14159265/
      I=CMPLX( 0.0 , 1.0 )
C
      IF(QEQ(GLAYER,0.0)) THEN
C.... Complex refractive index
        G1CMPX=CMPLX(1.0,0.0)
        G2CMPX=CMPLX(GN,-GK)
        GSNTH1=CMPLX( GSNTH , 0.0 )
C.... Snell's law
        GSNTH2 = GSNTH1*G1CMPX/G2CMPX
        GCSTH1 = CSQRT(1.0-GSNTH1**2)
        GCSTH2 = CSQRT(1.0-GSNTH2**2)
C.... Assume planeparallel layer of gas(=1) | bulk(=2) | gas(=1)
        RS1  = AMPLTD(1,G1CMPX,G2CMPX,GCSTH1,GCSTH2)
        RP1  = AMPLTD(3,G1CMPX,G2CMPX,GCSTH1,GCSTH2)
        RS12 = CABS(RS1)**2
        RP12 = CABS(RP1)**2
        T2   = EXP(-4.*PI*GK*GLNGTH*GWAVNO/REAL(GCSTH2))
C.... S-polarized
        R123 = 1.0-(RS12*T2)**2
        RS   = RS12+((1.-RS12)*T2)**2*RS12/R123
        TS   = T2*(1.0-RS12)**2/R123
C.... P-polarized
        R123 = 1.0-(RP12*T2)**2
        RP   = RP12+((1.-RP12)*T2)**2*RP12/R123
        TP   = T2*(1.0-RP12)**2/R123
      ELSE
C.... Complex refractive index
        G1CMPX = CMPLX(1.0,0.0)
        G2CMPX = CMPLX(GNS,-GKS)
        G3CMPX = CMPLX(GN,-GK)
        GSNTH1 = CMPLX(GSNTH,0.0)
C.... Snell's law
        GSNTH2 = GSNTH1*G1CMPX/G2CMPX
        GSNTH3 = GSNTH2*G2CMPX/G3CMPX
        GCSTH1 = CSQRT(1.0-GSNTH1**2)
        GCSTH2 = CSQRT(1.0-GSNTH2**2)
        GCSTH3 = CSQRT(1.0-GSNTH3**2)
C.... Complex phases
        D2 = 2.0*PI*GLAYER*GWAVNO*GCSTH2*G2CMPX
        D3 = 2.0*PI*GLNGTH*GWAVNO*GCSTH3*G3CMPX
C.... Complex amplitudes: assume planeparallel layer of
C       gas(=1) | coating(=2) | bulk(=3) | gas(=1)
        RS1 = AMPLTD(1,G1CMPX,G2CMPX,GCSTH1,GCSTH2)
        RS2 = AMPLTD(1,G2CMPX,G3CMPX,GCSTH2,GCSTH3)
        RS3 = AMPLTD(1,G3CMPX,G1CMPX,GCSTH3,GCSTH1)
        RP1 = AMPLTD(3,G1CMPX,G2CMPX,GCSTH1,GCSTH2)
        RP2 = AMPLTD(3,G2CMPX,G3CMPX,GCSTH2,GCSTH3)
        RP3 = AMPLTD(3,G3CMPX,G1CMPX,GCSTH3,GCSTH1)
        TS1 = AMPLTD(2,G1CMPX,G2CMPX,GCSTH1,GCSTH2)
        TS2 = AMPLTD(2,G2CMPX,G3CMPX,GCSTH2,GCSTH3)
        TS3 = AMPLTD(2,G3CMPX,G1CMPX,GCSTH3,GCSTH1)
        TP1 = AMPLTD(4,G1CMPX,G2CMPX,GCSTH1,GCSTH2)
        TP2 = AMPLTD(4,G2CMPX,G3CMPX,GCSTH2,GCSTH3)
        TP3 = AMPLTD(4,G3CMPX,G1CMPX,GCSTH3,GCSTH1)
C
        IF(GCOHL.GT.GLAYER) THEN
C-----------------------------------------------------------------------
C.... When coherence length greater than layer thickness,
C     interference in layer
          CD2  = CEXP(-I*D2)
          CD3  = CEXP(-I*D3)
          RD3  = CABS(CD3)
C
          SALP = RS1+RS2*CD2
          SBET = RS3*(CD2+RS1*RS2)*RD3
          SRHO = 1.0+RS1*RS2*CD2
          SSIG = RS1*RS3*(1.+CD2)*RD3
          SEPS = TS1*TS2*TS3*CD2*CD3
          PALP = RP1+RP2*CD2
          PBET = RP3*(CD2+RP1*RP2)*RD3
          PRHO = 1.0+RP1*RP2*CD2
          PSIG = RP1*RP3*(1.+CD2)*RD3
          PEPS = TP1*TP2*TP3*CD2*CD3
C.... Squared amplitudes, averaged over the phase in the bulk
          IF(CABS(SRHO).LE.CABS(SSIG)) CALL SET_ERR(451,
     *                'Error. WRONG INTEGRAL, SET COHL=0.',1)
C          IF(CABS(SRHO).LE.CABS(SSIG))STOP 'WRONG INTEGRAL, SET COHL=0.'
          RS=(CABS(SALP)**2+CABS(SBET)**2)/(CABS(SRHO)**2-CABS(SSIG)**2)
     1     -2.*REAL(SALP*CONJG(SBET)*SRHO*CONJG(SSIG))/
     1     ((CABS(SRHO)**2-CABS(SSIG)**2)*CABS(SRHO)**2)
          TS=CABS(SEPS)**2/(CABS(SRHO)**2-CABS(SSIG)**2)
C
          IF(CABS(PRHO).LE.CABS(PSIG)) CALL SET_ERR(452,
     *                'Error. WRONG INTEGRAL, SET COHL=0.',1)
C          IF(CABS(PRHO).LE.CABS(PSIG))STOP 'WRONG INTEGRAL, SET COHL=0.'
          RP=(CABS(PALP)**2+CABS(PBET)**2)/(CABS(PRHO)**2-CABS(PSIG)**2)
     1     -2.*REAL(PALP*CONJG(PBET)*SRHO*CONJG(PSIG))/
     1     ((CABS(PRHO)**2-CABS(PSIG)**2)*CABS(PRHO)**2)
          TP=CABS(PEPS)**2/(CABS(PRHO)**2-CABS(PSIG)**2)
        ELSE
C-----------------------------------------------------------------------
C.... Otherwise no interference in layer and bulk
          RS12 = CABS(RS1)**2
          RS23 = CABS(RS2)**2
          RS31 = CABS(RS3)**2
          RP12 = CABS(RP1)**2
          RP23 = CABS(RP2)**2
          RP31 = CABS(RP3)**2
          T2   = EXP(-4.*PI*GKS*GLAYER*GWAVNO/REAL(GCSTH2))
          T3   = EXP(-4.*PI*GK *GLNGTH*GWAVNO/REAL(GCSTH3))
C.... S-polarized
          R123 = 1.0-RS12*T2**2*RS23
          R231 = 1.0-RS23*T3**2*RS31
          T23  = T2*(1.0-RS23)*T3
          RS   = RS12 + (1.0-RS12)**2*T2**2*RS23/R123 +
     1    (1.-RS12)**2*T23**2*RS31/(R123*(R123*R231-RS12*T23**2*RS31))
          TS=(1.-RS12)*T23*(1.-RS31)/(R123*(R123*R231-RS12*T23**2*RS31))
C.... P-polarized
          R123 = 1.0-RP12*T2**2*RP23
          R231 = 1.0-RP23*T3**2*RP31
          T23  = T2*(1.0-RP23)*T3
          RP   = RP12 + (1.-RP12)**2*T2**2*RP23/R123 +
     1    (1.-RP12)**2*T23**2*RP31/(R123*(R123*R231-RP12*T23**2*RP31))
          TP=(1.-RP12)*T23*(1.-RP31)/(R123*(R123*R231-RP12*T23**2*RP31))
        ENDIF
      ENDIF
C
      GRES = 0.5*(RS+RP)
      GTES = 0.5*(TS+TP)
      GEM  = 1.0-GRES-GTES
C
      END
C-----------------------------------------------------------------------
      COMPLEX FUNCTION AMPLTD(IA,G1CMPX,G2CMPX,GCSTH1,GCSTH2)
      COMPLEX G1CMPX,G2CMPX,GCSTH1,GCSTH2
C
C.... IA=1: r_s = reflectivity s-polarized
C     IA=2: t_s = transmissivity s-polarized
C     IA=3: r_p = reflectivity p-polarized
C     IA=4: t_p = transmissivity p-polarized
C
      IF(IA.EQ.1) THEN
        AMPLTD = (GCSTH1*G1CMPX-GCSTH2*G2CMPX)/
     1           (GCSTH1*G1CMPX+GCSTH2*G2CMPX)
      ELSEIF(IA.EQ.2) THEN
        AMPLTD = 2.0*GCSTH1*G1CMPX/(GCSTH1*G1CMPX+GCSTH2*G2CMPX)
      ELSEIF(IA.EQ.3) THEN
        AMPLTD = (GCSTH1*G2CMPX-GCSTH2*G1CMPX)/
     1           (GCSTH1*G2CMPX+GCSTH2*G1CMPX)
      ELSEIF(IA.EQ.4) THEN
        AMPLTD = 2.0*GCSTH1*G1CMPX/(GCSTH2*G1CMPX+GCSTH1*G2CMPX)
      ENDIF
      END
C-----------------------------------------------------------------------
c
      SUBROUTINE GAUSSJ(A,N,NP,B,M,MP)
C.... From W H Press, B P Flannery, S A Teukolsky, and W T Vetterling:
C     'Numerical Recipes in FORTRAN'
      PARAMETER (JPTHMX=50)
      DIMENSION A(JPTHMX,JPTHMX),B(JPTHMX,JPTHMX)
      DIMENSION IPIV(JPTHMX),INDXR(JPTHMX),INDXC(JPTHMX)
      LOGICAL QEQ
      DO 11 J=1,N
        IPIV(J)=0
11    CONTINUE
      DO 22 I=1,N
        BIG=0.
        DO 13 J=1,N
          IF(IPIV(J).NE.1) THEN
            DO 12 K=1,N
              IF(IPIV(K).EQ.0) THEN
                IF(ABS(A(J,K)).GE.BIG) THEN
                  BIG=ABS(A(J,K))
                  IROW=J
                  ICOL=K
                ENDIF
              ELSEIF(IPIV(K).GT.1) THEN
                GOTO 999
              ENDIF
12          CONTINUE
          ENDIF
13      CONTINUE
        IPIV(ICOL)=IPIV(ICOL)+1
        IF(IROW.NE.ICOL) THEN
          DO 14 L=1,N
            DUM=A(IROW,L)
            A(IROW,L)=A(ICOL,L)
            A(ICOL,L)=DUM
14        CONTINUE
          DO 15 L=1,M
            DUM=B(IROW,L)
            B(IROW,L)=B(ICOL,L)
            B(ICOL,L)=DUM
15        CONTINUE
        ENDIF
        INDXR(I)=IROW
        INDXC(I)=ICOL
        IF(QEQ(A(ICOL,ICOL),0.0)) GOTO 999
        PIVINV=1./A(ICOL,ICOL)
        A(ICOL,ICOL)=1.
        DO 16 L=1,N
          A(ICOL,L)=A(ICOL,L)*PIVINV
16      CONTINUE
        DO 17 L=1,M
          B(ICOL,L)=B(ICOL,L)*PIVINV
17      CONTINUE
        DO 21 LL=1,N
          IF(LL.NE.ICOL) THEN
            DUM=A(LL,ICOL)
            A(LL,ICOL)=0.
            DO 18 L=1,N
              A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
18          CONTINUE
            DO 19 L=1,M
              B(LL,L)=B(LL,L)-B(ICOL,L)*DUM
19          CONTINUE
          ENDIF
21      CONTINUE
22    CONTINUE
      DO 24 L=N,1,-1
        IF(INDXR(L).NE.INDXC(L)) THEN
          DO 23 K=1,N
            DUM=A(K,INDXR(L))
            A(K,INDXR(L))=A(K,INDXC(L))
            A(K,INDXC(L))=DUM
23        CONTINUE
        ENDIF
24    CONTINUE
      RETURN
  999 CALL WRTCHS('Singular Matrix in GAUSSJ in gx2dvf.for')
        CALL SET_ERR(297,
     *       'Error. Singular Matrix in GAUSSJ in gx2dvf.for.',3)
C      CALL EAROUT(3)
      END
c