c

C.... FILE NAME GXSUN.FTN-------------------------------- 050416
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 '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:
C
      REAL NORML(3),VELIN(3),VELOUT(3),DIR(3),SUNDIR(3)
      PARAMETER (MAXDMN = 10)
      CHARACTER*12 COBNAM,COBNM2,BLOCKNAM,PATNM*8
      CHARACTER*68 CTIME,CVAL,LINE*256,COBTP*12
      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,LBSOR,L0SOR0,LBSOR3,
     1     L0SR30,L0WALSOLA,L0BLKSOLA,L0TRANS,L0TRANSW
      SAVE SUNDIR,QDIR,QDIF
      SAVE IQDIR,IQDIF,ISKY
      DATA DONE /.FALSE./
C
      NAMSUB='GXSUN'
C***********************************************************************
      IF(IGR==1.AND.ISC==2) THEN
        IF(.NOT.NULLPR) THEN
          CALL WRYT40('GXSUN of:                           120816')
          CALL WRITST
          CALL WRIT40('GXSUN of:                           120816')
        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
              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
                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)
                ELSEIF(ICUT>0) THEN
                  XPOS=CUTCELL%CEN(1,ICUT) ! coords of sunny centre
                  YPOS=CUTCELL%CEN(2,ICUT)
                  ZPOS=CUTCELL%CEN(3,ICUT)
                ELSE                 ! un-cut cell
                  XPOS=F(KXXG+IXX)     ! coords of cell centre
                  YPOS=F(KYYG+IYY)
                  ZPOS=F(KZZG+IZZ)
                  EPSX=F(KXDX+IXX)*0.001
                  EPSY=F(KYDY+IYY)*0.001
                  EPSZ=F(KZDZ+IZZ)*0.001
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,'SUNSOR')
        ENDIF
        IF(ITM3>0) THEN
          LBSOR3=LBNAME('#QS3'); L0SR30=0
          IF(LBSOR3==0) THEN
            CALL GXMAKE0(L0SR30,NX*NY*IZSUNL,'SUNSOR3')
          ENDIF
        ENDIF
C... get solar absorbtion factors for plates, and store them
        IF(NMWALL>0) THEN
          L0WALLSOLA=0; L0TRANSW=0
          DO IWL= 1,NMWALL ! loop over wall patches
            IPAT= NINT(ABS(F(L0WALL+IWL)))
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
          ENDDO
        ENDIF
C... must use global patch number in parallel
        NUMPT=ITWO(GD_NUMPAT,NUMPAT,LPAR)
        L0BLKSOLA=0; CALL GXMAKE0(L0TRANS,NUMPT,'TRANSPARENCY') ! allocate storage
        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=='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
        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)
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)
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) 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)
              IF(LBSR>0) THEN
                L0SOR=L0F(LBSR)
              ELSE
                L0SR0=ITWO(L0SOR0,L0SR30,INDVAR==ITEM1)
                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.
        CVAL=' ';CALL GETSDC('SUNLIGHT','TIME',CVAL)
        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
          CALL SPLTZZ(CTIME,WD,NWDS,NCHARS,LT)
          ID=IRDZZZ(3);IHOUR=IRDZZZ(4);IMINS=IRDZZZ(5);ISECS=IRDZZZ(6)
          ITIM=IHOUR*60*60+IMINS*60+ISECS
C... ITIM is time in sec set in Q1. Add elapsed time of run to get current time of day.
          ITIM=ITIM+TIM-DT/2. ! time is at mid-point of time step to get average
          IHOUR=ITIM/3600
          ITIM=ITIM-IHOUR*3600
          IMINS=ITIM/60
          ISECS=ITIM-IMINS*60
          IF(IHOUR>24) THEN
            ID=ID+1; IHOUR=IHOUR-24
          ENDIF
          WRITE(WD(3),'(I2.2)') ID
          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)
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)
  191     CONTINUE
        ENDIF
      ENDIF
C****************************************************************
      END
