c

C.... FILE NAME GXSUN.FTN-------------------------------- 220212
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   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------------------------------------------------------------------------------------
      SUBROUTINE GXSUN
      USE weather_data
      INCLUDE '/phoenics/d_includ/farray'
      INCLUDE '/phoenics/d_includ/objnam'
      INCLUDE '/phoenics/d_includ/patcmn'
      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,IF01(3),KYDY,KYYG,IF02(9),KZDZ,KZZG,IF03(287)
      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/SLDPRP/SOLPRP,PORPRP,VACPRP,SOLCRT,SOLSAV
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
      COMMON /FACES/L0FACE,L0FACZ 
      COMMON /VRMCMN/L0MSKZ
      COMMON/NAMFN/NAMFUN,NAMSUB
      CHARACTER*6 NAMFUN,NAMSUB,COBTP*12
      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
      SAVE LPAR,DONE
      SAVE L0LIT0,LBLIT,NUMSUN,ITM3,IZSUNL,LBSOR,L0SOR0,LBSOR3,
     1     L0SR30
      SAVE SUNDIR,QDIR,QDIF
      SAVE IQDIR,IQDIF,ISKY
      DATA DONE /.FALSE./
C
      NAMSUB='GXSUN'
C***********************************************************************
      IF(IGR.EQ.1.AND.ISC.EQ.2) THEN
