c

C.... FILE NAME GXSUN.FTN-------------------------------- 190318
C#### JCL 19.03.18 save 1st 2 elements of WD for CTIME, as get corrupted
C####              between visits. Also SAVE INLIT to prevent undefined value
C####              appearing in LIT field.
C#### JCL 09.03.18 get time right for passing to OSGCANSEE when going
C####              past midnight in transient mode.(E2018-0007)
C#### JCL 14.12.17 use global patch number in GET_SUN_LIT
C#### SCM 06.11.17 Use DATMAKER directory not PRELPRE to locate OSGcansee
C#### JCL 22.02.17 corrections to solar heat source for Sparsol
C#### JCL 04.01.17 move solid-cell skip above cut-cell section to avoid heat sources in
C####              'buried' cut cells - e.g DRG case of 04.01.17
C#### JCL 12.08.16 add missing declaration of LVRCEL
C#### JCL 05.04.16 restore lost IT1 index for Parsol/Xparsol. Echo date
C####              of GXSUN to screen / RESULT                
C#### JCL 21.03.16 protection against un-assigned CUTCELL%ICUTH
C#### JCL 18.08.15 allow for 0/1 transparency flag
C#### JCL 17.08.15 allow flow boundaries to be transparent                  
C#### JCL 05.09.14 allow for cut cells in SPARSOL      
C#### JCL 13.02.14 corrections to getting solar absorption factor
C#### JCL 28.01.13 add missing argument L0LIT0 to GET_SUN_LIT, otherwise
C####              L0LIT0 undefined inside GET_SUN_LIT
C#### JCL 06.12.12 must use global patch number in parallel when saving solar
C####              absorption factor
C#### JCL 06.11.12 move reading of illumination file into subroutine so that
C####              same code is used in steady and transient. 
C#### JCL 17.10.12 Allow plate objects to be illuminated. Add solar absorption
C####              factor for each object
C#### SCM 01.05.12 Declare LCUTCL as logical
C#### JCL 22.02.12 minor formatting tidying
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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/objnam'
      INCLUDE '/phoenics/d_includ/patcmn'
      INCLUDE '/phoeclos/d_includ/d_earth/parvar'
      INCLUDE '/phoenics/d_includ/cmndmn'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoeclos/d_includ/d_earth/pbcstr'
      INCLUDE '/phoenics/d_includ/parear'
!      INCLUDE '/phoenics/d_includ/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,CSAV(2)*20
      COMMON /FACES/L0FACE,L0FACZ 
      COMMON /VRMCMN/L0MSKZ
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB
C#### SCM 01.05.12 Declare LCUTCL as logical
      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
C#### JCL 12.08.16 add missing declaration of LVRCEL
      LOGICAL LVRCEL
      SAVE LPAR,DONE,ALLCELLS
      SAVE L0LIT0,LBLIT,NUMSUN,ITM3,IZSUNL,LBSOR,L0SOR0,LBSOR3,
     1     L0SR30,L0WALSOLA,L0BLKSOLA,L0TRANS,L0TRANSW,L0TRANSF
      SAVE SUNDIR,QDIR,QDIF
C#### JCL 19.03.18 add INITLIT and CSAV(2) to SAVE list
      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:                           190318')
          CALL WRITST
          CALL WRIT40('GXSUN of:                           190318')
        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
C#### JCL 05.09.14 allow for cut cells in SPARSOL      
C####              IF(LCUTCL(IJK)) THEN ! cell is cut, so must have some surface
              IF(ISG62==0) THEN
                IPBC = IPBSEQ(IDMN,IJK) ! get cut-cell number
              ELSE
C#### JCL 21.03.16 protection against un-assigned CUTCELL%ICUTH
                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)
C#### JCL 05.09.14 allow for cut cells in SPARSOL      
                NPBC=0; ICUT=0
                IF(ISG62==0) THEN
                  IF(LCUTCL(IJK)) THEN
                    NPBC=NINT(F(NPBIJK(IDMN)+IJK))
                  ENDIF
                ELSE
C#### JCL 21.03.16 protection against un-assigned CUTCELL%ICUTH
                  IF(ASSOCIATED(CUTCELL%ICUTH)) THEN
                    ICUT=CUTCELL%ICUTH(IJK)
                  ELSE
                    ICUT=0
                  ENDIF    
                ENDIF      
C####                IF(LCUTCL(IJK)) THEN ! cut cell
C####                  NPBC= NINT(F(NPBIJK(IDMN)+IJK));IT1=ISHPB(IDMN,NPBC)
                IF(NPBC>0) THEN
C#### JCL 05.04.16 restore lost IT1 index for Parsol/Xparsol                
                  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)
C#### JCL 18.08.15 offset for target location
                  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#### JCL 18.08.15 
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
C#### JCL 09.07.14 don't need +1
C####                    IF(SWP(IJ)) YPOS=F(KYYV+IYY+1)
                  IF(SWP(IJ)) YPOS=F(KYYV+IYY)+EPSY
                  IF(EWP(IJ)) XPOS=F(KXXU+IXX-1)-EPSX