C----------------------------------------------------------------------
      SUBROUTINE SET_SUN_SOR(LBLIT,L0LIT0,LBSOR,L0SOR0,LBSOR3,L0SR30,
     1      IZ2,SUNDIR,QDIR0,QDIF0,L0WALSOLA,L0BLKSOLA,L0TRANS,L0TRANSW)
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)
      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)
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
      lbsor2=lbname('#QS2'); l0sor2=0
      if(lbsor2>0) then
        do izz=1,nz ! initialise per unit area source stores to zero
          l0sor2=l0f(anyz(lbsor2,izz))
          do ij=1,nxny
            f(l0sor2+ij)=0.0
          enddo
        enddo
      endif
      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))
        ELSE
          L0SOR1=L0SOR0+(IZZ-1)*NXNY
        ENDIF
        IF(IST3) THEN ! store for T3 heat source
          IF(LBSOR3/=0) THEN
            L0SOR3=L0F(ANYZ(LBSOR3,IZZ))
          ELSE
            L0SOR3=L0SR30+(IZZ-1)*NXNY
          ENDIF
        ELSE
          L0SOR3=L0SOR1
        ENDIF
        IF(LBSOR2>0) THEN ! optional store for per-unit-area source
          l0sor2=l0f(anyz(lbsor2,izz))
        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(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
            IF(ISG62==0) THEN
              IF(LCUTCL(IJK)) THEN ! cell is CUT
                NPBC= IPBSEQ(IDMN,IJK)
              ENDIF
            ELSE
              IF(ASSOCIATED(CUTCELL%ICUTH)) THEN
                ICUT=CUTCELL%ICUTH(IJK)
              ELSE
                ICUT=0
              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
            IF(NPBC>0.OR.ICUT>0) THEN
              CALL NORM(NORML) ! unit-normal to cut surface
              ISOR=IJ
              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)
              SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
              L0SOR=ITWO(L0SOR3,L0SOR1,IST3)
              F(L0SOR+ISOR)=F(L0SOR+ISOR)+
     1              SOLA*AREA*(QDIR1*AMAX1(DOT(SUNDIR,NORML),0.0)+QDIF0)
              if(l0sor2>0) f(l0sor2+isor)= f(L0SOR+isor)/(area+tiny)
              if(lbsola>0) f(l0sola+isor)=sola
              CALL COPYVC(NORML,NORM0) ! save cut-cell normal
            ELSE
              CALL VECTOR(NORM0,0.,0.,0.)
            ENDIF
C
            IF(SLD(IJ)) CYCLE ! skip solid cells
            IF    (NF(IJ).OR.NWP(IJ)) THEN ! Blocked face on NORTH side
              IF(ICUT==0.OR.DOT(NORML,NORM0)<0.) THEN ! blocked cell is not 'behind' cut cell
              CALL VECTOR(NORML,0.,-1.,0.) ! set unit vector normal to surface
              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
              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
              SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
              QDIR=SOLA*QDIR1; 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,L0SOR2,TINY)
              if(lbsola>0) f(l0sola+ijn)=sola ! save absorption factor for plotting
              ENDIF
            ENDIF
C
            IF(SF(IJ).OR.SWP(IJ)) THEN ! Blocked face on SOUTH side
              CALL VECTOR(NORML,0.,1.,0.) ! set unit vector normal to surface
              IF(ICUT==0.OR.DOT(NORML,NORM0)<0.) THEN
              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
              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
              SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
              QDIR=SOLA*QDIR1; QDIF=SOLA*QDIF0
              CALL SETSOR(IJ,IJS,IJS,IJS,IJS,IST3,L0SOR3,L0SOR1,QDIR,
     1                           QDIF,AREA,SUNDIR,NORML,L0SOR2,TINY)
              if(lbsola>0) f(l0sola+ijs)=sola
              ENDIF
            ENDIF