C... set index for T3
        ITM3=LBNAME('T3'); LBSOR3=0
        IF(NUMDMN.GT.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.GT.1) CALL INDXDM
        NCELL=0
        LPAR=MIMD.AND.NPROC.GT.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.NE.' ') THEN
          CTIME=CVAL; LT=LENGZZ(CTIME)
          DO IC=1,LT
            IF(CTIME(IC:IC).EQ.'/') CTIME(IC:IC)=' '
            IF(CTIME(IC:IC).EQ.'.') 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.NE.' '.AND..NOT.STEADY) THEN
            CALL GETSDI('SUNLIGHT','NHEATS',NHEATS)
            NPHOUR=1; CALL GETSDI('SUNLIGHT','NPHOUR',NPHOUR)
            IF(NHEATS.GT.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.EQ.' '.OR.CVAL.EQ.'CONSTANT') THEN
              IQDIR=0; QDIR=1000.;CALL GETSDR('SUNLIGHT','DIRHEAT',QDIR)
            ELSEIF(CVAL.EQ.'ALTITUDE') THEN
              IQDIR=1
            ENDIF    
            CVAL=' '; CALL GETSDC('SUNLIGHT','DIFFUSE',CVAL)
            IF(CVAL.EQ.' '.OR.CVAL.EQ.'CONSTANT') THEN
              IQDIF=0; QDIF=100.; CALL GETSDR('SUNLIGHT','DIFHEAT',QDIF)
            ELSEIF(CVAL.EQ.'ALTITUDE') THEN
              IQDIF=1
              CVAL='CLEAR';CALL GETSDC('SUNLIGHT','SKY',CVAL)
              ISKY=ITWO(0,1,CVAL.EQ.'CLEAR')
            ENDIF    
          ENDIF  
          CALL WRITST
          CALL WRIT40(' Solar Heating data')
          CALL WRIT40(' ------------------')
          IF(EPWNAM.NE.' ') 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.EQ.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.EQ.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.EQ.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.EQ.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
        LBSRF=LBNAME('SRF'); L0SRF=0
        L0FACZ0=L0FACZ; IZ=1; NUMSUN=0; IREC=4
        IZSUNL=-(NZ+1) ! this will be largest IZ with a surface
        DO IZZ=1,NZ ! loop over whole domain looking for surfaces
          IZZM1=(IZZ-1)*NXNY; L0FACZ= L0FACE+IZZM1
          IF(LBSRF.GT.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(LCUTCL(IJK)) 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  
              ENDIF
              IF(IBLK.GT.0) THEN ! found some surface 
                IZSUNL=MAX(IZZ,IZSUNL)
                IF(LCUTCL(IJK)) THEN ! cut cell
                  NPBC= NINT(F(NPBIJK(IDMN)+IJK));IT1 = ISHPB(IDMN,NPBC)
                  XPOS=F(IT1+PBC_SNYX) ! coords of sunny centre
                  YPOS=F(IT1+PBC_SNYY)
                  ZPOS=F(IT1+PBC_SNYZ)
                ELSE ! un-cut cell
                  XPOS=F(KXXG+IXX) ! coords of cell centre
                  YPOS=F(KYYG+IYY)
                  ZPOS=F(KZZG+IZZ)
                ENDIF  
                WRITE(111,REC=IREC,FMT='(I9,3(1PE13.6))') 
     1                               IJK,XPOS,YPOS,ZPOS ! write to output file
                IREC=IREC+1
                IF(L0SRF.GT.0) F(L0SRF+IJ)=1. ! save surface flag if stored
                NUMSUN=NUMSUN+1 ! count number of surfaces
              ELSE ! no surface in this cell
                IF(L0SRF.GT.0) F(L0SRF+IJ)=0. ! save surface flag if stored
              ENDIF  
            ENDDO
          ENDDO
        ENDDO
        IZSUNL=MAX(0,MIN(IZSUNL+1,NZ))
        CLOSE(111,STATUS='KEEP')
        LBLIT=LBNAME('LIT'); L0LIT0=0 ! 3D store for illuminated cells
        IF(LBLIT.EQ.0) THEN
          CALL GXMAKE0(L0LIT0,NX*NY*NZ,'SUNLIT')
        ENDIF  
        LBSOR=LBNAME('#QS1'); L0SOR0=0
        IF(LBSOR.EQ.0) THEN
          CALL GXMAKE0(L0SOR0,NX*NY*IZSUNL,'SUNSOR')
        ENDIF  
        IF(ITM3.GT.0) THEN
          LBSOR3=LBNAME('#QS3'); L0SR30=0
          IF(LBSOR3.EQ.0) THEN
            CALL GXMAKE0(L0SR30,NX*NY*IZSUNL,'SUNSOR3')
          ENDIF  
        ENDIF  
        IF(.NOT.STEADY) RETURN ! for transient get lit cells at start of step
C... Invoke OSGCANSEE to get list of illuminated points        
        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
        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.GT.0) THEN ! ray hit something, check if still lit
            NHITS=IHIT
            HIT: DO IH=1,NHITS
              DO IPAT=1,NUMPAT
                IF(OBJNAM(IPAT).EQ.WD(IH+2)) THEN
                  COBTP=' '
                  CALL GETSDC('OBJTYP','^'//NAMPAT(IPAT),COBTP)
                  IF(COBTP.EQ.'BLOCKAGE') THEN
                    RMAT=-1.; CALL GETSDR(OBJNAM(IPAT),'MATERIAL',RMAT)
                    IF(RMAT.LT.SOLPRP) THEN ! blockage is fluid
                      IHIT=0; CYCLE HIT     ! might be lit
                    ELSE                    ! blockage is solid->opaque
                      IHIT=NHITS; EXIT HIT  ! not lit
                    ENDIF  
                  ELSEIF(COBTP.EQ.'PLATE') THEN
                    IHIT=0; CYCLE HIT ! object is plate. Might be lit
                  ELSE  ! not plate and not blockage - not lit
                    IHIT=NHITS; EXIT HIT
                  ENDIF
                ENDIF
              ENDDO  
            ENDDO HIT
          ENDIF    
          IF(LBLIT.GT.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.EQ.0.AND.QGT(SUNDIR(3),0.0))
          ELSE  
            F(L0LIT0+IJK)=ITWO(1,0,IHIT.EQ.0.AND.QGT(SUNDIR(3),0.0))
          ENDIF  
        ENDDO  
        CALL GET_SUN_HEAT(SUNDIR(3),IQDIR,IQDIF,ISKY,QDIR,QDIF,LUPR1)
        CALL SET_SUN_SOR(LBLIT,L0LIT0,LBSOR,L0SOR0,LBSOR3,L0SR30,IZSUNL,
     1                                                SUNDIR,QDIR,QDIF)
C ... loop back for next FGV domain
  101   CONTINUE
        IF(IDMN.LT.NUMDMN) THEN
          IDMN=IDMN+1
          GO TO 100
        ENDIF
C... reset indices for domain 1
        IF(NUMDMN.GT.1) 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.EQ.13.AND.ISC.EQ.13) THEN