C#### JCL 09.07.14 don't need +1
C####                    IF(WWP(IJ)) XPOS=F(KXXU+IXX+1)
                  IF(WWP(IJ)) XPOS=F(KXXU+IXX)+EPSX
                  IF(HWP(IJ)) ZPOS=F(KZZW+IZZ-1)-EPSZ
C#### JCL 09.07.14 don't need +1
C####                    IF(LWP(IJ)) ZPOS=F(KZZW+IZZ+1)
                  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')
C####        LBLIT=LBNAME('LIT'); L0LIT0=0 ! 3D store for illuminated cells
        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
C#### JCL 18.08.15 only create storage for solar absorption factor and
C####              transparency factor if non-default values found.
C####          CALL GXMAKE0(L0WALSOLA,NX*NY*NZ,'WALLSOL')
          L0WALLSOLA=0; L0TRANSW=0
          walls: DO IWL= 1,NMWALL ! loop over wall patches
            IPAT= NINT(ABS(F(L0WALL+IWL)))
            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
C#### JCL 09.07.14 +NZ should be +NXNY
C####              IZADD=(IZZ-1)*NZ
                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  
C#### JCL 18.08.15 Allow for transparency flag        
            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 walls 
        ENDIF
C... must use global patch number in parallel        
        NUMPT=ITWO(GD_NUMPAT,NUMPAT,LPAR)
C#### JCL 18.08.15 only allocate L0BLKSOLA if any factor/=1. Allow for transparency flag        
C####        CALL GXMAKE0(L0BLKSOLA,NUMPT,'BLKSOL')
        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
C#### JCL 18.08.15 Allow for transparency flag
          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 (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) 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
C#### JCL 09.03.18 get time right for passing to OSGCANSEE when going
C####              past midnight in transient mode.
          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)
C#### JCL 19.03.18 save 1st 2 elements of WD for later use
            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
C#### JCL 19.03.18 restore 1st two elements of WD, as corrupted between visits
          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   
C#### JCL 28.01.13 add missing argument L0LIT0 to GET_SUN_LIT
C#### JCL 18.08.15 Allow for transparency flag        
          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          
C#### JCL 18.08.15 Allow for transparency flag        
          CALL SET_SUN_SOR(LBLIT,L0LIT0,LBSOR,L0SOR0,LBSOR3,L0SR30,
     1     IZSUNL,SUNDIR,QDIR,QDIF,L0WALSOLA,L0BLKSOLA,L0TRANS,L0TRANSW,
     1                                                        L0TRANSF)
  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,
     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 '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/patnos'
      INCLUDE '/phoenics/d_includ/cmndmn'
      INCLUDE '/phoenics/d_includ/satear'
      INCLUDE '/phoenics/d_includ/grdloc'
      INCLUDE '/phoenics/d_includ/satgrd'
      INCLUDE '/phoenics/d_includ/grdear'
      INCLUDE '/phoeclos/d_includ/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)
      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
      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)) ! 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)
        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
C#### JCL 04.01.17 move solid-cell skip above cut-cell section to avoid heat sources in
C####              'buried' cut cells - e.g DRG case of 04.01.17
            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
            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%ICUTH(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 
              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)
              TRANSF=1.0; IF(L0TRANSF>0) TRANSF=F(L0TRANSF+IJK) 
              SOLA=SOLA*(1.0-TRANS) ! 1-TRANS is opaqueness
C####              L0SOR=ITWO(L0SOR3,L0SOR1,IST3)
              SOR=SOLA*AREA*(QDIR1*TRANSF*
     1                               AMAX1(DOT(SUNDIR,NORML),0.0)+QDIF0)
              IF(NPBC>0) THEN ! Parsol / Xparsol
