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