C------------------- SECTION 14 ------------------- value = GRND1
C
        IF(NPATCH.EQ.'SUN') THEN
          IF(INDVAR.EQ.ITEM1.OR.INDVAR.EQ.ITM3) THEN
            IF(IZ.GT.IZSUNL) THEN ! above highest surface, just set zero
              CALL FN1(VAL,0)
            ELSE ! copy source into VAL
              L0VAL=L0F(VAL)
              LBSR=ITWO(LBSOR,LBSOR3,INDVAR.EQ.ITEM1)
              IF(LBSR.GT.0) THEN
                L0SOR=L0F(LBSR)
              ELSE
                L0SR0=ITWO(L0SOR0,L0SR30,INDVAR.EQ.ITEM1)
                L0SOR=L0SR0+(IZ-1)*NXNY
              ENDIF  
              DO I=1,NXNY
                F(L0VAL+I)=F(L0SOR+I)
              ENDDO
            ENDIF      
          ENDIF
        ENDIF  
      ELSEIF(IGR.EQ.19.AND.ISC.EQ.1) THEN
C   * ------------------- SECTION 1 ---- Start of time step.
        CVAL=' ';CALL GETSDC('SUNLIGHT','TIME',CVAL)
        IF(CVAL.NE.' ') THEN
          CTIME=CVAL; LT=LENGZZ(CTIME)
          DO IC=1,LT
            IF(CTIME(IC:IC).EQ.'/') CTIME(IC:IC)=' '
            IF(CTIME(IC:IC).EQ.'.') 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.GT.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   
          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) 
          REWIND(111)
          READ(111,'(A)',END=191) 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
          LBLIT=LBNAME('LIT') ! 3D store for illuminated cells
          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.GT.0) THEN ! ray hit something, check if still lit
              NHITS=IHIT
              HIT1: DO IH=1,NHITS
                DO IPAT=1,NUMPAT
                  IF(OBJNAM(IPAT).EQ.WD(IH+2)) THEN
                    COBTP=' '
                    CALL GETSDC('OBJTYP','^'//NAMPAT(IPAT),COBTP)
                    IF(COBTP.EQ.'BLOCKAGE') THEN
                      RMAT=-1.;CALL GETSDR(OBJNAM(IPAT),'MATERIAL',RMAT)
                      IF(RMAT.LT.SOLPRP) THEN ! blockage is fluid
                        IHIT=0; CYCLE HIT1    ! might be lit
                      ELSE                    ! blockage is solid->opaque
                        IHIT=NHITS; EXIT HIT1 ! not lit
                      ENDIF  
                    ELSEIF(COBTP.EQ.'PLATE') THEN
                      IHIT=0; CYCLE HIT1      ! object is plate. Might be lit
                    ELSE  ! not plate and not blockage - not lit
                      IHIT=NHITS; EXIT HIT1
                    ENDIF
                  ENDIF
                ENDDO  
              ENDDO HIT1
            ENDIF    
            IF(LBLIT.GT.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.EQ.0.AND.QGT(SUNDIR(3),0.0))
            ELSE  
              F(L0LIT0+IJK)=ITWO(1,0,IHIT.EQ.0.AND.QGT(SUNDIR(3),0.0))
            ENDIF  
          ENDDO  
          CALL GET_SUN_HEAT(SUNDIR(3),IQDIR,IQDIF,ISKY,QDIR,QDIF,LUPR1)
          CALL SET_SUN_SOR(LBLIT,L0LIT0,LBSOR,L0SOR0,LBSOR3,L0SR30,
     1                                          IZSUNL,SUNDIR,QDIR,QDIF)
  191     CONTINUE
        ENDIF    
      ENDIF
C****************************************************************
      END
C----------------------------------------------------------------------
      SUBROUTINE SET_SUN_SOR(LBLIT,L0LIT0,LBSOR,L0SOR0,LBSOR3,L0SR30,
     1                                            IZ2,SUNDIR,QDIR0,QDIF)
      INCLUDE '/phoenics/d_includ/farray'
      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 
      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
      REAL SUNDIR(3),NORML(3)
C      
      IST3=LBSOR3.GT.0.OR.L0SR30.GT.0
      LBSOR=LBNAME('#QS1'); L0SOR=0
      ILIM=ITWO(NZ,IZ2,LBSOR.GT.0)
      DO IZZ=1,ILIM ! initialise TEM1 source stores to zero
        IF(LBSOR.GT.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.GT.0)
        DO IZZ=1,ILIM ! initialise T3 source stores to zero
          IF(LBSOR3.GT.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.gt.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  
      L0FACZ0=L0FACZ; L0SOR2=0
      DO IZZ=1,IZ2 ! loop over slabs containing surface
        IF(LBLIT.NE.0) THEN ! store containing lit flag
          L0LIT=L0F(ANYZ(LBLIT,IZZ))
        ELSE
          L0LIT=L0LIT0+(IZZ-1)*NXNY
        ENDIF  
        IF(LBSOR.NE.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.NE.0) THEN
            L0SOR3=L0F(ANYZ(LBSOR3,IZZ))
          ELSE  
            L0SOR3=L0SR30+(IZZ-1)*NXNY
          ENDIF
        ELSE  
          L0SOR3=L0SOR1
        ENDIF
        IF(LBSOR2.GT.0) THEN
          l0sor2=l0f(anyz(lbsor2,izz))
        ENDIF  
        IZZM1=(IZZ-1)*NXNY; L0FACZ= L0FACE+IZZM1
        IADDZ=ITWO(NXNY,NFM,LBSOR.EQ.0)
        IADD2=ITWO(NXNY,NFM,LBSOR2.EQ.0)
        IADD3=ITWO(NXNY,NFM,LBSOR3.EQ.0)
        DO IX=IXF,IXL
          DO IY=IYF,IYL
            IJ=(IX-1)*NY+IY
            IJK=IJ+IZZM1
            IF(F(L0LIT+IJ).GT.0.1) THEN ! the cell is illuminated
              QDIR=QDIR0
            ELSE
              QDIR=0.0
            ENDIF    
            IF(LCUTCL(IJK)) THEN ! cell is CUT
              NPBC= IPBSEQ(IDMN,IJK)
              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))
              CALL NORM(NORML) ! unit-normal to cut surface
              ISOR=IJ
              L0SOR=ITWO(L0SOR3,L0SOR1,IST3)
              F(L0SOR+ISOR)=F(L0SOR+ISOR)+
     1                     AREA*(QDIR*AMAX1(DOT(SUNDIR,NORML),0.0)+QDIF)
              if(l0sor2.gt.0) f(l0sor2+isor)= f(L0SOR+isor)/(area+tiny)                  
            ENDIF
            IF(SLD(IJ)) CYCLE ! skip solid cells
            IF    (NF(IJ).OR.NWP(IJ)) THEN ! Blocked face on NORTH side
              CALL VECTOR(NORML,0.,-1.,0.) ! set unit vector normal to surface
              NEXT=ITWO(1,0,IX.LT.NX)
              IJKN=IJK+NEXT              ! neighbour index
              IPBCN = IPBSEQ(IDMN,IJKN)  ! neighbour cut-cell index
              IF(PARSOL.AND.IPBCN.GT.0) 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
              ENDIF
              IJN=IJ+NEXT
              CALL SETSOR(IJ,IJN,IJN,IJN,IJN,IST3,L0SOR3,L0SOR1,QDIR,
     1                           QDIF,AREA,SUNDIR,NORML,L0SOR2,TINY)
            ENDIF
            IF(SF(IJ).OR.SWP(IJ)) THEN ! Blocked face on SOUTH side
              CALL VECTOR(NORML,0.,1.,0.) ! set unit vector normal to surface
              NEXT=ITWO(-1,0,IX.GT.1)
              IJKS=IJK+NEXT              ! neighbour index
              IPBCS = IPBSEQ(IDMN,IJKS)  ! neighbour cut-cell index
              IF(PARSOL.AND.IPBCS.GT.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
              ENDIF
              IJS=IJ+NEXT
              CALL SETSOR(IJ,IJS,IJS,IJS,IJS,IST3,L0SOR3,L0SOR1,QDIR,
     1                           QDIF,AREA,SUNDIR,NORML,L0SOR2,TINY)
            ENDIF
            IF(EF(IJ).OR.EWP(IJ)) THEN ! Blocked face on EAST side
              CALL VECTOR(NORML,-1.,0.,0.) ! set unit vector normal to surface
              NEXT=ITWO(NY,0,IY.LT.NY)
              IJKE=IJK+NEXT              ! neighbour index
              IPBCE = IPBSEQ(IDMN,IJKE)  ! neighbour cut-cell index
              IF(PARSOL.AND.IPBCE.GT.0) 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
              ENDIF
              IJE=IJ+NEXT
              CALL SETSOR(IJ,IJE,IJE,IJE,IJE,IST3,L0SOR3,L0SOR1,QDIR,
     1                           QDIF,AREA,SUNDIR,NORML,L0SOR2,TINY)
            ENDIF
            IF(WF(IJ).OR.WWP(IJ)) THEN ! Blocked face on WEST side
              CALL VECTOR(NORML,1.,0.,0.) ! set unit vector normal to surface
              NEXT=ITWO(-NY,0,IY.GT.1)
              IJKW=IJK+NEXT              ! neighbour index
              IPBCW = IPBSEQ(IDMN,IJKW)  ! neighbour cut-cell index
              IF(PARSOL.AND.IPBCW.GT.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
              ENDIF
              IJW=IJ+NEXT
              CALL SETSOR(IJ,IJW,IJW,IJW,IJW,IST3,L0SOR3,L0SOR1,QDIR,
     1                           QDIF,AREA,SUNDIR,NORML,L0SOR2,TINY)
            ENDIF
            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
              NEXT=ITWO(NXNY,0,IZ.LT.NZ)
              IJKH=IJK+NEXT              ! neighbour index
              IPBCH = IPBSEQ(IDMN,IJKH)  ! neighbour cut-cell index
              IF(PARSOL.AND.IPBCH.GT.0) 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
              ENDIF
              IJH=IJ+NEXT
              CALL SETSOR(IJ,IJH,IJ+IADDZ,IJ+IADD2,IJ+IADD3,IST3,L0SOR3,
     1            L0SOR1,QDIR,QDIF,AREA,SUNDIR,NORML,L0SOR2,TINY)
            ENDIF
            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
              NEXT=ITWO(-NXNY,0,IZZ.GT.1)
              IJKL=IJK+NEXT              ! neighbour index
              IPBCL = IPBSEQ(IDMN,IJKL)  ! neighbour cut-cell index
              IF(PARSOL.AND.IPBCL.GT.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
              ENDIF
              IJL=IJ+NEXT
              CALL SETSOR(IJ,IJL,IJ-IADDZ,IJ-IADD2,IJ-IADD3,IST3,L0SOR3,
     1            L0SOR1,QDIR,QDIF,AREA,SUNDIR,NORML,L0SOR2,TINY)
            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.EQ.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.EQ.2) THEN  
        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.EQ.1) THEN   ! Diffuse radiation from solar altitude
        IF(GTZ(X)) THEN ! sun above horizon
          IF(ISKY.EQ.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.EQ.2) THEN  
        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.EQ.1) THEN
        IF(ISKY.EQ.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 or PRPS 198/199 source in cell
          L0SOR=L0SOR1 ! T1 or PRPS 198/199 source in cell
          ISOR2=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.gt.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