SUBROUTINE BACK_INT_GRID3D (NX,NY,NZ, NA, NPTS, NCELLS, . GRIDPTR, NEIGHPTR, TREEPTR, CELLFLAGS, . GRIDPOS, XGRID, YGRID, ZGRID, SWEEPORD, . MU, PHI, TRANSMIN, BCFLAG, IPFLAG, . EXTINCT, NSTOKES, SOURCE, KANG, GRIDRAD) C Sweeps through the spatial grid computing radiance vectors for the C discrete ordinate direction (MU,PHI) by integrating the source C function and extinction backward from each point. The sweeping C order through the grid points is in SWEEPORD. If a point is C already done (nonnegative in GRIDRAD) then it is skipped. C For each grid point a ray opposite to the ordinate direction is traced C back to a face containing known (valid) radiances. As each cell is C crossed (usually only one) the integration of the source function C across the cell is done and added to the radiance. The radiance C at the end of the ray path is interpolated from the known grid point C radiances on the face. The resulting radiances are put in GRIDRAD. IMPLICIT NONE INTEGER NX, NY, NZ, NA, NPTS, NCELLS, NSTOKES, KANG INTEGER BCFLAG, IPFLAG INTEGER GRIDPTR(8,NCELLS), NEIGHPTR(6,NCELLS), TREEPTR(2,NCELLS) INTEGER SWEEPORD(NPTS,*) INTEGER*2 CELLFLAGS(NCELLS) REAL GRIDPOS(3,NPTS), XGRID(*), YGRID(*), ZGRID(NZ) REAL MU, PHI, TRANSMIN REAL EXTINCT(NPTS), SOURCE(NSTOKES,NA,NPTS) REAL GRIDRAD(NSTOKES,NPTS) INTEGER BITX, BITY, BITZ, IOCT, JOCT INTEGER IPCELL, ICELL, INEXTCELL, IPT, IFACE INTEGER I1, I2, I3, I4, IOPP, JFACE, KFACE, IC, IZ INTEGER GRIDFACE(4,6), OPPFACE(6) INTEGER IORDER, ICORNER, JOCTORDER3(8), JOCTORDER2(8) LOGICAL VALIDRAD, VALIDFACE, IPINX, IPINY, BTEST DOUBLE PRECISION PI, CX, CY, CZ, CXINV, CYINV, CZINV DOUBLE PRECISION XE, YE, ZE DOUBLE PRECISION SO, SOX, SOY, SOZ, EPS, path DOUBLE PRECISION U, V, F1, F2, F3, F4 DOUBLE PRECISION EXT, EXT0, EXT1 DOUBLE PRECISION SRC(NSTOKES), SRCEXT0(NSTOKES), SRCEXT1(NSTOKES) DOUBLE PRECISION EXT0P, SRCEXT0P(NSTOKES) DOUBLE PRECISION TAU, TRANSCELL, ABSCELL, TRANSMIT DOUBLE PRECISION RAD(NSTOKES), RAD0(NSTOKES) DATA GRIDFACE/1,3,5,7, 2,4,6,8, 1,2,5,6, 3,4,7,8, . 1,2,3,4, 5,6,7,8/, OPPFACE/2,1,4,3,6,5/ DATA JOCTORDER3/1,3,5,7,2,4,6,8/, JOCTORDER2/1,3,1,3,2,4,2,4/ EPS = 1.0E-3*(GRIDPOS(3,GRIDPTR(8,1))-GRIDPOS(3,GRIDPTR(1,1))) C Make the ray direction (opposite to the discrete ordinate direction) PI = ACOS(-1.0) CX = SQRT(1.0-MU**2)*COS(PHI+PI) CY = SQRT(1.0-MU**2)*SIN(PHI+PI) CZ = -MU IF (ABS(CX) .GT. 1.0E-5) THEN CXINV = 1.0D0/CX ELSE CX = 0.0 CXINV = 1.0E6 ENDIF IF (ABS(CY) .GT. 1.0E-5) THEN CYINV = 1.0D0/CY ELSE CY = 0.0 CYINV = 1.0E6 ENDIF CZINV = 1.0D0/CZ C Setup the sweeping direction and the gridpoint location C BITc is 0 for positive direction trace back ray, 1 for negative. IF (CX .LT. 0.0) THEN BITX = 1 ELSE BITX = 0 ENDIF IF (CY .LT. 0.0) THEN BITY = 1 ELSE BITY = 0 ENDIF IF (CZ .LT. -1.0E-3) THEN BITZ = 1 ELSE IF (CZ .GT. 1.0E-3) THEN BITZ = 0 ELSE STOP 'BACK_INT_GRID: Bad MU' ENDIF C IOCT is the octant the ray come from or the discrete ordinate goes to. IOCT = 1 + BITX + 2*BITY + 4*BITZ IF (BTEST(IPFLAG,1)) THEN JOCT = JOCTORDER2(IOCT) ELSE JOCT = JOCTORDER3(IOCT) ENDIF IZ = (NZ-1)*(1-BITZ)+1 C Sweep through all the grid points DO IORDER = 1, NPTS IPCELL = ISHFT(SWEEPORD(IORDER,JOCT),-3) ICORNER = IBITS(SWEEPORD(IORDER,JOCT),0,3)+1 IPT = GRIDPTR(ICORNER,IPCELL) IF (BTEST(BCFLAG,2) .OR. BTEST(BCFLAG,3)) THEN C For multiple processors get the subdomain boundary radiances C for each Z slab IF ((GRIDPOS(3,IPT) .LT. ZGRID(IZ) .AND. BITZ.EQ.0) .OR. . (GRIDPOS(3,IPT) .GT. ZGRID(IZ) .AND. BITZ.EQ.1)) THEN IZ = IZ + 2*BITZ-1 CALL CALC_BOUNDARY_RADIANCES (BCFLAG, IPFLAG, JOCT, IZ, . NX, NY, NZ, XGRID, YGRID, ZGRID, . NA, NPTS, NCELLS, GRIDPTR, . NEIGHPTR, TREEPTR, CELLFLAGS, GRIDPOS, . MU, PHI, EXTINCT, NSTOKES, SOURCE, . KANG, GRIDRAD) ENDIF ENDIF IF (GRIDRAD(1,IPT) .GE. 0.0) GOTO 500 C Start off at the grid point ICELL = IPCELL TRANSMIT = 1.0 EXT1 = EXTINCT(IPT) RAD(:) = 0.0 SRCEXT1(:) = EXT1*SOURCE(:,KANG,IPT) XE = GRIDPOS(1,IPT) YE = GRIDPOS(2,IPT) ZE = GRIDPOS(3,IPT) C Loop until finding a face with known radiances VALIDRAD = .FALSE. DO WHILE (.NOT. VALIDRAD) C Make sure current cell is valid IF (ICELL .LE. 0) THEN WRITE (6,*) 'BACK_INT_GRID: ICELL=0 ', . MU, PHI, IPT, INEXTCELL,xe,ye,ze, . gridpos(1,ipt),gridpos(2,ipt),gridpos(3,ipt) STOP ENDIF IPINX = BTEST(INT(CELLFLAGS(ICELL)),0) IPINY = BTEST(INT(CELLFLAGS(ICELL)),1) C Find boundaries of the current cell C Find the three possible intersection planes (X,Y,Z) C from the coordinates of the opposite corner grid point IOPP = GRIDPTR(9-IOCT,ICELL) C Get the distances to the 3 planes and select the closest IF (IPINX) THEN SOX = 1.0E20 ELSE SOX = (GRIDPOS(1,IOPP)-XE)*CXINV ENDIF IF (IPINY) THEN SOY = 1.0E20 ELSE SOY = (GRIDPOS(2,IOPP)-YE)*CYINV ENDIF SOZ = (GRIDPOS(3,IOPP)-ZE)*CZINV SO = MIN(SOX,SOY,SOZ) IF (SO .LT. -EPS) THEN WRITE (6,*) 'BACK_INT_GRID: SO<0 ', . MU,PHI,XE,YE,ZE,SOX,SOY,SOZ,SO,IPT,ICELL STOP ENDIF C Compute the coordinates of the cell exitting location XE = XE + SO*CX YE = YE + SO*CY ZE = ZE + SO*CZ C Get the intersection face number (i.e. neighptr index) IF (SOX .LE. SOZ .AND. SOX .LE. SOY) THEN IFACE = 2-BITX JFACE = 1 ELSE IF (SOY .LE. SOZ) THEN IFACE = 4-BITY JFACE = 2 ELSE IFACE = 6-BITZ JFACE = 3 ENDIF C Get the next cell to go to INEXTCELL = NEIGHPTR(IFACE,ICELL) IF (INEXTCELL .LT. 0) THEN CALL NEXT_CELL (XE, YE, ZE, IFACE, JFACE, ICELL, GRIDPOS, . GRIDPTR, NEIGHPTR, TREEPTR, CELLFLAGS, INEXTCELL) ENDIF C Get the face grid pointers IF (NEIGHPTR(IFACE,ICELL) .GE. 0) THEN C If going to same or larger face then use previous face KFACE = IFACE IC = ICELL ELSE C If going to smaller face then use next face (more accurate) KFACE = OPPFACE(IFACE) IC = INEXTCELL ENDIF I1 = GRIDPTR(GRIDFACE(1,KFACE),IC) I2 = GRIDPTR(GRIDFACE(2,KFACE),IC) I3 = GRIDPTR(GRIDFACE(3,KFACE),IC) I4 = GRIDPTR(GRIDFACE(4,KFACE),IC) C Compute the face interpolation factors IF (JFACE .EQ. 1) THEN U = (ZE-GRIDPOS(3,I1))/(GRIDPOS(3,I3)-GRIDPOS(3,I1)) IF (IPINY) THEN V = 0.5 ELSE V = (YE-GRIDPOS(2,I1))/(GRIDPOS(2,I2)-GRIDPOS(2,I1)) ENDIF ELSE IF (JFACE .EQ. 2) THEN U = (ZE-GRIDPOS(3,I1))/(GRIDPOS(3,I3)-GRIDPOS(3,I1)) IF (IPINX) THEN V = 0.5 ELSE V = (XE-GRIDPOS(1,I1))/(GRIDPOS(1,I2)-GRIDPOS(1,I1)) ENDIF ELSE IF (IPINY) THEN U = 0.5 ELSE U = (YE-GRIDPOS(2,I1))/(GRIDPOS(2,I3)-GRIDPOS(2,I1)) ENDIF IF (IPINX) THEN V = 0.5 ELSE V = (XE-GRIDPOS(1,I1))/(GRIDPOS(1,I2)-GRIDPOS(1,I1)) ENDIF ENDIF if (u .lt. -0.0001 .or. v. lt. -0.0001 .or. . u .gt. 1.0001 .or. v .gt. 1.0001) then print *, 'BACK_INT_GRID3D: u,v<0 or u,v>1: ', . mu,phi,xe,ye,ze,u,v,iface,jface,kface, . gridpos(1,i1),gridpos(2,i1),gridpos(3,i1) endif C Get the location coordinate (does the boundary wrapping) IF (INEXTCELL .GT. 0) THEN IF (JFACE .EQ. 1) THEN XE = GRIDPOS(1,GRIDPTR(IOCT,INEXTCELL)) ELSE IF (JFACE .EQ. 2) THEN YE = GRIDPOS(2,GRIDPTR(IOCT,INEXTCELL)) ELSE ZE = GRIDPOS(3,GRIDPTR(IOCT,INEXTCELL)) ENDIF ENDIF C Interpolate extinction and source function at face intersection F1 = (1-U)*(1-V) F2 = (1-U)*V F3 = U*(1-V) F4 = U*V EXT0 = F1*EXTINCT(I1) + F2*EXTINCT(I2) . + F3*EXTINCT(I3) + F4*EXTINCT(I4) C Correctly interpolate source using extinction*source SRCEXT0(:) = F1*SOURCE(:,KANG,I1)*EXTINCT(I1) . + F2*SOURCE(:,KANG,I2)*EXTINCT(I2) . + F3*SOURCE(:,KANG,I3)*EXTINCT(I3) . + F4*SOURCE(:,KANG,I4)*EXTINCT(I4) C Compute the cell radiance: integration of the source function EXT = 0.5*(EXT0+EXT1) TAU=EXT*SO IF (TAU .GE. 0.5) THEN TRANSCELL = EXP(-TAU) ABSCELL = 1.0 - TRANSCELL ELSE ABSCELL = TAU*(1.0-0.5*TAU* . (1.0-0.33333333333*TAU*(1-0.25*TAU))) TRANSCELL = 1.0 - ABSCELL ENDIF IF (TAU .LE. 2.0) THEN IF (EXT .EQ. 0.0) THEN SRC(:) = 0.0 ELSE C Linear extinction, linear source*extinction, to first order SRC(:) = ( 0.5*(SRCEXT0(:)+SRCEXT1(:)) + 0.08333333333 . *(EXT0*SRCEXT1(:)-EXT1*SRCEXT0(:))*SO )/EXT ENDIF ELSE C Combined first order expansion and constant extinction formula EXT0P = EXT0 SRCEXT0P(:) = SRCEXT0(:) IF (TAU .GT. 4.0) THEN EXT0P = EXT1 + (EXT0-EXT1)*4.0/TAU IF (EXT0 .GT. 0.0) SRCEXT0P(:) = SRCEXT0(:)*EXT0P/EXT0 ENDIF SRC(:) = 1.0/(EXT0P+EXT1) *( SRCEXT0P(:)+SRCEXT1(:) . + (EXT0P*SRCEXT1(:)-EXT1*SRCEXT0P(:))*2.0/(EXT0P+EXT1) . *(1-2/TAU+2*TRANSCELL/ABSCELL) ) ENDIF SRC(1) = MAX(SRC(1),0.0D0) C Add in the cell radiance and update the transmission. C If this is a boundary or the transmission is below the C cutoff and we have a valid radiance then set the flag C to stop the tracing and add in the interpolated face radiance. RAD(:) = RAD(:) + TRANSMIT*SRC(:)*ABSCELL TRANSMIT = TRANSMIT*TRANSCELL if (RAD(1) < -1.0E-5) then print '(A,3(1X,F6.4),8(1X,E12.5))', . 'BACK_INT_GRID3D: RAD<0 ', . xe,ye,ze,ext0,ext1,so,tau,transmit,src,abscell,rad(:) CALL ABORT_SHDOM_MPI ('The end') endif VALIDFACE=(GRIDRAD(1,I1).GE.-0.1 .AND.GRIDRAD(1,I2).GE.-0.1 . .AND. GRIDRAD(1,I3).GE.-0.1 .AND.GRIDRAD(1,I4).GE.-0.1) IF (INEXTCELL .LE. 0 .OR. . (TRANSMIT .LE. TRANSMIN .AND. VALIDFACE)) THEN IF (VALIDFACE) THEN VALIDRAD = .TRUE. RAD0(:) = F1*GRIDRAD(:,I1) + F2*GRIDRAD(:,I2) . + F3*GRIDRAD(:,I3) + F4*GRIDRAD(:,I4) RAD(:) = RAD(:) + TRANSMIT*RAD0(:) ELSE print *, 'BACK_INT_GRID3D: INEXTCELL=0: ', xe,ye,ze,icell print *, iz,joct,mu,phi print '(i7,4(1x,f7.4))', i1,GRIDRAD(1,I1), . gridpos(1,i1),gridpos(2,i1),gridpos(3,i1) print '(i7,4(1x,f7.4))', i2,GRIDRAD(1,I2), . gridpos(1,i2),gridpos(2,i2),gridpos(3,i2) print '(i7,4(1x,f7.4))', i3,GRIDRAD(1,I3), . gridpos(1,i3),gridpos(2,i3),gridpos(3,i3) print '(i7,4(1x,f7.4))', i4,GRIDRAD(1,I4), . gridpos(1,i4),gridpos(2,i4),gridpos(3,i4) print '(i7,8x,3(1x,f7.4))', ipt, . gridpos(1,ipt),gridpos(2,ipt),gridpos(3,ipt) STOP ENDIF ELSE EXT1 = EXT0 SRCEXT1(:) = SRCEXT0(:) ICELL = INEXTCELL ENDIF ENDDO GRIDRAD(:,IPT) = RAD(:) 500 CONTINUE ENDDO RETURN END SUBROUTINE BACK_INT_GRID3D_UNPOL (NX,NY,NZ, NA, NPTS, NCELLS, . GRIDPTR, NEIGHPTR, TREEPTR, CELLFLAGS, . GRIDPOS, XGRID, YGRID, ZGRID, SWEEPORD, . MU, PHI, TRANSMIN, BCFLAG, IPFLAG, . EXTINCT, SOURCE, KANG, GRIDRAD) C Sweeps through the spatial grid computing radiances for the C discrete ordinate direction (MU,PHI) by integrating the source C function and extinction backward from each point. The sweeping C order through the grid points is in SWEEPORD. If a point is C already done (nonnegative in GRIDRAD) then it is skipped. C For each grid point a ray opposite to the ordinate direction is traced C back to a face containing known (valid) radiances. As each cell is C crossed (usually only one) the integration of the source function C across the cell is done and added to the radiance. The radiance C at the end of the ray path is interpolated from the known grid point C radiances on the face. The resulting radiances are put in GRIDRAD. INTEGER NX, NY, NZ, NA, NPTS, NCELLS, KANG, BCFLAG, IPFLAG INTEGER GRIDPTR(8,NCELLS), NEIGHPTR(6,NCELLS), TREEPTR(2,NCELLS) INTEGER SWEEPORD(NPTS,*) INTEGER*2 CELLFLAGS(NCELLS) REAL GRIDPOS(3,NPTS), XGRID(*), YGRID(*), ZGRID(NZ) REAL MU, PHI, TRANSMIN REAL EXTINCT(NPTS), SOURCE(NA,NPTS), GRIDRAD(NPTS) INTEGER BITX, BITY, BITZ, IOCT, JOCT INTEGER IPCELL, ICELL, INEXTCELL, IPT, IFACE INTEGER I1, I2, I3, I4, IOPP, JFACE, KFACE, IC, IZ INTEGER GRIDFACE(4,6), OPPFACE(6) INTEGER IORDER, ICORNER, JOCTORDER3(8), JOCTORDER2(8) LOGICAL VALIDRAD, VALIDFACE, IPINX, IPINY, BTEST DOUBLE PRECISION PI, CX, CY, CZ, CXINV, CYINV, CZINV DOUBLE PRECISION XE, YE, ZE DOUBLE PRECISION SO, SOX, SOY, SOZ, EPS, path DOUBLE PRECISION U, V, F1, F2, F3, F4 DOUBLE PRECISION EXT, EXT0, EXT1, SRC, SRCEXT0, SRCEXT1 DOUBLE PRECISION EXT0P, SRCEXT0P DOUBLE PRECISION TAU, TRANSCELL, ABSCELL, RAD, RAD0, TRANSMIT DATA GRIDFACE/1,3,5,7, 2,4,6,8, 1,2,5,6, 3,4,7,8, . 1,2,3,4, 5,6,7,8/, OPPFACE/2,1,4,3,6,5/ DATA JOCTORDER3/1,3,5,7,2,4,6,8/, JOCTORDER2/1,3,1,3,2,4,2,4/ EPS = 1.0E-3*(GRIDPOS(3,GRIDPTR(8,1))-GRIDPOS(3,GRIDPTR(1,1))) C Make the ray direction (opposite to the discrete ordinate direction) PI = ACOS(-1.0) CX = SQRT(1.0-MU**2)*COS(PHI+PI) CY = SQRT(1.0-MU**2)*SIN(PHI+PI) CZ = -MU IF (ABS(CX) .GT. 1.0E-5) THEN CXINV = 1.0D0/CX ELSE CX = 0.0 CXINV = 1.0E6 ENDIF IF (ABS(CY) .GT. 1.0E-5) THEN CYINV = 1.0D0/CY ELSE CY = 0.0 CYINV = 1.0E6 ENDIF CZINV = 1.0D0/CZ C Setup the sweeping direction and the gridpoint location C BITc is 0 for positive direction trace back ray, 1 for negative. IF (CX .LT. 0.0) THEN BITX = 1 ELSE BITX = 0 ENDIF IF (CY .LT. 0.0) THEN BITY = 1 ELSE BITY = 0 ENDIF IF (CZ .LT. -1.0E-3) THEN BITZ = 1 ELSE IF (CZ .GT. 1.0E-3) THEN BITZ = 0 ELSE STOP 'BACK_INT_GRID: Bad MU' ENDIF C IOCT is the octant the ray come from or the discrete ordinate goes to. IOCT = 1 + BITX + 2*BITY + 4*BITZ IF (BTEST(IPFLAG,1)) THEN JOCT = JOCTORDER2(IOCT) ELSE JOCT = JOCTORDER3(IOCT) ENDIF IZ = (NZ-1)*(1-BITZ)+1 C Sweep through all the grid points DO IORDER = 1, NPTS IPCELL = ISHFT(SWEEPORD(IORDER,JOCT),-3) ICORNER = IBITS(SWEEPORD(IORDER,JOCT),0,3)+1 IPT = GRIDPTR(ICORNER,IPCELL) IF (BTEST(BCFLAG,2) .OR. BTEST(BCFLAG,3)) THEN C For multiple processors get the subdomain boundary radiances C for each Z slab IF ((GRIDPOS(3,IPT) .LT. ZGRID(IZ) .AND. BITZ.EQ.0) .OR. . (GRIDPOS(3,IPT) .GT. ZGRID(IZ) .AND. BITZ.EQ.1)) THEN IZ = IZ + 2*BITZ-1 CALL CALC_BOUNDARY_RADIANCES (BCFLAG, IPFLAG, JOCT, IZ, . NX, NY, NZ, XGRID, YGRID, ZGRID, . NA, NPTS, NCELLS, GRIDPTR, . NEIGHPTR, TREEPTR, CELLFLAGS, GRIDPOS, . MU, PHI, EXTINCT, 1, SOURCE, . KANG, GRIDRAD) ENDIF ENDIF IF (GRIDRAD(IPT) .GE. 0.0) GOTO 500 C Start off at the grid point ICELL = IPCELL RAD = 0.0 TRANSMIT = 1.0 EXT1 = EXTINCT(IPT) SRCEXT1 = EXT1*SOURCE(KANG,IPT) XE = GRIDPOS(1,IPT) YE = GRIDPOS(2,IPT) ZE = GRIDPOS(3,IPT) C Loop until finding a face with known radiances VALIDRAD = .FALSE. DO WHILE (.NOT. VALIDRAD) C Make sure current cell is valid IF (ICELL .LE. 0) THEN WRITE (6,*) 'BACK_INT_GRID: ICELL=0 ', . MU, PHI, IPT, INEXTCELL,xe,ye,ze, . gridpos(1,ipt),gridpos(2,ipt),gridpos(3,ipt) STOP ENDIF IPINX = BTEST(INT(CELLFLAGS(ICELL)),0) IPINY = BTEST(INT(CELLFLAGS(ICELL)),1) C Find boundaries of the current cell C Find the three possible intersection planes (X,Y,Z) C from the coordinates of the opposite corner grid point IOPP = GRIDPTR(9-IOCT,ICELL) C Get the distances to the 3 planes and select the closest IF (IPINX) THEN SOX = 1.0E20 ELSE SOX = (GRIDPOS(1,IOPP)-XE)*CXINV ENDIF IF (IPINY) THEN SOY = 1.0E20 ELSE SOY = (GRIDPOS(2,IOPP)-YE)*CYINV ENDIF SOZ = (GRIDPOS(3,IOPP)-ZE)*CZINV SO = MIN(SOX,SOY,SOZ) IF (SO .LT. -EPS) THEN WRITE (6,*) 'BACK_INT_GRID: SO<0 ', . MU,PHI,XE,YE,ZE,SOX,SOY,SOZ,SO,IPT,ICELL STOP ENDIF C Compute the coordinates of the cell exitting location XE = XE + SO*CX YE = YE + SO*CY ZE = ZE + SO*CZ C Get the intersection face number (i.e. neighptr index) IF (SOX .LE. SOZ .AND. SOX .LE. SOY) THEN IFACE = 2-BITX JFACE = 1 ELSE IF (SOY .LE. SOZ) THEN IFACE = 4-BITY JFACE = 2 ELSE IFACE = 6-BITZ JFACE = 3 ENDIF C Get the next cell to go to INEXTCELL = NEIGHPTR(IFACE,ICELL) IF (INEXTCELL .LT. 0) THEN CALL NEXT_CELL (XE, YE, ZE, IFACE, JFACE, ICELL, GRIDPOS, . GRIDPTR, NEIGHPTR, TREEPTR, CELLFLAGS, INEXTCELL) ENDIF C Get the face grid pointers IF (NEIGHPTR(IFACE,ICELL) .GE. 0) THEN C If going to same or larger face then use previous face KFACE = IFACE IC = ICELL ELSE C If going to smaller face then use next face (more accurate) KFACE = OPPFACE(IFACE) IC = INEXTCELL ENDIF I1 = GRIDPTR(GRIDFACE(1,KFACE),IC) I2 = GRIDPTR(GRIDFACE(2,KFACE),IC) I3 = GRIDPTR(GRIDFACE(3,KFACE),IC) I4 = GRIDPTR(GRIDFACE(4,KFACE),IC) C Compute the face interpolation factors IF (JFACE .EQ. 1) THEN U = (ZE-GRIDPOS(3,I1))/(GRIDPOS(3,I3)-GRIDPOS(3,I1)) IF (IPINY) THEN V = 0.5 ELSE V = (YE-GRIDPOS(2,I1))/(GRIDPOS(2,I2)-GRIDPOS(2,I1)) ENDIF ELSE IF (JFACE .EQ. 2) THEN U = (ZE-GRIDPOS(3,I1))/(GRIDPOS(3,I3)-GRIDPOS(3,I1)) IF (IPINX) THEN V = 0.5 ELSE V = (XE-GRIDPOS(1,I1))/(GRIDPOS(1,I2)-GRIDPOS(1,I1)) ENDIF ELSE IF (IPINY) THEN U = 0.5 ELSE U = (YE-GRIDPOS(2,I1))/(GRIDPOS(2,I3)-GRIDPOS(2,I1)) ENDIF IF (IPINX) THEN V = 0.5 ELSE V = (XE-GRIDPOS(1,I1))/(GRIDPOS(1,I2)-GRIDPOS(1,I1)) ENDIF ENDIF if (u .lt. -0.0001 .or. v. lt. -0.0001 .or. . u .gt. 1.0001 .or. v .gt. 1.0001) then print *, 'BACK_INT_GRID3D: u,v<0 or u,v>1: ', . mu,phi,xe,ye,ze,u,v,iface,jface,kface, . gridpos(1,i1),gridpos(2,i1),gridpos(3,i1) endif C Get the location coordinate (does the boundary wrapping) IF (INEXTCELL .GT. 0) THEN IF (JFACE .EQ. 1) THEN XE = GRIDPOS(1,GRIDPTR(IOCT,INEXTCELL)) ELSE IF (JFACE .EQ. 2) THEN YE = GRIDPOS(2,GRIDPTR(IOCT,INEXTCELL)) ELSE ZE = GRIDPOS(3,GRIDPTR(IOCT,INEXTCELL)) ENDIF ENDIF C Interpolate extinction and source function at face intersection F1 = (1-U)*(1-V) F2 = (1-U)*V F3 = U*(1-V) F4 = U*V EXT0 = F1*EXTINCT(I1) + F2*EXTINCT(I2) . + F3*EXTINCT(I3) + F4*EXTINCT(I4) C Correctly interpolate source using extinction*source SRCEXT0 = (F1*SOURCE(KANG,I1)*EXTINCT(I1) . + F2*SOURCE(KANG,I2)*EXTINCT(I2) . + F3*SOURCE(KANG,I3)*EXTINCT(I3) . + F4*SOURCE(KANG,I4)*EXTINCT(I4)) C Compute the cell radiance: integration of the source function EXT = 0.5*(EXT0+EXT1) TAU=EXT*SO IF (TAU .GE. 0.5) THEN TRANSCELL = EXP(-TAU) ABSCELL = 1.0 - TRANSCELL ELSE ABSCELL = TAU*(1.0-0.5*TAU* . (1.0-0.33333333333*TAU*(1-0.25*TAU))) TRANSCELL = 1.0 - ABSCELL ENDIF IF (TAU .LE. 2.0) THEN IF (EXT .EQ. 0.0) THEN SRC = 0.0 ELSE C Linear extinction, linear source*extinction, to first order SRC = ( 0.5*(SRCEXT0+SRCEXT1) . + 0.08333333333*(EXT0*SRCEXT1-EXT1*SRCEXT0)*SO )/EXT ENDIF ELSE C Combined first order expansion and constant extinction formula c SRC = 0.5/EXT *( SRCEXT0+SRCEXT1 c . + (EXT0*SRCEXT1-EXT1*SRCEXT0)*SO c . *(1-2/TAU+2*TRANSCELL/ABSCELL)/TAU ) EXT0P = EXT0 SRCEXT0P = SRCEXT0 IF (TAU .GT. 4.0) THEN EXT0P = EXT1 + (EXT0-EXT1)*4.0/TAU IF (EXT0 .GT. 0.0) SRCEXT0P = SRCEXT0*EXT0P/EXT0 ENDIF SRC = 1.0/(EXT0P+EXT1) *( SRCEXT0P+SRCEXT1 . + (EXT0P*SRCEXT1-EXT1*SRCEXT0P)*2.0/(EXT0P+EXT1) . *(1-2/TAU+2*TRANSCELL/ABSCELL) ) ENDIF SRC = MAX(SRC,0.0D0) C Add in the cell radiance and update the transmission. C If this is a boundary or the transmission is below the C cutoff and we have a valid radiance then set the flag C to stop the tracing and add in the interpolated face radiance. RAD = RAD + TRANSMIT*SRC*ABSCELL TRANSMIT = TRANSMIT*TRANSCELL if (RAD < -1.0E-5) then print '(A,3(1X,F6.4),8(1X,E12.5))', . 'BACK_INT_GRID3D: RAD<0 ', . xe,ye,ze,ext0,ext1,so,tau,transmit,src,abscell,rad CALL ABORT_SHDOM_MPI ('The end') endif VALIDFACE = (GRIDRAD(I1).GE.-0.1 .AND. GRIDRAD(I2).GE.-0.1 . .AND. GRIDRAD(I3).GE.-0.1 .AND. GRIDRAD(I4).GE.-0.1) IF (INEXTCELL .LE. 0 .OR. . (TRANSMIT .LE. TRANSMIN .AND. VALIDFACE)) THEN IF (VALIDFACE) THEN VALIDRAD = .TRUE. RAD0 = F1*GRIDRAD(I1) + F2*GRIDRAD(I2) . + F3*GRIDRAD(I3) + F4*GRIDRAD(I4) RAD = RAD + TRANSMIT*RAD0 ELSE print *, 'BACK_INT_GRID3D: INEXTCELL=0: ', xe,ye,ze,icell print *, iz,joct,mu,phi print '(i7,4(1x,f7.4))', i1,GRIDRAD(I1), . gridpos(1,i1),gridpos(2,i1),gridpos(3,i1) print '(i7,4(1x,f7.4))', i2,GRIDRAD(I2), . gridpos(1,i2),gridpos(2,i2),gridpos(3,i2) print '(i7,4(1x,f7.4))', i3,GRIDRAD(I3), . gridpos(1,i3),gridpos(2,i3),gridpos(3,i3) print '(i7,4(1x,f7.4))', i4,GRIDRAD(I4), . gridpos(1,i4),gridpos(2,i4),gridpos(3,i4) print '(i7,8x,3(1x,f7.4))', ipt, . gridpos(1,ipt),gridpos(2,ipt),gridpos(3,ipt) STOP ENDIF ELSE EXT1 = EXT0 SRCEXT1 = SRCEXT0 ICELL = INEXTCELL ENDIF ENDDO GRIDRAD(IPT) = RAD 500 CONTINUE ENDDO RETURN END SUBROUTINE BACK_INT_GRID2D (NX,NY,NZ, NA, NPTS, NCELLS, . GRIDPTR, NEIGHPTR, TREEPTR, CELLFLAGS, . GRIDPOS, SWEEPORD, MU, PHI, TRANSMIN, . EXTINCT, NSTOKES, SOURCE, KANG, GRIDRAD) C Sweeps through the spatial grid computing radiances for the C discrete ordinate direction (MU,PHI) by integrating the source C function and extinction backward from each point. If a point is C already done (nonnegative in GRIDRAD) then it is skipped. C For each grid point a ray opposite to the ordinate direction is traced C back to a face containing known (valid) radiances. As each cell is C crossed (usually only one) the integration of the source function C across the cell is done and added to the radiance. The radiance C ad the end of the ray path is interpolated from the known grid point C radiances on the face. The resulting radiances are put in GRIDRAD. IMPLICIT NONE INTEGER NX, NY, NZ, NA, NPTS, NCELLS, KANG, NSTOKES INTEGER GRIDPTR(8,NCELLS), NEIGHPTR(6,NCELLS), TREEPTR(2,NCELLS) INTEGER SWEEPORD(NPTS,*) INTEGER*2 CELLFLAGS(NCELLS) REAL GRIDPOS(3,NPTS) REAL MU, PHI, TRANSMIN REAL EXTINCT(NPTS), SOURCE(NSTOKES,NA,NPTS) REAL GRIDRAD(NSTOKES,NPTS) INTEGER BITX, BITZ, IOCT, JOCT INTEGER IPCELL, ICELL, INEXTCELL, IPT, IFACE INTEGER I1, I2, IOPP, JFACE, KFACE, IC, IZ INTEGER GRIDFACE(2,6), OPPFACE(6) INTEGER IORDER, ICORNER, JOCTORDER(8) LOGICAL VALIDRAD, VALIDFACE, IPINX, BTEST DOUBLE PRECISION PI, CX, CZ, CXINV, CZINV DOUBLE PRECISION XE, YE, ZE DOUBLE PRECISION SO, SOX, SOZ, EPS DOUBLE PRECISION U, F1, F2 DOUBLE PRECISION EXT, EXT0, EXT1 DOUBLE PRECISION SRC(NSTOKES), SRCEXT0(NSTOKES), SRCEXT1(NSTOKES) DOUBLE PRECISION EXT0P, SRCEXT0P(NSTOKES) DOUBLE PRECISION TAU, TRANSCELL, ABSCELL, TRANSMIT DOUBLE PRECISION RAD(NSTOKES), RAD0(NSTOKES) DATA GRIDFACE/1,5, 2,6, 0,0, 0,0, 1,2, 5,6/ DATA OPPFACE/2,1,4,3,6,5/ DATA JOCTORDER/1,3,1,3,2,4,2,4/ EPS = 1.0E-3*(GRIDPOS(3,GRIDPTR(8,1))-GRIDPOS(3,GRIDPTR(1,1))) C Make the ray direction (opposite to the discrete ordinate direction) PI = ACOS(-1.0) CX = SQRT(1.0-MU**2)*COS(PHI+PI) CZ = -MU IF (ABS(CX) .GT. 1.0E-5) THEN CXINV = 1.0D0/CX ELSE CX = 0.0 CXINV = 1.0E6 ENDIF CZINV = 1.0D0/CZ C Setup the sweeping direction and the gridpoint location C BITc is 0 for positive direction trace back ray, 1 for negative. IF (CX .LT. 0.0) THEN BITX = 1 ELSE BITX = 0 ENDIF IF (CZ .LT. -1.0E-3) THEN BITZ = 1 ELSE IF (CZ .GT. 1.0E-3) THEN BITZ = 0 ELSE STOP 'BACK_INT_GRID: Bad MU' ENDIF C IOCT is the octant the ray come from or the discrete ordinate goes to. IOCT = 1 + BITX + 4*BITZ JOCT = JOCTORDER(IOCT) C Sweep through all the grid points DO IORDER = 1, NPTS IPCELL = ISHFT(SWEEPORD(IORDER,JOCT),-3) ICORNER = IBITS(SWEEPORD(IORDER,JOCT),0,3)+1 IPT = GRIDPTR(ICORNER,IPCELL) IF (GRIDRAD(1,IPT) .GE. 0.0) GOTO 500 C Start off at the grid point ICELL = IPCELL TRANSMIT = 1.0 EXT1 = EXTINCT(IPT) RAD(:) = 0.0 SRCEXT1(:) = EXT1*SOURCE(:,KANG,IPT) XE = GRIDPOS(1,IPT) YE = GRIDPOS(2,IPT) ZE = GRIDPOS(3,IPT) C Loop until finding a face with known radiances VALIDRAD = .FALSE. DO WHILE (.NOT. VALIDRAD) C Make sure current cell is valid IF (ICELL .LE. 0) THEN WRITE (6,*) 'BACK_INT_GRID: ICELL=0 ', . MU, PHI, IPT, INEXTCELL STOP ENDIF IPINX = BTEST(INT(CELLFLAGS(ICELL)),0) C Find boundaries of the current cell C Find the three possible intersection planes (X,Y,Z) C from the coordinates of the opposite corner grid point IOPP = GRIDPTR(9-IOCT,ICELL) C Get the distances to the 3 planes and select the closest IF (IPINX) THEN SOX = 1.0E20 ELSE SOX = (GRIDPOS(1,IOPP)-XE)*CXINV ENDIF SOZ = (GRIDPOS(3,IOPP)-ZE)*CZINV SO = MIN(SOX,SOZ) IF (SO .LT. -EPS) THEN WRITE (6,*) 'BACK_INT_GRID2D: SO<0 ', . MU,PHI,XE,ZE,SOX,SOZ,SO,IPT,ICELL STOP ENDIF C Compute the coordinates of the cell exitting location XE = XE + SO*CX ZE = ZE + SO*CZ C Get the intersection face number (i.e. neighptr index) IF (SOX .LE. SOZ) THEN IFACE = 2-BITX JFACE = 1 ELSE IFACE = 6-BITZ JFACE = 3 ENDIF C Get the next cell to go to INEXTCELL = NEIGHPTR(IFACE,ICELL) IF (INEXTCELL .LT. 0) THEN CALL NEXT_CELL (XE, YE, ZE, IFACE, JFACE, ICELL, GRIDPOS, . GRIDPTR, NEIGHPTR, TREEPTR, CELLFLAGS, INEXTCELL) ENDIF C Get the face grid pointers IF (NEIGHPTR(IFACE,ICELL) .GE. 0) THEN C If going to same or larger face then use previous face KFACE = IFACE IC = ICELL ELSE C If going to smaller face then use next face (more accurate) KFACE = OPPFACE(IFACE) IC = INEXTCELL ENDIF I1 = GRIDPTR(GRIDFACE(1,KFACE),IC) I2 = GRIDPTR(GRIDFACE(2,KFACE),IC) C Compute the face interpolation factors IF (JFACE .EQ. 1) THEN U = (ZE-GRIDPOS(3,I1))/(GRIDPOS(3,I2)-GRIDPOS(3,I1)) ELSE IF (IPINX) THEN U = 0.5 ELSE U = (XE-GRIDPOS(1,I1))/(GRIDPOS(1,I2)-GRIDPOS(1,I1)) ENDIF ENDIF C Get the location coordinate (does the boundary wrapping) IF (INEXTCELL .GT. 0) THEN IF (JFACE .EQ. 1) THEN XE = GRIDPOS(1,GRIDPTR(IOCT,INEXTCELL)) ELSE ZE = GRIDPOS(3,GRIDPTR(IOCT,INEXTCELL)) ENDIF ENDIF C Interpolate extinction and source function at face intersection F1 = 1-U F2 = U EXT0 = F1*EXTINCT(I1) + F2*EXTINCT(I2) C Correctly interpolate source using extinction*source SRCEXT0(:) = (F1*SOURCE(:,KANG,I1)*EXTINCT(I1) . + F2*SOURCE(:,KANG,I2)*EXTINCT(I2)) C Compute the cell radiance: integration of the source function EXT = 0.5*(EXT0+EXT1) TAU=EXT*SO IF (TAU .GE. 0.5) THEN TRANSCELL = EXP(-TAU) ABSCELL = 1.0 - TRANSCELL ELSE ABSCELL = TAU*(1.0-0.5*TAU* . (1.0-0.33333333333*TAU*(1-0.25*TAU))) TRANSCELL = 1.0 - ABSCELL ENDIF IF (TAU .LE. 2.0) THEN IF (EXT .EQ. 0.0) THEN SRC(:) = 0.0 ELSE C Linear extinction, linear source*extinction, to first order SRC(:) = ( 0.5*(SRCEXT0(:)+SRCEXT1(:)) + 0.08333333333 . *(EXT0*SRCEXT1(:)-EXT1*SRCEXT0(:))*SO )/EXT ENDIF ELSE C Combined first order expansion and constant extinction formula EXT0P = EXT0 SRCEXT0P(:) = SRCEXT0(:) IF (TAU .GT. 4.0) THEN EXT0P = EXT1 + (EXT0-EXT1)*4.0/TAU IF (EXT0 .GT. 0.0) SRCEXT0P(:) = SRCEXT0(:)*EXT0P/EXT0 ENDIF SRC(:) = 1.0/(EXT0P+EXT1) *( SRCEXT0P(:)+SRCEXT1(:) . + (EXT0P*SRCEXT1(:)-EXT1*SRCEXT0P(:))*2.0/(EXT0P+EXT1) . *(1-2/TAU+2*TRANSCELL/ABSCELL) ) ENDIF SRC(1) = MAX(SRC(1),0.0D0) C Add in the cell radiance and update the transmission. C If this is a boundary or the transmission is below the C cutoff and we have a valid radiance then set the flag C to stop the tracing and add in the interpolated face radiance. RAD(:) = RAD(:) + TRANSMIT*SRC(:)*ABSCELL TRANSMIT = TRANSMIT*TRANSCELL VALIDFACE=(GRIDRAD(1,I1) .GE.-0.1 .AND.GRIDRAD(1,I2).GE.-0.1) IF (INEXTCELL .LE. 0 .OR. . (TRANSMIT .LE. TRANSMIN .AND. VALIDFACE)) THEN IF (VALIDFACE) THEN VALIDRAD = .TRUE. RAD0(:) = F1*GRIDRAD(:,I1) + F2*GRIDRAD(:,I2) RAD(:) = RAD(:) + TRANSMIT*RAD0(:) ELSE print *, 'BACK_INT_GRID2D: INEXTCELL=0: ', xe,ye,ze,icell print *, iz,joct,mu,phi print '(i7,4(1x,f7.4))', i1,GRIDRAD(1,I1), . gridpos(1,i1),gridpos(2,i1),gridpos(3,i1) print '(i7,4(1x,f7.4))', i2,GRIDRAD(1,I2), . gridpos(1,i2),gridpos(2,i2),gridpos(3,i2) print '(i7,8x,3(1x,f7.4))', ipt, . gridpos(1,ipt),gridpos(2,ipt),gridpos(3,ipt) STOP ENDIF ELSE EXT1 = EXT0 SRCEXT1(:) = SRCEXT0(:) ICELL = INEXTCELL ENDIF ENDDO GRIDRAD(:,IPT) = RAD(:) 500 CONTINUE ENDDO RETURN END SUBROUTINE BACK_INT_GRID1D (NX,NY,NZ, NA, NPTS, NCELLS, . GRIDPTR, NEIGHPTR, TREEPTR, CELLFLAGS, . GRIDPOS, SWEEPORD, MU, PHI, TRANSMIN, . EXTINCT, NSTOKES, SOURCE, KANG, GRIDRAD) C Sweeps through the spatial grid computing radiances for the C discrete ordinate direction (MU,PHI) by integrating the source C function and extinction backward from each point. If a point is C already done (nonnegative in GRIDRAD) then it is skipped. C For each grid point a ray opposite to the ordinate direction is traced C back to a face containing known (valid) radiances. As each cell is C crossed (usually only one) the integration of the source function C across the cell is done and added to the radiance. The radiance C ad the end of the ray path is interpolated from the known grid point C radiances on the face. The resulting radiances are put in GRIDRAD. IMPLICIT NONE INTEGER NX, NY, NZ, NA, NPTS, NCELLS, KANG, NSTOKES INTEGER GRIDPTR(8,NCELLS), NEIGHPTR(6,NCELLS), TREEPTR(2,NCELLS) INTEGER SWEEPORD(NPTS,*) INTEGER*2 CELLFLAGS(NCELLS) REAL GRIDPOS(3,NPTS) REAL MU, PHI, TRANSMIN REAL EXTINCT(NPTS), SOURCE(NSTOKES,NA,NPTS) REAL GRIDRAD(NSTOKES,NPTS) INTEGER BITZ, IOCT, JOCT INTEGER IPCELL, ICELL, INEXTCELL, IPT, IFACE, I1, IOPP INTEGER GRIDFACE(6) INTEGER IORDER, ICORNER, JOCTORDER(8) LOGICAL VALIDRAD DOUBLE PRECISION PI, CZ, CZINV DOUBLE PRECISION XE, YE, ZE DOUBLE PRECISION SO, EPS DOUBLE PRECISION EXT, EXT0, EXT1 DOUBLE PRECISION SRC(NSTOKES), SRCEXT0(NSTOKES), SRCEXT1(NSTOKES) DOUBLE PRECISION EXT0P, SRCEXT0P(NSTOKES) DOUBLE PRECISION TAU, TRANSCELL, ABSCELL, TRANSMIT DOUBLE PRECISION RAD(NSTOKES), RAD0(NSTOKES) DATA GRIDFACE/0, 0, 0, 0, 1, 5/ DATA JOCTORDER/1,1,1,1,2,2,2,2/ EPS = 1.0E-3*(GRIDPOS(3,GRIDPTR(8,1))-GRIDPOS(3,GRIDPTR(1,1))) C Make the ray direction (opposite to the discrete ordinate direction) PI = ACOS(-1.0) CZ = -MU CZINV = 1.0D0/CZ C Setup the sweeping direction and the gridpoint location C BITc is 0 for positive direction trace back ray, 1 for negative. IF (CZ .LT. -1.0E-3) THEN BITZ = 1 ELSE IF (CZ .GT. 1.0E-3) THEN BITZ = 0 ELSE STOP 'BACK_INT_GRID3D: Bad MU' ENDIF C IOCT is the octant the ray come from or the discrete ordinate goes to. IOCT = 1 + 4*BITZ JOCT = JOCTORDER(IOCT) C Sweep through all the grid points DO IORDER = 1, NPTS IPCELL = ISHFT(SWEEPORD(IORDER,JOCT),-3) ICORNER = IBITS(SWEEPORD(IORDER,JOCT),0,3)+1 IPT = GRIDPTR(ICORNER,IPCELL) IF (GRIDRAD(1,IPT) .GE. 0.0) GOTO 500 C Start off at the grid point ICELL = IPCELL TRANSMIT = 1.0 EXT1 = EXTINCT(IPT) RAD(:) = 0.0 SRCEXT1(:) = EXT1*SOURCE(:,KANG,IPT) XE = GRIDPOS(1,IPT) YE = GRIDPOS(2,IPT) ZE = GRIDPOS(3,IPT) C Loop until finding a face with known radiances VALIDRAD = .FALSE. DO WHILE (.NOT. VALIDRAD) C Make sure current cell is valid IF (ICELL .LE. 0) THEN WRITE (6,*) 'BACK_INT_GRID: ICELL=0 ', . MU, PHI, IPT, INEXTCELL STOP ENDIF C Find boundaries of the current cell C Find the three possible intersection planes (X,Y,Z) C from the coordinates of the opposite corner grid point IOPP = GRIDPTR(9-IOCT,ICELL) C Get the distances to the 3 planes and select the closest SO = (GRIDPOS(3,IOPP)-ZE)*CZINV IF (SO .LT. -EPS) THEN WRITE (6,*) 'BACK_INT_GRID1D: SO<0 ', . MU,PHI,ZE,SO,IPT,ICELL STOP ENDIF C Compute the coordinates of the cell exitting location ZE = ZE + SO*CZ C Get the intersection face number (i.e. neighptr index) IFACE = 6-BITZ C Get the next cell to go to INEXTCELL = NEIGHPTR(IFACE,ICELL) C Get the face grid pointers I1 = GRIDPTR(GRIDFACE(IFACE),ICELL) C Get the location coordinate (does the boundary wrapping) IF (INEXTCELL .GT. 0) THEN ZE = GRIDPOS(3,GRIDPTR(IOCT,INEXTCELL)) ENDIF C Get extinction and source function at face intersection EXT0 = EXTINCT(I1) SRCEXT0(:) = EXTINCT(I1)*SOURCE(:,KANG,I1) C Compute the cell radiance: integration of the source function EXT = 0.5*(EXT0+EXT1) TAU=EXT*SO IF (TAU .GE. 0.5) THEN TRANSCELL = EXP(-TAU) ABSCELL = 1.0 - TRANSCELL ELSE ABSCELL = TAU*(1.0-0.5*TAU* . (1.0-0.33333333333*TAU*(1-0.25*TAU))) TRANSCELL = 1.0 - ABSCELL ENDIF IF (TAU .LE. 2.0) THEN IF (EXT .EQ. 0.0) THEN SRC(:) = 0.0 ELSE C Linear extinction, linear source*extinction, to first order SRC(:) = ( 0.5*(SRCEXT0(:)+SRCEXT1(:)) + 0.08333333333 . *(EXT0*SRCEXT1(:)-EXT1*SRCEXT0(:))*SO )/EXT ENDIF ELSE C Combined first order expansion and constant extinction formula EXT0P = EXT0 SRCEXT0P(:) = SRCEXT0(:) IF (TAU .GT. 4.0) THEN EXT0P = EXT1 + (EXT0-EXT1)*4.0/TAU IF (EXT0 .GT. 0.0) SRCEXT0P(:) = SRCEXT0(:)*EXT0P/EXT0 ENDIF SRC(:) = 1.0/(EXT0P+EXT1) *( SRCEXT0P(:)+SRCEXT1(:) . + (EXT0P*SRCEXT1(:)-EXT1*SRCEXT0P(:))*2.0/(EXT0P+EXT1) . *(1-2/TAU+2*TRANSCELL/ABSCELL) ) ENDIF SRC(1) = MAX(SRC(1),0.0D0) C See if there are valid radiances for the intersection face: C If so, set flag, interpolate radiance, and add to C radiance from source function in cell C If not, only add in cell radiance and prepare for next cell IF (GRIDRAD(1,I1) .GE. -0.1) THEN VALIDRAD = .TRUE. RAD0(:) = GRIDRAD(:,I1) RAD(:) = RAD(:) +TRANSMIT*(RAD0(:)*TRANSCELL+SRC(:)*ABSCELL) ELSE RAD(:) = RAD(:) + TRANSMIT*SRC(:)*ABSCELL SRCEXT1(:) = SRCEXT0(:) TRANSMIT = TRANSMIT*TRANSCELL EXT1 = EXT0 ICELL = INEXTCELL ENDIF ENDDO GRIDRAD(:,IPT) = RAD(:) 500 CONTINUE ENDDO RETURN END SUBROUTINE NEXT_CELL (XE, YE, ZE, IFACE, JFACE, ICELL, GRIDPOS, . GRIDPTR, NEIGHPTR, TREEPTR, CELLFLAGS, INEXT) C Gets the next cell (INEXT) to go to that is adjacent to this face. C IFACE (1-6) has which face (-X,+X,-Y,+Y,-Z,+Z) of the current cell C (ICELL) we are at, and JFACE (1-3) is the direction (X,Y,Z). C The current location (XE,YE,ZE) on the face is used to find the cell C if there are several adjacent to the current cell. IMPLICIT NONE INTEGER IFACE, JFACE, ICELL, INEXT INTEGER GRIDPTR(8,*), NEIGHPTR(6,*), TREEPTR(2,*) INTEGER*2 CELLFLAGS(*) REAL GRIDPOS(3,*) DOUBLE PRECISION XE, YE, ZE INTEGER IC, IC1, DIR C Get the pointer to the neighbor cell INEXT = NEIGHPTR(IFACE,ICELL) IF (INEXT .LT. 0) THEN C If neighbor pointer is negative then we are going into a cell C that has subcells, and the absolute value of the pointer is the C cell to start the tree search. IC = ABS(INEXT) C Go down the tree until there is no children DO WHILE (TREEPTR(2,IC) .GT. 0) C Get the direction cell is split DIR = IBITS(INT(CELLFLAGS(IC)),2,2) IF (DIR .EQ. 0) STOP 'NEXT_CELL: No split direction' C If split is in direction of face then choose the near child cell C which is the negative side cell if this is a positive face. IC1 = TREEPTR(2,IC) IF (DIR .EQ. JFACE) THEN IC = IC1 + 1 - MOD(IFACE-1,2) ELSE C If split is not in direction of face then find which of two cells C has a face that contains the current location. C Compare the appropriate coordinate (X,Y,Z) of the grid point C on the split line with the current location. IC = IC1 IF (DIR .EQ. 1) THEN IF (XE .GT. GRIDPOS(1,GRIDPTR(8,IC1))) IC = IC + 1 ELSE IF (DIR .EQ. 2) THEN IF (YE .GT. GRIDPOS(2,GRIDPTR(8,IC1))) IC = IC + 1 ELSE IF (ZE .GT. GRIDPOS(3,GRIDPTR(8,IC1))) IC = IC + 1 ENDIF ENDIF ENDDO INEXT = IC ENDIF RETURN END