C
            IF(EF(IJ).OR.EWP(IJ)) THEN ! Blocked face on EAST side
              CALL VECTOR(NORML,-1.,0.,0.) ! set unit vector normal to surface
              IF(ICUT==0.OR.DOT(NORML,NORM0)<0.) THEN
              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
              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
              SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
              QDIR=SOLA*QDIR1; QDIF=SOLA*QDIF0
              CALL SETSOR(IJ,IJE,IJE,IJE,IJE,IST3,L0SOR3,L0SOR1,QDIR,
     1                           QDIF,AREA,SUNDIR,NORML,L0SOR2,TINY)
              if(lbsola>0) f(l0sola+ije)=sola
              ENDIF
            ENDIF
C
            IF(WF(IJ).OR.WWP(IJ)) THEN ! Blocked face on WEST side
              CALL VECTOR(NORML,1.,0.,0.) ! set unit vector normal to surface
              IF(ICUT==0.OR.DOT(NORML,NORM0)<0.) THEN
              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
              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
              SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
              QDIR=SOLA*QDIR1; QDIF=SOLA*QDIF0
              CALL SETSOR(IJ,IJW,IJW,IJW,IJW,IST3,L0SOR3,L0SOR1,QDIR,
     1                           QDIF,AREA,SUNDIR,NORML,L0SOR2,TINY)
              if(lbsola>0) f(l0sola+ijw)=sola
              ENDIF
            ENDIF
C
            IF(HF(IJ).OR.HWP(IJ)) THEN ! Blocked face on HIGH side 198/199
              CALL VECTOR(NORML,0.,0.,-1.) ! set unit vector normal to surface
              IF(ICUT==0.OR.DOT(NORML,NORM0)<0.) THEN
              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
              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
              SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
              QDIR=SOLA*QDIR1; QDIF=SOLA*QDIF0
              CALL SETSOR(IJ,IJH,IJ+IADDZ,IJ+IADD2,IJ+IADD3,IST3,L0SOR3,
     1            L0SOR1,QDIR,QDIF,AREA,SUNDIR,NORML,L0SOR2,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)) THEN ! Blocked face on LOW side 198/199
              CALL VECTOR(NORML,0.,0.,1.) ! set unit vector normal to surface
              IF(ICUT==0.OR.DOT(NORML,NORM0)<0.) THEN
              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)+TINY)
              ENDIF
              IPAT=NINT(F(L0PAT+IJKL))
              SOLA=1.0; TRANS=0.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
              SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
              QDIR=SOLA*QDIR1; QDIF=SOLA*QDIF0
              CALL SETSOR(IJ,IJL,IJ-IADDZ,IJ-IADD2,IJ-IADD3,IST3,L0SOR3,
     1            L0SOR1,QDIR,QDIF,AREA,SUNDIR,NORML,L0SOR2,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)
      INCLUDE 'farray'
      INCLUDE 'objnam'
      INCLUDE 'patcmn'
      INCLUDE 'satcfile'
      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 /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
      REAL SUNDIR(3),NORML(3)
      INTEGER ANYZ
      LOGICAL QGT,QEQ,QNE
C
      LPR=LENGZZ(PRELPRE)
      LINE=PRELPRE(1:LPR)//'dlls/osgcansee.exe '//
     1             'sunpos.txt facetdat -nographics'
      CALL WRYTPB('Start calculating solar shading')
      CALL WIN_EXEC_AND_WAIT_FOR(LINE,0,IHD)
      CALL WRYTPB('Done calculating solar shading')
C... open output file from OSGCANSEE
      OPEN(111,FILE='sunpos.lit',STATUS='OLD',IOSTAT=IOS)
      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')
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
      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,NUMPAT
              IF(OBJNAM(IPAT)==WD(IH+2)) THEN
                COBTP=' '
                CALL GETSDC('OBJTYP','^'//NAMPAT(IPAT),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(OBJNAM(IPAT),'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(NAMPAT(IPAT)(1:4)=='TREE') THEN
                    IHIT=IHIT+1; EXIT HIT ! found foliage object
                  ELSE
                    CYCLE HIT ! object is user_defined. Might be lit
                  ENDIF
                ELSEIF(COBTP=='INLET'.OR.COBTP=='OUTLET'.OR.COBTP=='FAN'
     1                                                            ) THEN
                  IHIT=IHIT+1
                  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