c
C.... FILE NAME GXSUN.FTN-------------------------------- 070322
C This subroutine uses the OSGCANSEE utility to determine which cells are lit
C by the sun and applies the resulting heat sources..
C
C At the start of the Earth run there is a loop over all cells searching for
C cells with a solid surface or plate next to them, or which contain a cut cell.
C These are put into a list which is sent to the OSGCANSEE program. This program
C uses OpenSceneGraph to determine, given the direction of the sun and the
C geometry (read from facetdat), which cells in the list can see the sun. A new
C list is created with the illuminated cells marked. This is read back into Earth
C and based on the list heat sources are assigned to the cells.
C
C The file sent to OSGCANSEE is sunpos.txt, and the file sent back is sunpos.lit.
C
C Osgcansee.exe operates in the background with no visible signs. For a large case
C there is just a longer pause as Earth starts up.
C
C The optional STOREd variable SRF contains a 1 for all cells identified as either
C containing a surface or being adjacent to a blockage or plate. They are candidates
C for being illuminated by the sun.
C
C The optional STOREd variable LIT contains a 1 for surface cells which see the sun,
C and 0 for all others. A surface contour on all blockages with averaging off shows
C where the sun is shining.
C
C The optional STOREd variable #QS1 stores the actual TEM1 heat source for each cell,
C and #QS2 the heat source per unit area.
C
C When IMMERSOL is off, the heat is added to TEM1 in the participating solid cells, or
C in fluid cells next to the solid for 198 solids.
C When IMMERSOL is on, it is added to T3 in the solid cells or TEM1 in fluid cells.
C When IMMERSOL is on #QS3 stores the T3 heat source.
C
C SRF, LIT, #QS1/2/3 do not have to be stored, they are primarily debug to check that all
C is well.
C
C The sun direction is dot-producted with the normal to the surface to get the actual
C direct radiation heat flux in each coordinate direction. The diffuse radiation flux
C is added to all surface cells whether illuminated or not. This procedure is done
C once for steady state, and once per time step for transient. It happens in the
C subroutine SET_SUN_SOR.
C
C The direct radiation flux is multiplied by an absorption factor (defaulted to 1.0)
C before use. The absorption factor can be set individually for each BLOCKAGE, PLATE
C and THINPLT object. The optional store #SOL contains the absoption factors.
C
C BLOCKAGE, PLATE and THINPLT object types are assumed to be opaque, and INLET and OUTLET
C objects are assumed to be transparent. These settings can be reversed for each object.
C
C In transient cases, the sun position is calculated at the start of each time step.
C Only the first line of sunposlit is overwritten with the new sun position so the
C search for candidate cells is only done once.
C
C A single whole-domain, all-timestep PATCH named SUN with type CELL applies the
C heat sources via a COVAL(SUN, TEM1, FIXFLU, GRND1). If IMMERSOL is active, the
C heat sources for IMMERSOL are set via COVAL(SUN, T3, FIXFLU, GRND1). The heat sources
C are taken from the 3D store created earlier and are set into VAL.
C
C The following SPEDAT commands are used to pass data into the subroutine:
C
C SPEDAT(SET,SUNLIGHT,TIME,C,TIME_OF_DAY)
C where TIME_OF_DAY is the data and time in the format
C DAY/Month/day_of_month/hour.minute.second/year for example
C DAY/JUL/10/09.00.00/2011
C
C SPEDAT(SET,SUNLIGHT,SUNNORTH,R,ANGLE_TO_NORTH)
C where ANGLE_TO_NORTH is the angle between the grid y axis and North in degrees
C
C SPEDAT(SET,SUNLIGHT,SUNLATIT,R,LATITUDE)
C where LATITUDE is the local latitude in degrees (+ve North, -ve South)
C
C SPEDAT(SET,SUNLIGHT,DIRECT,C,FLAG)
C where FLAG can be CONSTANT for constant direct radiation or
C ALTITUDE for direct radiation as a function of solar altitude
C SPEDAT(SET,SUNLIGHT,DIRHEAT,R,QDIR)
C where QDIR is the constant direct radiation in W/m^2 when FLAG is CONSTANT
C or the DIRECT spedat is missing
C
C SPEDAT(SET,SUNLIGHT,DIFFUSE,C,FLAG)
C where FLAG can be CONSTANT for constant diffuse radiation or
C ALTITUDE for diffuse radiation as afunction of solar altitude
C SPEDAT(SET,SUNLIGHT,DIFHEAT,R,QDIF)
C where QDIF is the constant diffuse radiation in W/m^2 when FLAG is CONSTANT
C or the DIFFUSE spedat is missing
C
C SPEDAT(SET,SUNLIGHT,SKY,C,FLAG)
C where FLAG can be CLEAR for a clear sky or
C CLOUDY for a cloudy sky.
C The SKY flag is only used when DIFFUSE is set to ALTITUDE.
C
C The 'ALTITUDE' settings obtain the direct and diffuse solar radiation from
C polynomial fits to tables A2.24 and A2.25 of The CIBSE Guide, Volume A Design Data.
C These give the following values in W/m^2:
C !---------------!------------------!--------!----------!
C !Solar Altitude ! Direct Radiation ! Diffuse Radiation !
C ! ! Normal to Sun ! Clear ! Cloudy !
C !---------------!------------------!--------!----------!
C ! 5 ! 210 ! 25 ! 25 !
C ! 10 ! 390 ! 40 ! 50 !
C ! 20 ! 620 ! 65 ! 100 !
C ! 25 ! 690 ! 70 ! 125 !
C ! 30 ! 740 ! 75 ! 150 !
C ! 35 ! 780 ! 80 ! 175 !
C ! 40 ! 815 ! 85 ! 200 !
C ! 45 ! 840 ! 90 ! 225 !
C ! 50 ! 860 ! 95 ! 250 !
C ! 60 ! 895 ! 100 ! 300 !
C ! 70 ! 910 ! 105 ! 355 !
C ! 80 ! 920 ! 110 ! 405 !
C ! 90 ! 930 ! 115 ! 455 !
C !---------------!------------------!--------!----------!
C The solar altitude is obtained from the sun direction vector as calculated
C either at the start of the run, or at the start of each step in a transient.
C
C These tables are coded in GET_SUN_HEAT.
C
C If a weather data file has been specified in VR-Editor, the following settings
C are made.
C
C SPEDAT(SET,SUNLIGHT,WEATHERFILE,C,weatherfile.epw)
C where weatherfile.epw is the name of the file.
C If the calculation is steady-state the direct and diffuse solar radiation
C values obtained from the data file are passed as described above
C SPEDAT(SET,SUNLIGHT,DIRHEAT,R,QDIR)
C SPEDAT(SET,SUNLIGHT,DIFHEAT,R,QDIF)
C
C If the calculation is transient, a table of solar flux values is sent.
C SPEDAT(SET,SUNLIGHT,NHEATS,I,nheats)
C sets the number of lines in the table
C SPEDAT(SET,SUNLIGHT,NPHOUR,I,nphour)
C sets the number of entries per hour.
C SPEDAT(SET,SUNLIGHT,HEATn,C,qdir(n),qdif(n))
C sets the direct and diffuse heat loads for line n of the table
C The first line contains the values at elapsed time t=0, the start
C of time step 1.
C At the start of each time step the heat values at the mid-time of
C the step are obtained by linear interpolation in the table.
C
C SPEDAT(SET,object_name,SOLA,R,solar_absorption_factor)
C sets the absorption factor for a participating blockage
C SPEDAT(SET,patchname,SOLA,R,solar_absorption_factor)
C sets the absorption facotr for plate and thinplt objects - the patch
C name is used to allow different values for each side of the plate.
C
C SPEDAT(SET,object_name,TRANSPARENCY,R,transparency_value)
C sets the transparency of the named object. If not present, the object is
C assumed to be fully opaque - transparency=0. In this (18.98.15) implementation
C any non-zero transparency is taken to mean fully transparent.
C------------------------------------------------------------------------------------
SUBROUTINE GXSUN
USE weather_data
USE cut_cell_data
INCLUDE 'farray'
INCLUDE 'objnam'
INCLUDE 'patnos'
INCLUDE 'patcmn'
INCLUDE 'd_earth/parvar'
INCLUDE 'cmndmn'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'd_earth/pbcstr'
INCLUDE 'parear'
! INCLUDE 'satcfile'
COMMON /WORDI1/ NWDS,NCHARS(20),NSEMI,H,NLINES
INTEGER H
COMMON /WORDL1/ ERROR,MORWDS
COMMON /WORDC1/ WD(20),INLINE
LOGICAL ERROR,MORWDS
CHARACTER WD*20, INLINE*120
COMMON/F0/KXDX,KXXG,KXXU,IF01(2),KYDY,KYYG,IF02(6),KYYV,IF03(2),
1 KZDZ,KZZG,IF04(6),KZZW,IF05(280)
COMMON/GENI/NXNY,IGFIL1(8),NFM,IGF(21),IPRL,IBTAU,ILTLS,IGFIL(15),
1 ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,IRADZ,IVFOL
COMMON/DRHODP/ITEMP,IDEN/DVMOD/IDVCGR
COMMON /INTDMN/IDMN,NUMDMN,NXMAX,NYMAX,NZMAX,NXYMAX,NXYZMX,NCLTOT,
1 NDOMMX
COMMON/LRNTM3/L0UTAU,NMWALL,L0WALL,L0DSKN,IVPRST
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX USER SECTION STARTS:
COMMON/IEGWF3/ LBSOR,L0SOR0,L0T3,L0TWAL,LBSOR2,L0SOR2
C
REAL NORML(3),VELIN(3),VELOUT(3),DIR(3),SUNDIR(3)
PARAMETER (MAXDMN = 10)
CHARACTER*12 COBNAM,COBNM2,BLOCKNAM,PATNM*8,NAMPT*12
CHARACTER*68 CTIME,CVAL,LINE*256,COBTP*12,CSAV(2)*20
COMMON /FACES/L0FACE,L0FACZ
COMMON /VRMCMN/L0MSKZ
COMMON/NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
LOGICAL SLD,QEQ,QNE,FLUID1,DOIT,GRN,LINVRO,MASKPT,QLE,LPTDMN,LEZ,
1 LPAR,NF,SF,EF,WF,HF,LF,NWP,SWP,EWP,WWP,HWP,LWP,DONE,QGT,
1 VAC,POR,LSOLID,LCUTCL,L2DVRO,ALLCELLS,MASKED
LOGICAL LVRCEL
SAVE LPAR,DONE,ALLCELLS
SAVE L0LIT0,LBLIT,NUMSUN,ITM3,IZSUNL,LBSOR3,
1 L0SR30,L0WALSOLA,L0BLKSOLA,L0TRANS,L0TRANSW,L0TRANSF
SAVE SUNDIR,QDIR,QDIF
SAVE IQDIR,IQDIF,ISKY,ITIM,IDAY,ITIM0,NDAY0,INILIT
SAVE CSAV
DATA DONE /.FALSE./
C
NAMSUB='GXSUN'
C***********************************************************************
IF(IGR==1.AND.ISC==2) THEN
IF(.NOT.NULLPR) THEN
CALL WRYT40('GXSUN of: 070322')
CALL WRITST
CALL WRIT40('GXSUN of: 070322')
ENDIF
C... set index for T3
ITM3=LBNAME('T3'); LBSOR3=0
IF(NUMDMN>MAXDMN) THEN
CALL WRITBL
CALL WRITST
CALL WRIT40('Please increase MAXDMN in GXSUN')
CALL WRIT2I('Current ',MAXDMN,'; Needed',NUMDMN)
CALL WRITST
CALL SET_ERR(557,'MAXDMN too small in GXSUN',1)
RETURN
ENDIF
C.... Loop over all FGV domains
100 CONTINUE
IF(.NOT.DONE) THEN
CALL FILE_DEFAULTS ! ensure prefix/config have been read
DONE=.TRUE.
ENDIF
IF(IDMN>1) CALL INDXDM
NCELL=0
LPAR=MIMD.AND.NPROC>1 ! parallel flag
IF(LPAR) CALL BCAST_FILE('facetdat') ! send entire FACETDAT file
C... Get SUN settings
CVAL=' ';CALL GETSDC('SUNLIGHT','TIME',CVAL) ! start time
IF(CVAL/=' ') THEN
CTIME=CVAL; LT=LENGZZ(CTIME)
DO IC=1,LT
IF(CTIME(IC:IC)=='/') CTIME(IC:IC)=' '
IF(CTIME(IC:IC)=='.') CTIME(IC:IC)=':'
ENDDO
GNORTH=0.0; CALL GETSDR('SUNLIGHT','SUNNORTH',GNORTH) ! North direction
GLATIT=51.; CALL GETSDR('SUNLIGHT','SUNLATIT',GLATIT) ! Latitude
EPWNAM=' ';CALL GETSDC('SUNLIGHT','WEATHERFILE',EPWNAM) ! Weather file
NHEATS=0
IF(EPWNAM/=' '.AND..NOT.STEADY) THEN
CALL GETSDI('SUNLIGHT','NHEATS',NHEATS)
NPHOUR=1; CALL GETSDI('SUNLIGHT','NPHOUR',NPHOUR)
IF(NHEATS>0) THEN
ALLOCATE(QDR(NHEATS),QDF(NHEATS),STAT=IERR)
DO I=1,NHEATS
WRITE(LINE,'(''HEAT'',I4)') I
LL=8; CALL REMSPC(LINE,LL)
CALL GETSDC('SUNLIGHT',LINE(1:LL),CVAL)
LL=LENGZZ(CVAL)
CALL SPLTZZ(CVAL,WD,NWDS,NCHARS,LL)
QDR(I)=RRDZZZ(1); QDF(I)=RRDZZZ(2)
ENDDO
IQDIR=2; IQDIF=2
ENDIF
ELSE
CVAL=' '; CALL GETSDC('SUNLIGHT','DIRECT',CVAL)
IF(CVAL==' '.OR.CVAL=='CONSTANT') THEN
IQDIR=0; QDIR=1000.;CALL GETSDR('SUNLIGHT','DIRHEAT',QDIR)
ELSEIF(CVAL=='ALTITUDE') THEN
IQDIR=1
ENDIF
CVAL=' '; CALL GETSDC('SUNLIGHT','DIFFUSE',CVAL)
IF(CVAL==' '.OR.CVAL=='CONSTANT') THEN
IQDIF=0; QDIF=100.; CALL GETSDR('SUNLIGHT','DIFHEAT',QDIF)
ELSEIF(CVAL=='ALTITUDE') THEN
IQDIF=1
CVAL='CLEAR';CALL GETSDC('SUNLIGHT','SKY',CVAL)
ISKY=ITWO(0,1,CVAL=='CLEAR')
ENDIF
ENDIF
CALL WRITST
CALL WRIT40(' Solar Heating data')
CALL WRIT40(' ------------------')
IF(EPWNAM/=' ') THEN
LL=LENGZZ(EPWNAM)
WRITE(LUPR1,'('' Using weather data file '',A)')
1 EPWNAM(1:LL)
LL=LENGZZ(CTIME)
WRITE(LUPR1,'('' Time of day'',A)') CTIME(4:LL)
IF(NHEATS==0) THEN
WRITE(LUPR1,
1 '('' Direct solar heating '',1PE13.6,'' W/m^2'')') QDIR
WRITE(LUPR1,
1 '('' Diffuse solar heating '',1PE13.6,'' W/m^2'')') QDIF
ELSE
WRITE(LUPR1,'('' Time Direct Diffuse'')')
DO I=1,NHEATS
WRITE(LUPR1,'(I4,'' '',1PE13.6,'' '',1PE13.6)') I-1,
1 QDR(I),QDF(I)
ENDDO
ENDIF
ELSE
WRITE(LUPR1,'('' Latitude '',F6.2,A)') GLATIT,CHAR(176)
LL=LENGZZ(CTIME)
WRITE(LUPR1,'('' Time of day'',A)') CTIME(4:LL)
IF(IQDIR==0) THEN
WRITE(LUPR1,
1 '('' Constant direct solar heating '',1PE13.6,'' W/m^2'')')QDIR
ELSE
WRITE(LUPR1,
1 '('' Direct solar heating from solar altitude'')')
ENDIF
IF(IQDIF==0) THEN
WRITE(LUPR1,
1 '('' Constant diffuse solar heating '',1PE13.6,'' W/m^2'')')QDIF
ELSE
WRITE(LUPR1,
1 '('' Diffuse solar heating from solar altitude'')')
IF(ISKY==0) THEN
WRITE(LUPR1,'('' Sky assumed clear'')')
ELSE
WRITE(LUPR1,'('' Sky assumed cloudy'')')
ENDIF
ENDIF
ENDIF
CALL WRITST
OPEN(111,FILE='sunpos.txt',STATUS='UNKNOWN')
CLOSE(111,STATUS='DELETE')
OPEN(111,FILE='sunpos.txt',STATUS='NEW',FORM='FORMATTED',
1 ACCESS='DIRECT',CARRIAGECONTROL='LIST',RECL=80)
WRITE(111,REC=1,FMT='(A)') CTIME
WRITE(111,REC=2,FMT='(''Latitude '',1PE13.6)') GLATIT
WRITE(111,REC=3,FMT='(''North '',1PE13.6)') GNORTH
ELSE
RETURN ! time not set so do nothing
ENDIF
C
C... Loop over cells looking for potentially-lit surfaces. Write coordinates
C... to SUNPOS.TXT for OSGCANSEE to check
LBSRF=LBNAME('SRF'); L0SRF=0
LBLIT=LBNAME('LIT'); L0LIT0=0 ! 3D store for illuminated cells
IF(LBLIT>0.AND.STEADY) THEN !if LIT store exists
IWRIT=ITWO(-1,0,QEQ(FIINIT(LBLIT),READFI)) ! only write 1 line if RESTRT
ELSE ! LIT not present or transient, always write all lines of sunpos.txt
IWRIT=0 ! always write all lines of sunpos.txt
ENDIF
L0FACZ0=L0FACZ; IZ=1; NUMSUN=0; IREC=4
IZSUNL=-(NZ+1) ! this will be largest IZ with a surface
L0CELL=L0F(EASP1) ! store for cells to look at
ALLCELLS=.FALSE.; CALL GETSDL('SUNLIGHT','ALLCELLS',ALLCELLS)
DO IZZ=1,NZ ! loop over whole domain looking for surfaces
IZZM1=(IZZ-1)*NXNY; L0FACZ= L0FACE+IZZM1;ncell=0
IF(LBSRF>0) L0SRF=L0F(ANYZ(LBSRF,IZZ)) ! store for surfaces (debug)
DO IXX=1,NX; DO IYY=1,NY
IJ=(IXX-1)*NY+IYY; IJK=(IXX-1)*NY+IYY+IZZM1; IBLK=0
IF(.NOT.LSOLID(IJK))THEN ! fluid cell
IF(ISG62==0) THEN
IPBC = IPBSEQ(IDMN,IJK) ! get cut-cell number
ELSE
IF(ASSOCIATED(CUTCELL%ICUTH)) THEN
IPBC = CUTCELL%ICUTH(IJK)
ELSE
IPBC=0
ENDIF
ENDIF
IF(IPBC>0) THEN ! cell is cut, so must have some surface
IBLK=IBLK+1
ELSE ! check all 6 neighbours
IF(NF(IJ)) IBLK=IBLK+1; IF(SF(IJ)) IBLK=IBLK+1
IF(EF(IJ)) IBLK=IBLK+1; IF(WF(IJ)) IBLK=IBLK+1
IF(HF(IJ)) IBLK=IBLK+1; IF(LF(IJ)) IBLK=IBLK+1
IF(NWP(IJ)) IBLK=IBLK+1; IF(SWP(IJ)) IBLK=IBLK+1
IF(EWP(IJ)) IBLK=IBLK+1; IF(WWP(IJ)) IBLK=IBLK+1
IF(HWP(IJ)) IBLK=IBLK+1; IF(LWP(IJ)) IBLK=IBLK+1
ENDIF
L0PAT=L0PATNO(IDMN)+(IZZ-1)*NXNY ! index for patch-number store
CALL GET_SURFACE(IXX,IYY,IZZ,0,L0PAT,CAREA,NORML,IPLUS)
IF(CAREA>0.0) IBLK=IBLK+1
F(L0CELL+IJ)=IBLK; NCELL=NCELL+IBLK
ELSE ! solid cell, cannot be illuminated
F(L0CELL+IJ)=0
ENDIF
ENDDO; ENDDO
C... Search for patches from FOLIAGE object
DO IPT=1,NUMPAT
IF(NAMPAT(IPT)(1:4)=='TREE') THEN
MASKED=MASKPT(IPT)
CALL GETPAT(IPT,IR,TYP,JXF,JXL,JYF,JYL,JZF,JZL,JTF,JTL)
IF(IZZJZL) CYCLE ! not in range of patch
IF(MASKED) THEN
L0MSKZ=L0PVAR(22,IPT,IZZ) ! Set current address of VRMASK
IPW=0
ENDIF
DO IXX=JXF,JXL; DO IYY=JYF,JYL
IJ=(IXX-1)*NY+IYY
IF(MASKED) THEN
IPW=IPW+1
IF(LVRCEL(IPW)) THEN
F(L0CELL+IJ)=1.; NCELL=NCELL+1
ENDIF
ELSE
F(L0CELL+IJ)=1.; NCELL=NCELL+1
ENDIF
ENDDO; ENDDO
ENDIF
ENDDO
C... Having found the cells, may need to adjust the coordinates
IF(NCELL>0.OR.IZSUNL<0) THEN
DO IXX=1,NX; DO IYY=1,NY
IJ=(IXX-1)*NY+IYY
IF(SLD(IJ)) CYCLE
IBLK=F(L0CELL+IJ)
IJK=(IXX-1)*NY+IYY+IZZM1
IF(IBLK>0) THEN ! found some surface
IZSUNL=MAX(IZZ,IZSUNL)
NPBC=0; ICUT=0
IF(ISG62==0) THEN
IF(LCUTCL(IJK)) THEN
NPBC=NINT(F(NPBIJK(IDMN)+IJK))
ENDIF
ELSE
IF(ASSOCIATED(CUTCELL%ICUTH)) THEN
ICUT=CUTCELL%ICUTH(IJK)
ELSE
ICUT=0
ENDIF
ENDIF
C.... Offset for target location
EPSX=F(KXDX+IXX)*0.05
EPSY=F(KYDY+IYY)*0.05
EPSZ=F(KZDZ+IZZ)*0.05
IF(NPBC>0) THEN
IT1=ISHPB(IDMN,NPBC)
XPOS=F(IT1+PBC_SNYX) ! coords of sunny centre
YPOS=F(IT1+PBC_SNYY)
ZPOS=F(IT1+PBC_SNYZ)
C.... Move point away from surface in case it is still inside facets.
CALL VECTOR(NORML,F(PBC_CTVEC(1)+IT1),
1 F(PBC_CTVEC(2)+IT1),F(PBC_CTVEC(3)+IT1))
CALL NORM(NORML)
XPOS=XPOS+NORML(1)*EPSX
YPOS=YPOS+NORML(2)*EPSY
ZPOS=ZPOS+NORML(2)*EPSZ
ELSEIF(ICUT>0) THEN
XPOS=CUTCELL%CEN(1,ICUT) ! coords of sunny centre
YPOS=CUTCELL%CEN(2,ICUT)
ZPOS=CUTCELL%CEN(3,ICUT)
C.... Move point away from surface in case it is still inside facets.
CALL VECTOR(NORML,CUTCELL%NORML(1,ICUT),
1 CUTCELL%NORML(2,ICUT),CUTCELL%NORML(3,ICUT))
CALL NORM(NORML)
XPOS=XPOS+NORML(1)*EPSX
YPOS=YPOS+NORML(2)*EPSY
ZPOS=ZPOS+NORML(2)*EPSZ
ELSE ! un-cut cell
XPOS=F(KXXG+IXX) ! coords of cell centre
YPOS=F(KYYG+IYY)
ZPOS=F(KZZG+IZZ)
C... on face of PLATE, move test points to outer edge of cell, as in facetdat
C... plates are expanded to fill the cells on either side.
C... Move point to just outside cell face. This way originating plate
C... does not appear in list of hits.
IF(NWP(IJ)) YPOS=F(KYYV+IYY-1)-EPSY
IF(SWP(IJ)) YPOS=F(KYYV+IYY)+EPSY
IF(EWP(IJ)) XPOS=F(KXXU+IXX-1)-EPSX
IF(WWP(IJ)) XPOS=F(KXXU+IXX)+EPSX
IF(HWP(IJ)) ZPOS=F(KZZW+IZZ-1)-EPSZ
IF(LWP(IJ)) ZPOS=F(KZZW+IZZ)+EPSZ
ENDIF
IF(L0SRF>0) F(L0SRF+IJ)=1. ! save surface flag if stored
NUMSUN=NUMSUN+1 ! count number of surfaces
ELSEIF(ALLCELLS) THEN ! testing all cells
XPOS=F(KXXG+IXX) ! coords of cell centre
YPOS=F(KYYG+IYY)
ZPOS=F(KZZG+IZZ)
NUMSUN=NUMSUN+1 ! count number of 'surfaces'
IF(L0SRF>0) F(L0SRF+IJ)=0. ! save surface flag if stored
ELSE ! only looking for surface cells
IF(L0SRF>0) F(L0SRF+IJ)=0. ! save surface flag if stored
CYCLE
ENDIF
IF(IWRIT<=0) WRITE(111,REC=IREC,FMT='(I9,3(1PE13.6))')
1 IJK,XPOS,YPOS,ZPOS ! write to output file
IF(IWRIT<0) IWRIT=1
IREC=IREC+1
ENDDO; ENDDO
ENDIF
ENDDO
IZSUNL=MAX(0,MIN(IZSUNL+1,NZ))
CLOSE(111,STATUS='KEEP')
IF(LBLIT==0) THEN
CALL GXMAKE0(L0LIT0,NX*NY*NZ,'SUNLIT')
ENDIF
LBSOR=LBNAME('#QS1'); L0SOR0=0
IF(LBSOR==0) THEN
CALL GXMAKE0(L0SOR0,NX*NY*IZSUNL,'QS1')
ENDIF
LBSOR2=LBNAME('#QS2'); L0SOR2=0
IF(LBSOR2==0) THEN
CALL GXMAKE0(L0SOR2,NX*NY*IZSUNL,'QS2')
ENDIF
IF(ITM3>0) THEN
LBSOR3=LBNAME('#QS3'); L0SR30=0
IF(LBSOR3==0) THEN
CALL GXMAKE0(L0SR30,NX*NY*IZSUNL,'SUNSOR3')
ENDIF
ENDIF
C... must use global patch number in parallel
NUMPT=ITWO(GD_NUMPAT,NUMPAT,LPAR)
C... get solar absorbtion factors for plates, and store them
L0WALLSOLA=0; L0TRANSW=0
WALLS: DO IP=1,NUMPT
IF(LPAR) THEN ! get global object name
COBNAM=GD_OBJNAM(IP); IF(COBNAM==' ') CYCLE
NAMPT=GD_NAMPAT(IP)(2:)
ELSE ! get local object name
COBNAM=OBJNAM(IP); IF(COBNAM=='NOTSET') CYCLE
NAMPT=NAMPAT(IP)
ENDIF
COBTP=' '; CALL GETSDC('OBJTYP','^'//NAMPT,COBTP) ! get object type
IF(COBTP=='PLATE') THEN ! process PLATE objects
IF(LPAR) THEN
IPAT=GD_INDPAT(IP,1) ! get local patch number. -ve means not on this proc.
IF(IPAT<0) CYCLE
ELSE
IPAT=IP
ENDIF
C... only create storage for solar absorption factor and
C transparency factor if non-default values found.
IF(.NOT.LPTDMN(IPAT)) CYCLE
C... get parent patch type and limits (IPAT is local patch)
CALL GETPAT(IPAT,IREG,TYP,IX1,IX2,IY1,IY2,IZ1,IZ2,IT1,IT2)
SOLA=1.0; CALL GETSDR(NAMPAT(IPAT),'SOLA',SOLA)
IF(QNE(SOLA,1.0)) THEN ! factor /= 1.0
IF(L0WALSOLA==0) THEN ! not stored yet
CALL GXMAKE(L0WALSOLA,NX*NY*NZ,'WALLSOL') ! allocate storage
DO II=1,NX*NY*NZ; F(L0WALSOLA+II)=1.0; ENDDO ! initialise to 1.0
ENDIF
DO IZZ=IZ1,IZ2
IZADD=(IZZ-1)*NXNY
IF(MASKPT(IPAT))THEN
L0MSKZ=L0PVAR(22,IPAT,IZZ); IPW=0
DO IXX=IX1,IX2; DO IYY=IY1,IY2
IPW=IPW+1; IJK=(IXX-1)*NY+IYY+IZADD
IF(L2DVRO(IPW)) THEN
F(L0WALSOLA+IJK)=SOLA
ENDIF
ENDDO; ENDDO
ELSE
DO IXX=IX1,IX2; DO IYY=IY1,IY2
IJK=(IXX-1)*NY+IYY+IZADD
F(L0WALSOLA+IJK)=SOLA
ENDDO; ENDDO
ENDIF
ENDDO
ENDIF
TRANS=0.0; CALL GETSDR(OBJNAM(IPAT),'TRANSPARENCY',TRANS) ! get transparency factor
IF(QNE(TRANS,0.0)) THEN ! /= 0.0 means fully transparent
IF(L0TRANSW==0) THEN ! not stored yet
CALL GXMAKE0(L0TRANSW,NX*NY*NZ,'WALLTRANSP') ! allocate storage
ENDIF
DO IZZ=IZ1,IZ2
IZADD=(IZZ-1)*NXNY
IF(MASKPT(IPAT))THEN
L0MSKZ=L0PVAR(22,IPAT,IZZ); IPW=0
DO IXX=IX1,IX2; DO IYY=IY1,IY2
IPW=IPW+1; IJK=(IXX-1)*NY+IYY+IZADD
IF(L2DVRO(IPW)) THEN
F(L0TRANSW+IJK)=TRANS
ENDIF
ENDDO; ENDDO
ELSE
DO IXX=IX1,IX2; DO IYY=IY1,IY2
IJK=(IXX-1)*NY+IYY+IZADD
F(L0TRANSW+IJK)=TRANS
ENDDO; ENDDO
ENDIF
ENDDO
ENDIF
ENDIF
ENDDO walls
C... only allocate L0BLKSOLA if any factor/=1. Allow for transparency flag
L0BLKSOLA=0; CALL GXMAKE0(L0TRANS,NUMPT,'TRANSPARENCY') ! allocate storage
L0TRANSF=0
DO IP=1,NUMPT
COBTP=' '
IF(LPAR) THEN ! get global object and patch names for parallel
IPAT=GD_INDPAT(IP,1) ! get local patch number. -ve means not on this proc.
IF(IPAT<0) CYCLE
COBNAM=GD_OBJNAM(IP); PATNM=GD_NAMPAT(IP)
ELSE
IF(.NOT.LPTDMN(IP)) CYCLE ! not this domain
COBNAM=OBJNAM(IP)
IF(MASKPT(IP)) THEN
PATNM='^'//NAMPAT(IP)
ELSE
PATNM='!'//NAMPAT(IP)
ENDIF
ENDIF
CALL GETSDC('OBJTYP',PATNM,COBTP) ! get object type
IF(COBTP=='BLOCKAGE') THEN ! if it is a blockage
SOLA=1.0; CALL GETSDR(COBNAM,'SOLA',SOLA) ! get solar absorption factor
IF(QNE(SOLA,1.0)) THEN ! factor /= 1
IF(L0BLKSOLA==0) THEN ! not stored yet
CALL GXMAKE(L0BLKSOLA,NUMPT,'BLKSOL') ! allocate storge
DO II=1,NUMPT; F(L0BLKSOLA+II)=1.0; ENDDO ! initialise to 1.0
ENDIF
F(L0BLKSOLA+IP)=SOLA ! and save it
ENDIF
ENDIF
IF(COBTP=='INLET'.OR.COBTP=='OUTLET'.OR.COBTP=='OPENING'.OR.
1 COBTP=='APERTURE'.OR.COBTP=='FAN') THEN
TRAN0=1.0
ELSE
TRAN0=0.0
ENDIF
TRANS=TRAN0; CALL GETSDR(COBNAM,'TRANSPARENCY',TRANS) ! get transparency factor
F(L0TRANS+IP)=TRANS ! save value for this object (save global patch number for parallel)
ENDDO
INILIT=ITWO(1,0,ALLCELLS)
IF(.NOT.STEADY) RETURN ! for transient get lit cells at start of step
C... Invoke OSGCANSEE to get list of illuminated points
CALL GET_SUN_LIT(NUMSUN,LBLIT,L0LIT,L0LIT0,SUNDIR,INILIT,
1 L0TRANS,L0TRANSF)
C... Get solar heating from solar altitude or weather file
CALL GET_SUN_HEAT(SUNDIR(3),IQDIR,IQDIF,ISKY,QDIR,QDIF,LUPR1)
C... Set solar heat sources for use in Group 13
CALL SET_SUN_SOR(LBLIT,L0LIT0,LBSOR,L0SOR0,LBSOR3,L0SR30,IZSUNL,
1 SUNDIR,QDIR,QDIF,L0WALSOLA,L0BLKSOLA,L0TRANS,L0TRANSW,
1 L0TRANSF)
C ... loop back for next FGV domain
! 101 CONTINUE
IF(IDMN1) THEN
IDMN=1
CALL INDXDM
ENDIF
C*****************************************************************
C
C--- GROUP 13. Boundary conditions and special sources
C Index for Coefficient - CO
C Index for Value - VAL
ELSEIF(IGR==13.AND.ISC==13) THEN
C------------------- SECTION 14 ------------------- value = GRND1
C
IF(NPATCH=='SUN') THEN
IF(INDVAR==ITEM1.OR.INDVAR==ITM3.OR.INDVAR==H1) THEN
IF(IZ>IZSUNL) THEN ! above highest surface, just set zero
CALL FN1(VAL,0.0)
ELSE ! copy source into VAL
L0VAL=L0F(VAL)
LBSR=ITWO(LBSOR,LBSOR3,INDVAR==ITEM1.OR.INDVAR==H1)
IF(LBSR>0) THEN
L0SOR=L0F(LBSR)
ELSE
L0SR0=ITWO(L0SOR0,L0SR30,INDVAR==ITEM1.OR.INDVAR==H1)
L0SOR=L0SR0+(IZ-1)*NXNY
ENDIF
DO I=1,NXNY
F(L0VAL+I)=F(L0SOR+I)
ENDDO
ENDIF
ENDIF
ENDIF
ELSEIF(IGR==19.AND.ISC==1) THEN
C * ------------------- SECTION 1 ---- Start of time step.
C... must set INDVAR to 0 here to ensure calls to GET_SURFACE do not
C enter the INDVAR==W1 section and return false area
INDVAR0=INDVAR; INDVAR=0
CVAL=' ';CALL GETSDC('SUNLIGHT','TIME',CVAL)
IF(CVAL/=' ') THEN
IF(ISTEP==FSTEP) THEN
CTIME=CVAL; LT=LENGZZ(CTIME)
DO IC=1,LT
IF(CTIME(IC:IC)=='/') CTIME(IC:IC)=' '
IF(CTIME(IC:IC)=='.') CTIME(IC:IC)=' '
ENDDO
CALL SPLTZZ(CTIME,WD,NWDS,NCHARS,LT)
CSAV(1)=WD(1); CSAV(2)=WD(2)
IDAY=IRDZZZ(3);IHOUR=IRDZZZ(4);IMINS=IRDZZZ(5)
ISECS=IRDZZZ(6)
ITIM0=IHOUR*60*60+IMINS*60+ISECS
NDAY0=0
ENDIF
C... ITIM is time in sec set in Q1. Add elapsed time of run to get current time of day.
ITIM=ITIM0+TIM-DT/2. ! time is at mid-point of time step to get average
IHOUR=ITIM/3600
ITIM1=ITIM-IHOUR*3600
IMINS=ITIM1/60
ISECS=ITIM1-IMINS*60
IF(IHOUR>24) THEN
NDAYS=IHOUR/24
IF(NDAY0/=NDAYS) THEN
NDAY0=NDAYS; IDAY=IDAY+1
ENDIF
IHOUR=IHOUR-NDAYS*24
ENDIF
WD(1)=CSAV(1); NCHARS(1)=LENGZZ(WD(1))
WD(2)=CSAV(2); NCHARS(2)=LENGZZ(WD(2))
WRITE(WD(3),'(I2.2)') IDAY
WRITE(WD(4),'(I2.2)') IHOUR
WRITE(WD(5),'(I2.2)') IMINS
WRITE(WD(6),'(I2.2)') ISECS
CTIME=WD(1)(1:NCHARS(1))//' '//WD(2)(1:NCHARS(2))//' '
1//WD(3)(1:2)//' '//WD(4)(1:2)//':'//WD(5)(1:2)//':'//WD(6)(1:2)
1//' '//WD(7)
C... overwrite first line of sunpos with updated time.
OPEN(111,FILE='sunpos.txt',STATUS='OLD',FORM='FORMATTED',
1 ACCESS='DIRECT',CARRIAGECONTROL='LIST',RECL=80)
WRITE(111,REC=1,FMT='(A)') CTIME(1:LENGZZ(CTIME))
CLOSE(111,STATUS='KEEP')
C... Invoke OSGCANSEE to get new list of illuminated points
CALL GET_SUN_LIT(NUMSUN,LBLIT,L0LIT,L0LIT0,SUNDIR,INILIT,
1 L0TRANS,L0TRANSF)
C... Get solar heating from solar altitude or weather file
CALL GET_SUN_HEAT(SUNDIR(3),IQDIR,IQDIF,ISKY,QDIR,QDIF,LUPR1)
C... Set solar heat sources for use in Group 13
CALL SET_SUN_SOR(LBLIT,L0LIT0,LBSOR,L0SOR0,LBSOR3,L0SR30,
1 IZSUNL,SUNDIR,QDIR,QDIF,L0WALSOLA,L0BLKSOLA,L0TRANS,L0TRANSW,
1 L0TRANSF)
ENDIF
INDVAR=INDVAR0 ! restore original INDVAR in case it is used later
ENDIF
C****************************************************************
END
C----------------------------------------------------------------------
SUBROUTINE SET_SUN_SOR(LBLIT,L0LIT0,LBSOR,L0SOR0,LBSOR3,L0SR30,
1 IZ2,SUNDIR,QDIR0,QDIF0,L0WALSOLA,L0BLKSOLA,L0TRANS,L0TRANSW,
1 L0TRANSF)
C... Subroutine to fill the array(s) containing the solar heat sources.
C It is called once in steady state, and once at the start of each time
C step in transient.
C If Immersol is active, heat sources are added to T3 on the inside of
C participating solids, and to TEM1 in cells adjacent to non-participating
C blockages, plates and thin plates.
C If Immersol is not active, all sources are added to TEM1.
C
C The optional STOREs #QS1,#QS2 and #QS3 contain the TEM1 source per cell,
C the total heat source per-unit-area and the T3 source per cell respectively.
C The optional STORE #SOL contains the solar absorption coefficient.
USE cut_cell_data
INCLUDE 'farray'
INCLUDE 'patnos'
INCLUDE 'cmndmn'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
INCLUDE 'd_earth/pbcstr'
COMMON /GEODMN0/ I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
COMMON/GENI/NXNY,IGFIL1(8),NFM,IGF(21),IPRL,IBTAU,ILTLS,IGFIL(15),
1 ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,IRADZ,IVFOL
COMMON /INTDMN/IDMN,NUMDMN,NXMAX,NYMAX,NZMAX,NXYMAX,NXYZMX,NCLTOT,
1 NDOMMX
COMMON /FACES/L0FACE,L0FACZ
1 /SODAI/IFILL1(6),KGEOM,IFILL2(8)
1 /F0/ KXDX,LB1(4),KYDY,KYYG,KYR,KYR2,KYRV,LB2(6),KZDZ,
1 LB3(283)
COMMON/IEGWF3/ IFIL(4),LBSOR2,L0SOR2
LOGICAL SLD,QEQ,QNE,FLUID1,DOIT,GRN,LINVRO,MASKPT,QLE,LPTDMN,LEZ,
1 LPAR,NF,SF,EF,WF,HF,LF,NWP,SWP,EWP,WWP,HWP,LWP,DONE,QGT,
1 VAC,POR,IST3,LCUTCL
REAL SUNDIR(3),NORML(3),NORM0(3)
INTEGER IMAXDIR(1)
C
IST3=LBSOR3>0.OR.L0SR30>0
! LBSOR=LBNAME('#QS1')
L0SOR=0; ILIM=ITWO(NZ,IZ2,LBSOR>0)
DO IZZ=1,ILIM ! initialise TEM1 source stores to zero
IF(LBSOR>0) THEN
L0SOR=L0F(ANYZ(LBSOR,IZZ))
ELSE
L0SOR=L0SOR0+(IZZ-1)*NXNY
ENDIF
DO IJ=1,NXNY
F(L0SOR+IJ)=0.0
ENDDO
ENDDO
! LBSOR3=LBNAME('#QS3')
IF(IST3) THEN
ILIM=ITWO(NZ,IZ2,LBSOR3>0)
DO IZZ=1,ILIM ! initialise T3 source stores to zero
IF(LBSOR3>0) THEN
L0SOR=L0F(ANYZ(LBSOR3,IZZ))
ELSE
L0SOR=L0SR30+(IZZ-1)*NXNY
ENDIF
DO IJ=1,NXNY
F(L0SOR+IJ)=0.0
ENDDO
ENDDO
ENDIF
ILIM=ITWO(NZ,IZ2,LBSOR2>0)
DO IZZ=1,ILIM ! initialise per unit area source stores to zero
IF(LBSOR2>0) THEN
L0SR2=L0F(ANYZ(LBSOR2,IZZ))
ELSE
L0SR2=L0SOR2+(IZZ-1)*NXNY
ENDIF
DO IJ=1,NXNY
F(L0SR2+IJ)=0.0
ENDDO
ENDDO
LBSOLA=LBNAME('#SOL'); L0SOLA=0
IF(LBSOLA>0) THEN
DO IZZ=1,NZ ! initialise solar absorbtion stores to zero
L0SOLA=L0F(ANYZ(LBSOLA,IZZ))
DO IJ=1,NXNY
F(L0SOLA+IJ)=0.0
ENDDO
ENDDO
ENDIF
L0FACZ0=L0FACZ !; L0SOR2=0
DO IZZ=1,IZ2 ! loop over slabs containing surface
IF(LBLIT/=0) THEN ! store containing lit flag
L0LIT=L0F(ANYZ(LBLIT,IZZ))
ELSE
L0LIT=L0LIT0+(IZZ-1)*NXNY
ENDIF
IF(LBSOR/=0) THEN ! store for TEM1 heat source
L0SOR1=L0F(ANYZ(LBSOR,IZZ)) ! in 3D store
ELSE
L0SOR1=L0SOR0+(IZZ-1)*NXNY ! in local 3D store
ENDIF
IF(IST3) THEN ! store for T3 heat source
IF(LBSOR3/=0) THEN
L0SOR3=L0F(ANYZ(LBSOR3,IZZ)) ! in 3D store
ELSE
L0SOR3=L0SR30+(IZZ-1)*NXNY ! in local 3D store
ENDIF
ELSE
L0SOR3=L0SOR1
ENDIF
!!! L0SOR=ITWO(L0SOR3,L0SOR1,IST3)
L0SOR=L0SOR1
IF(LBSOR2>0) THEN ! optional store for per-unit-area source
L0SR2=L0F(ANYZ(LBSOR2,IZZ))
ELSE
L0SR2=L0SOR2+(IZZ-1)*NXNY
ENDIF
IF(LBSOLA>0) THEN ! optional store for solar absorption factor
L0SOLA=L0F(ANYZ(LBSOLA,IZZ))
ENDIF
IZZM1=(IZZ-1)*NXNY; L0FACZ= L0FACE+IZZM1
IADDZ=ITWO(NXNY,NFM,LBSOR==0)
IADD2=ITWO(NXNY,NFM,LBSOR2==0)
IADD3=ITWO(NXNY,NFM,LBSOR3==0)
L0PAT=L0PATNO(IDMN) ! index for patch-number store
IF(NPOR>0) L0NPOR=L0F(ANYZ(NPOR,IZZ))
IF(EPOR>0) L0EPOR=L0F(ANYZ(EPOR,IZZ))
IF(HPOR>0) L0HPOR=L0F(ANYZ(HPOR,IZZ))
DO IX=IXF,IXL
DO IY=IYF,IYL
IJ=(IX-1)*NY+IY
IJK=IJ+IZZM1
IF(SLD(IJ)) CYCLE ! skip solid cells
IF(NINT(F(L0LIT+IJ))>0) THEN ! the cell is illuminated
QDIR1=QDIR0 ! use full direct heat
ELSE
QDIR1=0.0 ! cell not illuminated, set direct heat to zero
ENDIF
SOLA=1.0 ! initialise solar absorption factor
TRANS=0.0 ! initialise transparency factor
C
NPBC=0; ICUT=0; CAREA=0.0; IFACE=0
IF(ISG62==0) THEN ! Parsol
IF(LCUTCL(IJK)) THEN ! cell is CUT
NPBC= IPBSEQ(IDMN,IJK)
ENDIF
ELSE ! Sparsol
IF(ASSOCIATED(CUTCELL%ICUTH)) THEN
ICUT=CUTCELL%ICUTM(IJK)
ENDIF
ENDIF
C
IF(NPBC>0) THEN ! PARSOL
IT1 = ISHPB(IDMN,NPBC)
AREA=F(PBC_CTAR+IT1)
CALL VECTOR(NORML,F(PBC_CTVEC(1)+IT1),
1 F(PBC_CTVEC(2)+IT1),F(PBC_CTVEC(3)+IT1))
ELSEIF(ICUT>0) THEN ! SPARSOL
AREA=CUTCELL%AREA(ICUT)
CALL VECTOR(NORML,CUTCELL%NORML(1,ICUT),
1 CUTCELL%NORML(2,ICUT),CUTCELL%NORML(3,ICUT))
ENDIF
C
IF(NPBC>0.OR.ICUT>0) THEN ! Parsol or Sparsol cut cell
CALL NORM(NORML) ! unit-normal to cut surface
ISOR=IJ; ISOR2=IJ
IF(NPBC>0) THEN ! Parsol / Xparsol
IPAT=NINT(F(L0PAT+IJK))
IF(L0BLKSOLA>0.AND.IPAT>0) SOLA=F(L0BLKSOLA+IPAT)
IF(L0TRANS>0.AND.IPAT>0) TRANS=F(L0TRANS+IPAT)
TRANSF=1.0; IF(L0TRANSF>0) TRANSF=F(L0TRANSF+IJK)
SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
SOR=SOLA*AREA*(QDIR1*TRANSF*
1 AMAX1(DOT(SUNDIR,NORML),0.0)+QDIF0)
IF(L0SR2>0) THEN ! recover summed area
AREA0=F(L0SOR+ISOR)/(F(L0SR2+IJ)+TINY)
AREA=AREA+AREA0 ! sum current area
ENDIF
F(L0SOR+ISOR)=F(L0SOR+ISOR)+SOR
IF(L0SR2>0) F(L0SR2+IJ)=F(L0SOR+ISOR)/(AREA+TINY)
IF(LBSOLA>0) F(L0SOLA+IJ)=SOLA
IMAXDIR=MAXLOC(ABS(NORML)) ! biggest normal
IF(IMAXDIR(1)==1) THEN ! mainly X, E
IFACE=3
ELSEIF(IMAXDIR(1)==2) THEN ! mainly Y, N
IFACE=1
ELSE ! mainly Z, H
IFACE=5
ENDIF
IF(NORML(IMAXDIR(1))>0.0) THEN ! switch to W, S, L
IFACE=IFACE+1
ENDIF
CAREA=0.0
ELSE ! Sparsol
IMAXDIR=MAXLOC(ABS(NORML)) ! biggest normal
IF(IMAXDIR(1)==1) THEN ! mainly X, E
IAD1=NY; IAD2=IAD1; IAD3=IAD1; IADS=IAD1;IFACE=3
ELSEIF(IMAXDIR(1)==2) THEN ! mainly Y, N
IAD1=1; IAD2=IAD1; IAD3=IAD1; IADS=IAD1; IFACE=1
ELSE ! mainly Z, H
IAD1=IADDZ
IAD2=ITWO(NXNY,NFM,LBSOR2==0)
IAD3=NXNY; IADS=NFM; IFACE=5
ENDIF
IF(NORML(IMAXDIR(1))>0.0) THEN ! switch to W, S, L
IAD1=-IAD1; IAD2=-IAD2; IAD3=-IAD3; IADS=-IADS
IFACE=IFACE+1
ENDIF
DZ=F(KZDZ+IZZ); DX=F(KXDX+IX); DY=F(KYDY+IY) ! cell sizes
IF(.NOT.CARTES.AND.NX>1) DX=2.0*SIN(0.5*DX)
IF(IFACE==1) THEN ! N Recover true face area
AREA=DX*F(KYRV+IY)*DZ
ELSEIF(IFACE==2) THEN ! S
AREA=DX*F(KYRV+IY-1)*DZ
ELSEIF(IFACE<5) THEN ! E & W
AREA=DY*DZ
ELSE ! H & L
AREA=DX*DY*F(KYR+IY)
ENDIF
IPAT=NINT(F(L0PAT+IJK+IAD3)) ! need to add offset into neighbour
IF(L0BLKSOLA>0.AND.IPAT>0) SOLA=F(L0BLKSOLA+IPAT)
IF(L0TRANS>0.AND.IPAT>0) TRANS=F(L0TRANS+IPAT)
TRANSF=1.0; IF(L0TRANSF>0) TRANSF=F(L0TRANSF+IJK)
SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
SOR=SOLA*(AREA*(QDIR1*TRANSF*
1 AMAX1(DOT(SUNDIR,NORML),0.0)+QDIF0))
IF(SLD(IJ+IAD3)) THEN !
IF(POR(IJ+IAD3).OR.VAC(IJ+IAD3)) THEN ! if neighbour is 198/199 keep source in cur
ISOR=IJ; ISOR2=IJ; IAD1=0; IAD2=0; IADS=0
ELSE ! if neighbour is participating, set source in solid
ISOR=IJ+IAD1; ISOR2=IJ+IAD2
ENDIF
ENDIF
IF(SLD(IJ+IAD3)) THEN ! neighbour is solid of some kind
IF(L0SR2>0) THEN ! recover summed area
AREA0=F(L0SOR+ISOR)/(F(L0SR2+ISOR2)+TINY)
AREA=AREA+AREA0 ! sum current area
ENDIF
F(L0SOR+ISOR)=F(L0SOR+ISOR)+SOR
IF(L0SR2>0) F(L0SR2+ISOR2)=F(L0SOR+ISOR)/(AREA+TINY)
IF(LBSOLA>0) F(L0SOLA+IJ+IADS)=SOLA
ENDIF
CAREA=0.0
ENDIF
CALL COPYVC(NORML,NORM0) ! save cut-cell normal
ELSEIF(ISG62==1) THEN
CALL VECTOR(NORM0,0.,0.,0.)
L0PATIZZ=L0PATNO(IDMN)+(IZZ-1)*NXNY ! index for patch-number store
CALL GET_SURFACE(IX,IY,IZZ,0,L0PATIZZ,CAREA,NORM0,IPLUS)
ELSE
CALL VECTOR(NORM0,0.,0.,0.)
ENDIF
C
IF ((NF(IJ).OR.NWP(IJ).OR.CAREA>0.0).AND.IFACE/=1) THEN ! Blocked face on NORTH side
CALL VECTOR(NORML,0.,-1.,0.) ! set unit vector normal to surface
IF(DOT(NORML,NORM0)>0..OR.CAREA==0.) THEN ! blocked cell is not 'behind' cut cell
IF(DOT(NORML,NORM0)>0..AND.CAREA>0.0) CAREA=0.0
NEXT=ITWO(1,0,IY0) THEN ! neighbour is cut-cell
ICUTIDXN = ISHPB(IDMN,IPBCN)
AREA=F(PBC_DAKAR(2)+ICUTIDXN) ! surface area for cut-cell
ELSE ! neighbour is fully blocked
AREA=F(I3DANY+IJK) ! surface area from cell area
C... divide area by porosity if on face of porosity
IF(NWP(IJ).AND.NPOR>0) AREA=AREA/(F(L0NPOR+IJ)+TINY)
ENDIF
IPAT=NINT(F(L0PAT+IJKN)) ! patch number of blockage affecting this cell
SOLA=1.0; TRANS=0.0; TRANSF=1.0
IF(NF(IJ).AND.IPAT>0) THEN ! face is blocked by blockage
IF(L0BLKSOLA>0) SOLA=F(L0BLKSOLA+IPAT) ! get absorption factor for blockage
IF(L0TRANS>0) TRANS=F(L0TRANS+IPAT) ! get transparency factor for blockage
ELSE ! face is blocked by plate
IF(L0WALSOLA>0) SOLA=F(L0WALSOLA+IJK); IJN=IJ ! get absoption factor for plate
IF(L0TRANSW>0) TRANS=F(L0TRANSW+IJK) ! get transparency factor for plate
ENDIF
IF(L0TRANSF>0) TRANSF=F(L0TRANSF+IJK)
SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
QDIR=SOLA*QDIR1*TRANSF; QDIF=SOLA*QDIF0 ! multiply radiation by absorption factor
C... Set the source into either the neighbour cell or the current cell
CALL SETSOR(IJ,IJN,IJN,IJN,IJN,IST3,L0SOR3,L0SOR1,QDIR,
1 QDIF,AREA,SUNDIR,NORML,L0SR2,TINY)
IF(LBSOLA>0) F(L0SOLA+IJN)=SOLA ! save absorption factor for plotting
ENDIF
ENDIF
C
IF((SF(IJ).OR.SWP(IJ).OR.CAREA>0.0).AND.IFACE/=2) THEN ! Blocked face on SOUTH side
CALL VECTOR(NORML,0.,1.,0.) ! set unit vector normal to surface
IF(DOT(NORML,NORM0)>0..OR.CAREA==0.) THEN ! blocked cell is not 'behind' cut cell
IF(DOT(NORML,NORM0)>0..AND.CAREA>0.0) CAREA=0.0
NEXT=ITWO(-1,0,IY>1)
IJKS=IJK+NEXT ! neighbour index
IPBCS = IPBSEQ(IDMN,IJKS) ! neighbour cut-cell index
IJS=IJ+NEXT
IF(PARSOL.AND.IPBCS>0) THEN ! neighbour is cut-cell
ICUTIDXS = ISHPB(IDMN,IPBCS)
AREA=F(PBC_DAKAR(1)+ICUTIDXS) ! surface area for cut-cell
ELSE ! neighbour is fully blocked
AREA=F(I3DANY+IJKS) ! surface area from cell area
IF(SWP(IJ).AND.NPOR>0) AREA=AREA/(F(L0NPOR+IJS)+TINY)
ENDIF
IPAT=NINT(F(L0PAT+IJKS))
SOLA=1.0; TRANS=0.0; TRANSF=1.0
IF(SF(IJ).AND.IPAT>0) THEN
IF(L0BLKSOLA>0) SOLA=F(L0BLKSOLA+IPAT)
IF(L0TRANS>0) TRANS=F(L0TRANS+IPAT) ! get transparency factor for blockage
ELSE
IF(L0WALSOLA>0) SOLA=F(L0WALSOLA+IJK); IJS=IJ
IF(L0TRANSW>0) TRANS=F(L0TRANSW+IJK) ! get transparency factor for plate
ENDIF
IF(L0TRANSF>0) TRANSF=F(L0TRANSF+IJK)
SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
QDIR=SOLA*QDIR1*TRANSF; QDIF=SOLA*QDIF0
CALL SETSOR(IJ,IJS,IJS,IJS,IJS,IST3,L0SOR3,L0SOR1,QDIR,
1 QDIF,AREA,SUNDIR,NORML,L0SR2,TINY)
IF(LBSOLA>0) F(L0SOLA+IJS)=SOLA
ENDIF
ENDIF
C
IF((EF(IJ).OR.EWP(IJ).OR.CAREA>0.0).AND.IFACE/=3) THEN ! Blocked face on EAST side
CALL VECTOR(NORML,-1.,0.,0.) ! set unit vector normal to surface
IF(DOT(NORML,NORM0)>0..OR.CAREA==0.) THEN ! blocked cell is not 'behind' cut cell
IF(DOT(NORML,NORM0)>0..AND.CAREA>0.0) CAREA=0.0
NEXT=ITWO(NY,0,IX0) THEN ! neighbour is cut-cell
ICUTIDXE = ISHPB(IDMN,IPBCE)
AREA=F(PBC_DAKAR(4)+ICUTIDXE) ! surface area for cut-cell
ELSE ! neighbour is fully blocked
AREA=F(I3DAEX+IJK) ! surface area from cell area
IF(EWP(IJ).AND.EPOR>0) AREA=AREA/(F(L0EPOR+IJ)+TINY)
ENDIF
IPAT=NINT(F(L0PAT+IJKE))
SOLA=1.0; TRANS=0.0; TRANSF=1.0
IF(EF(IJ).AND.IPAT>0) THEN
IF(L0BLKSOLA>0) SOLA=F(L0BLKSOLA+IPAT)
IF(L0TRANS>0) TRANS=F(L0TRANS+IPAT) ! get transparency factor for blockage
ELSE
IF(L0WALSOLA>0) SOLA=F(L0WALSOLA+IJK); IJE=IJ
IF(L0TRANSW>0) TRANS=F(L0TRANSW+IJK) ! get transparency factor for plate
ENDIF
IF(L0TRANSF>0) TRANSF=F(L0TRANSF+IJK)
SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
QDIR=SOLA*QDIR1*TRANSF; QDIF=SOLA*QDIF0
CALL SETSOR(IJ,IJE,IJE,IJE,IJE,IST3,L0SOR3,L0SOR1,QDIR,
1 QDIF,AREA,SUNDIR,NORML,L0SR2,TINY)
IF(LBSOLA>0) F(L0SOLA+IJE)=SOLA
ENDIF
ENDIF
C
IF((WF(IJ).OR.WWP(IJ).OR.CAREA>0.0).AND.IFACE/=4) THEN ! Blocked face on WEST side
CALL VECTOR(NORML,1.,0.,0.) ! set unit vector normal to surface
IF(DOT(NORML,NORM0)>0..OR.CAREA==0.) THEN ! blocked cell is not 'behind' cut cell
IF(DOT(NORML,NORM0)>0..AND.CAREA>0.0) CAREA=0.0
NEXT=ITWO(-NY,0,IX>1)
IJKW=IJK+NEXT ! neighbour index
IPBCW = IPBSEQ(IDMN,IJKW) ! neighbour cut-cell index
IJW=IJ+NEXT
IF(PARSOL.AND.IPBCW>0) THEN ! neighbour is cut-cell
ICUTIDXW = ISHPB(IDMN,IPBCW)
AREA=F(PBC_DAKAR(3)+ICUTIDXW) ! surface area for cut-cell
ELSE ! neighbour is fully blocked
AREA=F(I3DAEX+IJKW) ! surface area from cell area
IF(WWP(IJ).AND.EPOR>0) AREA=AREA/(F(L0EPOR+IJW)+TINY)
ENDIF
IPAT=NINT(F(L0PAT+IJKW))
SOLA=1.0; TRANS=0.0; TRANSF=1.0
IF(WF(IJ).AND.IPAT>0) THEN
IF(L0BLKSOLA>0) SOLA=F(L0BLKSOLA+IPAT)
IF(L0TRANS>0) TRANS=F(L0TRANS+IPAT) ! get transparency factor for blockage
ELSE
IF(L0WALSOLA>0) SOLA=F(L0WALSOLA+IJK); IJW=IJ
IF(L0TRANSW>0) TRANS=F(L0TRANSW+IJK) ! get transparency factor for plate
ENDIF
IF(L0TRANSF>0) TRANSF=F(L0TRANSF+IJK)
SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
QDIR=SOLA*QDIR1*TRANSF; QDIF=SOLA*QDIF0
CALL SETSOR(IJ,IJW,IJW,IJW,IJW,IST3,L0SOR3,L0SOR1,QDIR,
1 QDIF,AREA,SUNDIR,NORML,L0SR2,TINY)
IF(LBSOLA>0) F(L0SOLA+IJW)=SOLA
ENDIF
ENDIF
C
IF((HF(IJ).OR.HWP(IJ).OR.CAREA>0.0).AND.IFACE/=5) THEN ! Blocked face on HIGH side 198/1
CALL VECTOR(NORML,0.,0.,-1.) ! set unit vector normal to surface
IF(DOT(NORML,NORM0)>0..OR.CAREA==0.) THEN ! blocked cell is not 'behind' cut cell
IF(DOT(NORML,NORM0)>0..AND.CAREA>0.0) CAREA=0.0
NEXT=ITWO(NXNY,0,IZZ0) THEN ! neighbour is cut-cell
ICUTIDXH = ISHPB(IDMN,IPBCH)
AREA=F(PBC_DAKAR(6)+ICUTIDXH) ! surface area for cut-cell
ELSE ! neighbour is fully blocked
AREA=F(I3DAHZ+IJK) ! surface area from cell area
IF(HWP(IJ).AND.HPOR>0) AREA=AREA/(F(L0HPOR+IJ)+TINY)
ENDIF
IPAT=NINT(F(L0PAT+IJKH))
SOLA=1.0; TRANS=0.0; TRANSF=1.0
IF(HF(IJ).AND.IPAT>0) THEN
IF(L0BLKSOLA>0) SOLA=F(L0BLKSOLA+IPAT); IADD=NFM
IF(L0TRANS>0) TRANS=F(L0TRANS+IPAT) ! get transparency factor for blockage
ELSE
IF(L0WALSOLA>0) SOLA=F(L0WALSOLA+IJK); IJH=IJ; IADD=0
IF(L0TRANSW>0) TRANS=F(L0TRANSW+IJK) ! get transparency factor for plate
ENDIF
IF(L0TRANSF>0) TRANSF=F(L0TRANSF+IJK)
SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
QDIR=SOLA*QDIR1*TRANSF; QDIF=SOLA*QDIF0
CALL SETSOR(IJ,IJH,IJ+IADDZ,IJ+IADD2,IJ+IADD3,IST3,
1 L0SOR3,L0SOR1,QDIR,QDIF,AREA,SUNDIR,NORML,L0SR2,TINY)
IF(LBSOLA>0) THEN
IF(IJ==IJH) THEN
F(L0SOLA+IJ)=SOLA
ELSE
F(L0SOLA+IJ+IADD)=SOLA
ENDIF
ENDIF
ENDIF
ENDIF
C
IF((LF(IJ).OR.LWP(IJ).OR.CAREA>0.0).AND.IFACE/=6) THEN ! Blocked face on LOW side 198/19
CALL VECTOR(NORML,0.,0.,1.) ! set unit vector normal to surface
IF(DOT(NORML,NORM0)>0..OR.CAREA==0.) THEN ! blocked cell is not 'behind' cut cell
IF(DOT(NORML,NORM0)>0..AND.CAREA>0.0) CAREA=0.0
NEXT=ITWO(-NXNY,0,IZZ>1)
IJKL=IJK+NEXT ! neighbour index
IPBCL = IPBSEQ(IDMN,IJKL) ! neighbour cut-cell index
IJL=IJ+NEXT
IF(PARSOL.AND.IPBCL>0) THEN ! neighbour is cut-cell
ICUTIDXL = ISHPB(IDMN,IPBCL)
AREA=F(PBC_DAKAR(5)+ICUTIDXL) ! surface area for cut-cell
ELSE ! neighbour is fully blocked
AREA=F(I3DAHZ+IJKL) ! surface area from cell area
IF(LWP(IJ).AND.HPOR>0) AREA=AREA/(F(L0HPOR+IJ-NFM)+
1 TINY)
ENDIF
IPAT=NINT(F(L0PAT+IJKL))
SOLA=1.0; TRANS=0.0; TRANSF=1.0
IF(LF(IJ).AND.IPAT>0) THEN
IF(L0BLKSOLA>0) SOLA=F(L0BLKSOLA+IPAT); IADD=-NFM
IF(L0TRANS>0) TRANS=F(L0TRANS+IPAT) ! get transparency factor for blockage
ELSE
IF(L0WALSOLA>0) SOLA=F(L0WALSOLA+IJK); IJL=IJ; IADD=0
IF(L0TRANSW>0) TRANS=F(L0TRANSW+IJK) ! get transparency factor for plate
ENDIF
IF(L0TRANSF>0) TRANSF=F(L0TRANSF+IJK)
SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
QDIR=SOLA*QDIR1*TRANSF; QDIF=SOLA*QDIF0
CALL SETSOR(IJ,IJL,IJ-IADDZ,IJ-IADD2,IJ-IADD3,IST3,
1 L0SOR3,L0SOR1,QDIR,QDIF,AREA,SUNDIR,NORML,L0SR2,TINY)
IF(LBSOLA>0) THEN
IF(IJ==IJL) THEN
F(L0SOLA+IJ)=SOLA
ELSE
F(L0SOLA+IJ+IADD)=SOLA
ENDIF
ENDIF
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
L0FACZ=L0FACZ0
END
C----------------------------------------------------------------------
SUBROUTINE GET_SUN_HEAT(SUNALT,IQDIR,IQDIF,ISKY,QDIR,QDIF,LUPR1)
USE weather_data
COMMON/RGE/RFIL1(2),DT,TIM,RFIL2(11)
1 /LDATA/ LDT1(18),STEADY,LDT2(65)
REAL ACON(6),BCON(6),CCON(2)
LOGICAL LDT1,STEADY,LDT2,GTZ,QGT,QLE
DATA ACON/4.7359E-7, -1.5206E-4, 1.9734E-2, -1.3427E+0, 5.1158E+1,
1 -4.8624E+0/
DATA BCON/1.3872E-7, -3.6705E-5, 3.7394E-3, -1.8993E-1, 5.7408E+0,
1 3.5992E-3/
DATA CCON/5.0591E+0, -1.2638E+0/
X=ASIN(SUNALT)*45./ATAN(1.)
IF(IQDIR==1) THEN ! Direct radiation from solar altitude
IF(GTZ(X)) THEN ! sun above horizon
C... 5th order polynomial fit to Table A2.24 (Idn) CIBSE Guide Volume A
QDIR=ACON(1)*X**5+ACON(2)*X**4+ACON(3)*X**3+ACON(4)*X**2+
1 ACON(5)*X+ACON(6)
ELSE
QDIR=0.0
ENDIF
ELSEIF(IQDIR==2) THEN ! Direct radiation sent from weather file
TIM0=TIM-0.5*DT; QDIR=0.
IF(GTZ(X)) THEN
DO I=1,NHEATS
TIM1=(I-1)*3600/NPHOUR; TIM2=I*3600/NPHOUR
IF(QGT(TIM0,TIM1).AND.QLE(TIM0,TIM2)) THEN
QDIR=QDR(I)+(QDR(I+1)-QDR(I))*(TIM0-TIM1)/(TIM2-TIM1)
EXIT
ENDIF
ENDDO
ELSE
QDIR=0.0
ENDIF
ENDIF
IF(IQDIF==1) THEN ! Diffuse radiation from solar altitude
IF(GTZ(X)) THEN ! sun above horizon
IF(ISKY==0) THEN ! Clear sky
C... 5th order polynomial fit to Table A2.25 (Clear) CIBSE Guide Volume A
QDIF=BCON(1)*X**5+BCON(2)*X**4+BCON(3)*X**3+BCON(4)*X**2+
1 BCON(5)*X+BCON(6)
ELSE ! Cloudy sky
C... Linear fit to Table A2.25 (Cloudy) CIBSE Guide Volume A
QDIF=CCON(1)*X+CCON(2)
ENDIF
ELSE
QDIF=0.0
ENDIF
ELSEIF(IQDIF==2) THEN ! Diffuse radiation sent from weather file
IF(GTZ(X)) THEN ! sun above horizon
TIM0=TIM-0.5*DT; QDIF=0.
DO I=1,NHEATS
TIM1=(I-1)*3600/NPHOUR; TIM2=I*3600/NPHOUR
IF(QGT(TIM0,TIM1).AND.QLE(TIM0,TIM2)) THEN
QDIF=QDF(I)+(QDF(I+1)-QDF(I))*(TIM0-TIM1)/(TIM2-TIM1)
EXIT
ENDIF
ENDDO
ELSE
QDIF=0.0
ENDIF
ENDIF
CALL WRITST
IHR=INT(TIM/3600); TIM1=TIM-IHR*3600; IMN=TIM1/60
IF(.NOT.STEADY)
1 WRITE(LUPR1,'('' Elapsed time: '',I4,'' hrs '',I2,'' mins'')')
1 IHR,IMN
WRITE(LUPR1,'('' Solar altitude: '',1PE13.6,A)') X,CHAR(176)
WRITE(LUPR1,'('' Direct solar radiation: '',1PE13.6,'' W/m^2'')'
1 ) QDIR
WRITE(LUPR1,'('' Diffuse solar radiation: '',1PE13.6,'' W/m^2'')'
1 ) QDIF
IF(IQDIF==1) THEN
IF(ISKY==0) THEN
WRITE(LUPR1,'('' Sky: Clear'')')
ELSE
WRITE(LUPR1,'('' Sky: Cloudy'')')
ENDIF
ENDIF
CALL WRITST
END
C----------------------------------------------------------------------
SUBROUTINE SETSOR(IJ,IJN,IJN1,IJN2,IJN3,IST3,L0SOR3,L0SOR1,QDIR,
1 QDIF,AREA,SUNDIR,NORML,L0SOR2,TINY)
INCLUDE 'farray'
COMMON/RGE/RFIL1(2),DT,TIM,RFIL2(11)
REAL SUNDIR(3),NORML(3)
LOGICAL IST3,SLD,VAC,POR,NEZ
IF(SLD(IJN)) THEN ! neighbour is solid
IF(.NOT.(VAC(IJN).OR.POR(IJN))) THEN
IF(IST3) THEN
L0SOR=L0SOR3
ISOR1=IJN3 ! T3 and neighbour PRPS< 198, source in neighbour
ISOR2=IJN2
ELSE
L0SOR=L0SOR1
ISOR1=IJN1 ! T1 and neighbour PRPS<198, source in neighbour
ISOR2=IJN2
ENDIF
ELSE
ISOR1=IJ ; L0SOR=L0SOR1 ! T1 and PRPS 198/199 source in cell
ISOR2=IJ; IJN=IJ
ENDIF
ELSE
ISOR1=IJ; ISOR2=IJ; L0SOR=L0SOR1 ! Face blocked by plate, source in cell
ENDIF
QSOR=AREA*(QDIR*AMAX1(DOT(SUNDIR,NORML),0.0)+QDIF)
if(l0sor2>0) then
area0=f(l0sor+isor1)/(f(l0sor2+isor2)+tiny)
qtot=f(l0sor+isor1)+qsor; atot=area0+area
f(l0sor2+isor2)=qtot/(atot+tiny)
endif
F(L0SOR+ISOR1)=F(L0SOR+ISOR1)+QSOR
END
C----------------------------------------------------------------------
SUBROUTINE GET_SUN_LIT(NUMSUN,LBLIT,L0LIT,L0LIT0,SUNDIR,INILIT,
1 L0TRANS,L0TRANSF)
INCLUDE 'farray'
INCLUDE 'd_earth/parvar'
INCLUDE 'patnos'
INCLUDE 'objnam'
INCLUDE 'patcmn'
INCLUDE 'satcfile'
INCLUDE 'parear'
COMMON/GENI/NXNY,IGFIL1(8),NFM,IGF(21),IPRL,IBTAU,ILTLS,IGFIL(15),
1 ITEM1,ITEM2,ISPH1,ISPH2,ICON1,ICON2,IPRPS,IRADX,IRADY,IRADZ,IVFOL
COMMON /RDATA/RF1(28),READFI,RF2(56) /RDA8/FIINIT(150)
COMMON /IDATA/ NX,NY,NZ,LUPR1,IDF2(65),I0PAT,NUMPAT,IDF3(49)
COMMON /LDATA/LDATA1(18),STEADY,LDATA2(65)
LOGICAL LDATA1,STEADY,LDATA2
COMMON /INTDMN/IDMN,IDUM(8)
COMMON /SLDPRP/SOLPRP,PORPRP,VACPRP,SOLCRT,SOLSAV
COMMON /WORDI1/ NWDS,NCHARS(20),NSEMI,H,NLINES
INTEGER H
COMMON /WORDL1/ ERROR,MORWDS
COMMON /WORDC1/ WD(20),INLINE
LOGICAL ERROR,MORWDS
CHARACTER WD*20, INLINE*120
CHARACTER*68 CVAL,LINE*256,COBTP*12,COBNAM*12,NAMPT*12
REAL SUNDIR(3),NORML(3)
INTEGER ANYZ,GD_GLOBPAT
LOGICAL QGT,QEQ,QNE,LPAR
C
C
CALL WRYTPB('Start calculating solar shading')
LINE='sunpos.txt facetdat -nographics'
CALL MYRUN_OSGCANSEE(LINE)
CALL WRYTPB('Done calculating solar shading')
C... open output file from OSGCANSEE
OPEN(111,FILE='sunpos.lit',STATUS='OLD',IOSTAT=IOS)
IF(IOS/=0) THEN
WRITE(14,'(''Error number'',I4,'' opening sunpos.lit'')') IOS
ENDIF
READ(111,'(A)',END=101) LINE ! line 1 has sun direction
LOUT=LENGZZ(LINE); CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LOUT)
DO II=1,3
SUNDIR(II)=RRDZZZ(II+1)
ENDDO
CALL NORM(SUNDIR) ! normalise
IF(LBLIT>0) THEN
IF(STEADY.AND.QEQ(FIINIT(LBLIT),READFI)) THEN
CALL WRYTPB('Solar shading read from restart file')
CALL WRITST
CALL WRIT40('Solar shading read from restart file')
CALL WRITST
C... return from here as LIT will have been read from restart file
RETURN
ENDIF
C.... Initialise LIT store to 0 or 1 depending on setting of ALLCELS
DO IZZ=1,NZ
L0LIT=L0F(ANYZ(LBLIT,IZZ))
DO I=1,NXNY
F(L0LIT+I)=INILIT
ENDDO
ENDDO
ELSE
DO IZZ=1,NZ
DO I=1,NXNY
IJK=I+(IZZ-1)*NXNY
F(L0LIT0+IJK)=INILIT
ENDDO
ENDDO
ENDIF
L0PAT=L0PATNO(IDMN) ! index for patch-number store
C... must use global patch number in parallel
LPAR=MIMD.AND.NPROC>1
NUMPT=ITWO(GD_NUMPAT,NUMPAT,LPAR)
DO II=1,NUMSUN ! read results for each cell with surface
READ(111,'(A)') LINE
LOUT=LENGZZ(LINE); CALL SPLTZZ(LINE,WD,NWDS,NCHARS,LOUT)
C... IJK= first word; IHIT=second word. =0-> illuminated, >0 -> not
IJK=IRDZZZ(1); IHIT=IRDZZZ(2)
IF(IHIT>0) THEN ! ray hit something, check if still lit
NHITS=IHIT; IHIT=0
HIT: DO IH=1,NHITS
DO IPAT=1,NUMPT
IF(LPAR) THEN
COBNAM=GD_OBJNAM(IPAT); IF(COBNAM==' ') CYCLE
NAMPT=GD_NAMPAT(IPAT)(2:)
ELSE
COBNAM=OBJNAM(IPAT); IF(COBNAM=='NOTSET') CYCLE
NAMPT=NAMPAT(IPAT)
ENDIF
IF(COBNAM==WD(IH+2)) THEN
COBTP=' '; CALL GETSDC('OBJTYP','^'//NAMPT,COBTP)
IF(COBTP=='BLOCKAGE') THEN
IHIT=IHIT+1
IF(L0TRANS>0) THEN
IF(F(L0TRANS+IPAT)>0.0) THEN ! only allow fully transparent
IHIT=IHIT-1; CYCLE HIT ! might be lit
ENDIF
ENDIF
RMAT=-1.; CALL GETSDR(COBNAM,'MATERIAL',RMAT)
IF(RMATopaque
EXIT HIT ! not lit
ENDIF
ELSEIF(COBTP=='PLATE'.OR.COBTP=='THINPLT') THEN
IHIT=IHIT+1
IF(L0TRANS>0) THEN
IF(F(L0TRANS+IPAT)>0.0) THEN ! only allow fully transparent
IHIT=IHIT-1; CYCLE HIT ! might be lit
ENDIF
ENDIF
EXIT HIT ! not lit
ELSEIF((COBTP=='USER_DEFINED')) THEN
IF(NAMPT(1:4)=='TREE') THEN
IHIT=IHIT+1
IF(L0TRANSF==0) THEN
CALL GXMAKE(L0TRANSF,NX*NY*NZ,'L0TRANSF')
DO III=1,NX*NY*NZ; F(L0TRANSF+III)=1.0; ENDDO
ENDIF
TRANS=0.0; CALL GETSDR(COBNAM,'TRANSPARENCY',TRANS)
F(L0TRANSF+IJK)=TRANS
IF(TRANS>0.0) THEN ! only allow fully transparent
IHIT=IHIT-1; CYCLE HIT ! might be lit
ELSE
IHIT=IHIT+1; EXIT HIT ! found foliage object
ENDIF
ELSE
CYCLE HIT ! object is user_defined. Might be lit
ENDIF
C.... Allow flow boundaries to be transparent
ELSEIF(COBTP=='INLET'.OR.COBTP=='OUTLET'.OR.COBTP=='FAN'
1 ) THEN
IHIT=IHIT+1
C.... Allow for transparency flag
IF(L0TRANS>0) THEN
IF(F(L0TRANS+IPAT)<1.0) THEN ! only allow fully opaque
EXIT HIT ! not lit
ENDIF
ENDIF
IHIT=IHIT-1; CYCLE HIT ! object is INLET/OUTLET/FAN. Might be lit
ELSE ! not plate and not blockage - not lit
IHIT=IHIT+1
IF(L0TRANS>0) THEN
IF(F(L0TRANS+IPAT)>0.0) THEN ! only allow fully transparent
IHIT=IHIT-1; CYCLE HIT ! might be lit
ENDIF
ENDIF
EXIT HIT
ENDIF
ENDIF
ENDDO
ENDDO HIT
ENDIF
IF(LBLIT>0) THEN
IZZ=(IJK-1)/NXNY+1; IZAD=(IZZ-1)*NXNY; IJ=IJK-IZAD
L0LIT=L0F(ANYZ(LBLIT,IZZ))
F(L0LIT+IJ)=ITWO(1,0,IHIT<=0.AND.QGT(SUNDIR(3),0.0))
ELSE
F(L0LIT0+IJK)=ITWO(1,0,IHIT<=0.AND.QGT(SUNDIR(3),0.0))
ENDIF
ENDDO
101 CONTINUE
END
c