C####                F(L0SOR+ISOR)=F(L0SOR+ISOR)+
C####     1              SOLA*AREA*(QDIR1*AMAX1(DOT(SUNDIR,NORML),0.0)+QDIF0)
                IF(L0SOR2>0) THEN ! recover summed area
                  AREA0=F(L0SOR+ISOR)/(F(L0SOR2+IJ)+TINY)
                  AREA=AREA+AREA0 ! sum current area
                ENDIF
                F(L0SOR+ISOR)=F(L0SOR+ISOR)+SOR
                IF(L0SOR2>0) F(L0SOR2+IJ)=F(L0SOR+ISOR)/(AREA+TINY)
                IF(LBSOLA>0) F(L0SOLA+IJ)=SOLA
              ELSE            ! Sparsol
                IMAXDIR=MAXLOC(ABS(NORML)) ! biggest normal
                IF(IMAXDIR(1)==1) THEN     ! mainly X, E
                  IAD1=NY; IAD2=IAD1; IAD3=IAD1; IADS=IAD1
                ELSEIF(IMAXDIR(1)==2) THEN ! mainly Y, N
                  IAD1=1; IAD2=IAD1; IAD3=IAD1; IADS=IAD1
                ELSE                       ! mainly Z, H 
                  IAD1=ITWO(IADD3,IADDZ,IST3)
                  IAD2=ITWO(NXNY,NFM,LBSOR2==0)
                  IAD3=NXNY; IADS=NFM
                ENDIF
                IF(NORML(IMAXDIR(1))>0.0) THEN ! switch to W, S, L
                  IAD1=-IAD1; IAD2=-IAD2; IAD3=-IAD3; IADS=-IADS
                ENDIF
                ISOR=IJ+IAD1; ISOR2=IJ+IAD2
                IF(SLD(IJ+IAD3)) THEN ! nieghbour is solid
                  IF(L0SOR2>0) THEN ! recover summed area
                    AREA0=F(L0SOR+ISOR)/(F(L0SOR2+ISOR2)+TINY)
                    AREA=AREA+AREA0 ! sum current area
                  ENDIF
                  F(L0SOR+ISOR)=F(L0SOR+ISOR)+SOR
                  IF(L0SOR2>0) F(L0SOR2+ISOR2)=F(L0SOR+ISOR)/(AREA+TINY)
                  IF(LBSOLA>0) F(L0SOLA+IJ+IADS)=SOLA
                ENDIF
              ENDIF
              CALL COPYVC(NORML,NORM0) ! save cut-cell normal
            ELSE
              CALL VECTOR(NORM0,0.,0.,0.)  
            ENDIF
C
C#### JCL 04.01.17 move solid-cell skip above cut-cell section to avoid heat sources in
C####              'buried' cut cells - e.g DRG case of 04.01.17
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; 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,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; 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,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; 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,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; 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,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; 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,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; 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,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 '/phoenics/d_includ/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
C#### JCL 09.07.14 next line is duplicate
C####          L0SOR=L0SOR1 ! T1 or 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----------------------------------------------------------------------
C#### JCL 18.08.15 Allow for transparency flag        
      SUBROUTINE GET_SUN_LIT(NUMSUN,LBLIT,L0LIT,L0LIT0,SUNDIR,INILIT,
     1                                                 L0TRANS,L0TRANSF)      
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/patnos'
      INCLUDE '/phoenics/d_includ/objnam'
      INCLUDE '/phoenics/d_includ/patcmn'
      INCLUDE '/phoenics/d_includ/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 /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
      REAL SUNDIR(3),NORML(3)
      INTEGER ANYZ,GD_GLOBPAT
      LOGICAL QGT,QEQ,QNE
C
C#### SCM 06.11.17 Use DATMAKER directory not PRELPRE to locate OSGcansee
C####      LPR=LENGZZ(PRELPRE)
C####      LINE=PRELPRE(1:LPR)//'dlls/osgcansee.exe '//
      LLD=LENGZZ(DATMAKER); LLE=LENGZZ(EXESUF)
      LINE=DATMAKER(1:LLD)//'osgcansee'//EXESUF(:LLE)//
     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  
      L0PAT=L0PATNO(IDMN) ! index for patch-number store
      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
C#### JCL 14.12.17 recover global patch number
                IPTG=GD_GLOBPAT(IPAT) ! get global patch number (same as local for seq)
                COBTP=' '; CALL GETSDC('OBJTYP','^'//NAMPAT(IPAT),COBTP)
                IF(COBTP=='BLOCKAGE') THEN
                  IHIT=IHIT+1
C#### JCL 18.08.15 Allow for transparency flag        
                  IF(L0TRANS>0) THEN 
C#### JCL 14.12.17 use global patch number
                    IF(F(L0TRANS+IPTG)>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
C#### JCL 09.07.14 improve treatment of plates
C####               ELSEIF((COBTP=='PLATE'.OR.COBTP=='THINPLT').AND.
C####     1                                                    NHITS==1) THEN
C####                    IHIT=0; CYCLE HIT ! object is plate. Might be lit
                ELSEIF(COBTP=='PLATE'.OR.COBTP=='THINPLT') THEN
                  IHIT=IHIT+1
C#### JCL 18.08.15 Allow for transparency flag        
                  IF(L0TRANS>0) THEN 
C#### JCL 14.12.17 use global patch number
                    IF(F(L0TRANS+IPTG)>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
!                    IHIT=-(abs(IHIT)+1); cycle HIT ! found foliage object
                    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.5; CALL GETSDR(OBJNAM(IPAT),'TRANSPARENCY',
     1                                                            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#### JCL 17.08.15 allow flow boundaries to be transparent                  
                ELSEIF(COBTP=='INLET'.OR.COBTP=='OUTLET'.OR.COBTP=='FAN'
     1                                                            ) THEN
                  IHIT=IHIT+1
C#### JCL 18.08.15 Allow for transparency flag        
                  IF(L0TRANS>0) THEN 
C#### JCL 14.12.17 use global patch number
                    IF(F(L0TRANS+IPTG)<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 
C#### JCL 14.12.17 use global patch number
                    IF(F(L0TRANS+IPTG)>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