#$r002
   PHOTON USE
   p
 
 
 
 
 
   view x
   msg temperature contours
   con h1 x 1 fi;.001;gr x 1
   pause
   msg Press  and then  to END
   pause
 
   ENDUSE
TEXT(BFC RADIATION+HEAT SOURCE IN A SLAB:125
TITLE
  DISPLAY
   The problem considered is same as case 122 except that it is
   performed on a non-orthogonal BFC mesh.
  ENDDIS
REAL(EMWN,TWN,EWN,TA,GY,EG,EGC,EGW,ACON)
REAL(QRAD,QVOL,KRAD,OTHICK,LENGTH,QRADA,TWNA,QWNA,TGCA,TGWA)
CHAR(CH1);INTEGER(JJM1)
MESG( Enter optical thickness  0 < Kr*L << 10.0 (default 5)
READVDU(OTHICK,REAL,5.0)
LENGTH=1.0;KRAD=OTHICK/LENGTH
SCATT=0.0;EMISS=KRAD-SCATT
EMWN=1.0;TWN=1000.0;EWN=SIGMA*TWN**4;QRAD=EWN;QVOL=EWN/LENGTH
  ** analytical solution
QWNA=QVOL*0.5*LENGTH
EGC=EWN+QWNA*(3.*OTHICK/16.+1./EMWN-0.5+3./(4.*OTHICK))
EGW=EWN+QWNA*(1./EMWN-0.5+3./(4.*OTHICK))
TGCA=(EGC/SIGMA)**0.25;TGWA=(EGW/SIGMA)**0.25
ACON=QWNA/(EGC-EWN)
    GROUP 3,4,5. X,Y,Z-direction grid specification
GRDPWR(Y,10,0.5*LENGTH,1.0);GRDPWR(Z,10,1.0,1.0)
    GROUP 6. Body-fitted coordinates or grid distortion
BFC=T;NONORT=T;INTEGER(NI,NJ,NK);NZ=10;NY=10
NK=NZ+1;NJ=NY+1;NI=2
SETPT(1,6,1,0.,0.15,0.);SETPT(1,6,4,0.,0.15,0.3)
SETPT(NI,6,1,1.,0.15,0.);SETPT(NI,6,4,1.,0.15,0.3)
DOMAIN(1,NI,6,6,1,4);SETLIN(YC,YF)
DOMAIN(1,NI,1,6,1,1);SETLIN(YC,LNJ*(YL-YF)+YF)
DOMAIN(1,NI,1,6,4,4);SETLIN(YC,LNJ*(YL-YF)+YF)
DOMAIN(1,NI,6,NJ,1,1);SETLIN(YC,LNJ*(YL-YF)+YF)
DOMAIN(1,NI,6,NJ,4,4);SETLIN(YC,LNJ*(YL-YF)+YF)
SETPT(1,6,8,0.,0.35,0.7);SETPT(1,6,NK,0.,0.35,1.0)
SETPT(NI,6,8,1.,0.25,0.7);SETPT(NI,6,NK,1.,0.35,1.0)
DOMAIN(1,NI,6,6,8,NK);SETLIN(YC,YF)
DOMAIN(1,NI,6,6,4,8);SETLIN(YC,LNK*(YL-YF)+YF)
DOMAIN(1,NI,1,6,NK,NK);SETLIN(YC,LNJ*(YL-YF)+YF)
DOMAIN(1,NI,6,NJ,NK,NK);SETLIN(YC,LNJ*(YL-YF)+YF)
DOMAIN(1,NI,1,6,1,NK);MAGIC(T);DOMAIN(1,NI,6,NJ,1,NK);MAGIC(T)
    GROUP 7. Variables stored, solved & named
CP1=1.0
MESG( Enter required energy variable ? (TEM1 or H1)
READVDU(CH1,CHAR,H1)
IF(:CH1:.EQ.TEM1) THEN
+ MESG( TEM1 solution selected
ELSE
+ MESG( H1 solution selected
+ TMP1=LINH;TMP1B=1.0/CP1
ENDIF
RADIAT(RADI,EMISS,SCATT,:CH1:);STORE(EMPO)
    GROUP 8. Terms (in differential equations) & devices
TERMS(:CH1:,N,N,N,N,P,P)
    GROUP 11. Initialization of variable or porosity fields
FIINIT(:CH1:)=TWN;FIINIT(SRAD)=EWN
  ** analytical solution
STORE(HA);ACON=3.0*KRAD*QVOL/8.0
DO JJ=1,NY
+PATCH(IN:JJ:,INIVAL,1,NX,JJ,JJ,1,NZ,1,1)
+GY=0.5*YFRAC(JJ)
IF(JJ.NE.1) THEN
+JJM1=JJ-1;GY=YFRAC(JJM1)+0.5*(YFRAC(JJ)-YFRAC(JJM1))
ENDIF
+GY=GY*YVLAST;EG=EGC-ACON*GY*GY
+TA=(EG/SIGMA)**0.25;INIT(IN:JJ:,HA,ZERO,TA)
ENDDO
    GROUP 13. Boundary conditions and special sources
   ** Net radiation flux from wall
PATCH(WALLRB,NORTH,1,NX,NY,NY,1,NZ,1,1)
COVAL(WALLRB,SRAD,2.*EMWN/(2.0-EMWN),EWN)
   ** uniformly-distributed volumetric heat source
PATCH(QHEAT,VOLUME,1,NX,1,NY,1,NZ,1,1)
COVAL(QHEAT,:CH1:,FIXFLU,QVOL)
    GROUP 15. Termination of sweeps
LSWEEP=300
    GROUP 16. Termination of iterations
SELREF=F;RESREF(:CH1:)=1.E-4*QRAD;RESREF(SRAD)=RESREF(:CH1:)
    GROUP 17. Under-relaxation devices
RELAX(:CH1:,LINRLX,0.3);RELAX(SRAD,LINRLX,0.6)
    GROUP 22. Spot-value print-out
IYMON=9;IZMON=7;NPLT=30;NYPRIN=1;ITABL=3
    GROUP 23. Field print-out and plot control
OUTPUT(:CH1:,Y,N,N,Y,Y,Y);PATCH(YWISE,PROFIL,1,1,1,NY,1,1,1,1)
PLOT(YWISE,:CH1:,0.0,0.0);PLOT(YWISE,SRAD,0.0,0.0)
    GROUP 24. Dumps for restarts
TSTSWP=-1