#$r002
   PHOTON USE
   AUTOPLOT
   FILE
   PHI 5
 
 
   d 1 h1 z 3;d 1 ha z 3;col3 1;blb4 2
   scale y 1000 1150;scale x 0 0.5;redr
   msg    temperature profile; press  to continue
   ENDUSE
TEXT(BFC RADIATION+HEAT SOURCE IN A TUBE:124
TITLE
  DISPLAY
   The problem considered is same as case 123 except that it is
   performed on an orthogonal BFC mesh with porosities present.
  ENDDIS
REAL(EMWN,TWN,EWN,TA,GY,EG,EGC,EGW,ACON)
REAL(QRAD,QVOL,KRAD,OTHICK,DIAM,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)
DIAM=1.0;KRAD=OTHICK/DIAM
SCATT=0.0;EMISS=KRAD-SCATT
EMWN=1.0;TWN=1000.0;EWN=SIGMA*TWN**4;QRAD=EWN;QVOL=EWN/DIAM
  ** analytical solution
QWNA=QVOL*0.25*DIAM
EGC=EWN+QWNA*(3.*OTHICK/16.+1./EMWN-0.5+9./(8.*OTHICK))
EGW=EWN+QWNA*(1./EMWN-0.5+9./(8.*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
CARTES=F;XULAST=0.1
NREGY=2;INTEGER(NYG,NYB);NYG=50;NYB=5
IREGY=1;GRDPWR(Y,NYG,0.5*DIAM,1.0)
IREGY=2;GRDPWR(Y,NYB,0.5*DIAM,1.0)
GRDPWR(Z,5,1.0,1.0)
    GROUP 6. Body-fitted coordinates or grid distortion
BFC=T;NONORT=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
   ** Deactive conduction & any built-in sources
TERMS(:CH1:,N,N,N,N,P,P)
    GROUP 11. Initialization of variable or porosity fields
FIINIT(:CH1:)=TWN;FIINIT(SRAD)=EWN
CONPOR(BLOCKV,0.0,CELL,1,NX,NYG+1,NY,1,NZ)
  ** extra CONPOR now required to zero south face
     of blockage, otherwise earth blows up
CONPOR(BLOCKN,0.0,NORTH,1,NX,NYG,NYG,1,NZ)
  ** analytical solution
STORE(HA);ACON=0.75*KRAD*QWNA/DIAM
DO JJ=1,NYG
+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,NYG,NYG,1,NZ,1,1)
COVAL(WALLRB,SRAD,2.*EMWN/(2.0-EMWN),EWN)
   ** uniformly-distributed volumetric heat source
PATCH(QHEAT,VOLUME,1,NX,1,NYG,1,NZ,1,1)
COVAL(QHEAT,:CH1:,FIXFLU,QVOL)
    GROUP 15. Termination of sweeps
LSWEEP=100
    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 18. Limits on variables or increments to them
VARMIN(:CH1:)=0.5*TWN
    GROUP 19. Data communicated by satellite to GROUND
    GROUP 20. Preliminary print-out
    GROUP 21. Print-out of variables
    GROUP 22. Spot-value print-out
IYMON=1;NPLT=10;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