init.f90

FDS为火灾动力学模拟软件源代码,该软件为开源项目,代码语言主要为FORTRAN,可在WINDOWS和LINUX下编译运行,详细说明可参考http://fire.nist.gov/fds/官方网址

N_EDGES = 0
 
! Arguments for DEFINE_EDGE(I,J,K,IOR,IEC,NM,I_OBST)
 
DO K=0,KBAR
   DO J=0,JBAR
      IF (J>0) CALL DEFINE_EDGE(   0,J,K, 1,2,NM,0)
      IF (J>0) CALL DEFINE_EDGE(IBAR,J,K,-1,2,NM,0)
      IF (K>0) CALL DEFINE_EDGE(   0,J,K, 1,3,NM,0)
      IF (K>0) CALL DEFINE_EDGE(IBAR,J,K,-1,3,NM,0)
   ENDDO
ENDDO
DO K=0,KBAR
   DO I=0,IBAR
      IF (I>0) CALL DEFINE_EDGE(I,   0,K, 2,1,NM,0)
      IF (I>0) CALL DEFINE_EDGE(I,JBAR,K,-2,1,NM,0)
      IF (K>0) CALL DEFINE_EDGE(I,   0,K, 2,3,NM,0)
      IF (K>0) CALL DEFINE_EDGE(I,JBAR,K,-2,3,NM,0)
   ENDDO
ENDDO
DO J=0,JBAR
   DO I=0,IBAR
      IF (I>0) CALL DEFINE_EDGE(I,J,   0, 3,1,NM,0)
      IF (I>0) CALL DEFINE_EDGE(I,J,KBAR,-3,1,NM,0)
      IF (J>0) CALL DEFINE_EDGE(I,J,   0, 3,2,NM,0)
      IF (J>0) CALL DEFINE_EDGE(I,J,KBAR,-3,2,NM,0)
   ENDDO
ENDDO

OBST_LOOP_3: DO N=1,N_OBST
   OB => OBSTRUCTION(N)
   DO K=OB%K1,OB%K2
      DO J=OB%J1,OB%J2
         IF (J>OB%J1) CALL DEFINE_EDGE(OB%I1,J,K,-1,2,NM,N)
         IF (J>OB%J1) CALL DEFINE_EDGE(OB%I2,J,K, 1,2,NM,N)
         IF (K>OB%K1) CALL DEFINE_EDGE(OB%I1,J,K,-1,3,NM,N)
         IF (K>OB%K1) CALL DEFINE_EDGE(OB%I2,J,K, 1,3,NM,N)
      ENDDO
   ENDDO
   DO K=OB%K1,OB%K2
      DO I=OB%I1,OB%I2
         IF (I>OB%I1) CALL DEFINE_EDGE(I,OB%J1,K,-2,1,NM,N)
         IF (I>OB%I1) CALL DEFINE_EDGE(I,OB%J2,K, 2,1,NM,N)
         IF (K>OB%K1) CALL DEFINE_EDGE(I,OB%J1,K,-2,3,NM,N)
         IF (K>OB%K1) CALL DEFINE_EDGE(I,OB%J2,K, 2,3,NM,N)
      ENDDO
   ENDDO
   DO J=OB%J1,OB%J2
      DO I=OB%I1,OB%I2
         IF (I>OB%I1) CALL DEFINE_EDGE(I,J,OB%K1,-3,1,NM,N)
         IF (I>OB%I1) CALL DEFINE_EDGE(I,J,OB%K2, 3,1,NM,N)
         IF (J>OB%J1) CALL DEFINE_EDGE(I,J,OB%K1,-3,2,NM,N)
         IF (J>OB%J1) CALL DEFINE_EDGE(I,J,OB%K2, 3,2,NM,N)
      ENDDO
   ENDDO
ENDDO OBST_LOOP_3
 
END SUBROUTINE INITIALIZE_EDGES
 
 
SUBROUTINE INITIALIZE_POISSON_SOLVER
USE POIS, ONLY: H3CZIS,H2CZIS,H3CSIS,H2CYIS
REAL(EB) :: XLM,XMU
INTEGER  :: N,IERR
INTEGER, POINTER :: ITRN,JTRN,KTRN,LBC,MBC,NBC
INTEGER, POINTER, DIMENSION(:) :: NOC
TYPE (VENTS_TYPE), POINTER :: VT
 
!Allocate major arrays
 
ITRN =>M%ITRN 
JTRN =>M%JTRN
KTRN =>M%KTRN
LBC =>M%LBC
MBC =>M%MBC
NBC =>M%NBC
NOC=>TRANS(NM)%NOC
IF (NOC(1)==0 .AND. NOC(2)==0 .AND. NOC(3)==0) M%IPS=0
IF (NOC(1)/=0 .AND. NOC(2)==0 .AND. NOC(3)==0) M%IPS=1
IF (NOC(1)==0 .AND. NOC(2)/=0 .AND. NOC(3)==0) M%IPS=2
IF (NOC(1)==0 .AND. NOC(2)==0 .AND. NOC(3)/=0) M%IPS=3
IF (NOC(1)/=0 .AND. NOC(2)/=0 .AND. NOC(3)==0) M%IPS=4
IF (NOC(1)/=0 .AND. NOC(2)==0 .AND. NOC(3)/=0) M%IPS=5
IF (NOC(1)==0 .AND. NOC(2)/=0 .AND. NOC(3)/=0) M%IPS=6
IF (EVACUATION_ONLY(NM)                      ) M%IPS=7
IF (NOC(1)/=0 .AND. NOC(2)/=0 .AND. NOC(3)/=0) CALL SHUTDOWN('ERROR: Stretch at most 2 coordinate directions')
 
IF (M%IPS<=1 .OR. M%IPS==4) THEN
   ITRN = IBP1
   IF (JBAR>1) JTRN = JBP1
   IF (JBAR==1) JTRN = 1
   KTRN = KBP1
ENDIF
 
IF (M%IPS==2) THEN
   ITRN = JBP1
   JTRN = IBP1
   KTRN = KBP1
   ALLOCATE(M%BZST(JBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZST',IZERO)
   ALLOCATE(M%BZFT(JBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZFT',IZERO)
ENDIF
 
IF (M%IPS==3 .OR. M%IPS==6) THEN
   ITRN = KBP1
   IF (JBAR>1) JTRN = JBP1
   IF (JBAR==1) JTRN = 1
   KTRN = IBP1
   ALLOCATE(M%BXST(KBP1,JTRN),STAT=IZERO)
   CALL ChkMemErr('INIT','BXST',IZERO)
   ALLOCATE(M%BXFT(KBP1,JTRN),STAT=IZERO)
   CALL ChkMemErr('INIT','BXFT',IZERO)
   ALLOCATE(M%BYST(KBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BYST',IZERO)
   ALLOCATE(M%BYFT(KBP1,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BYFT',IZERO)
   ALLOCATE(M%BZST(JTRN,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZST',IZERO)
   ALLOCATE(M%BZFT(JTRN,IBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BZFT',IZERO)
ENDIF
 
IF (M%IPS==5) THEN
   ITRN = IBP1
   JTRN = KBP1
   KTRN = JBP1
   ALLOCATE(M%BXST(KBP1,JBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BXST',IZERO)
   ALLOCATE(M%BXFT(KBP1,JBP1),STAT=IZERO)
   CALL ChkMemErr('INIT','BXFT',IZERO)
ENDIF
 
IF (M%IPS==7) THEN
   ITRN = IBP1
   JTRN = JBP1
   KTRN = 1
ENDIF
 
IF (M%IPS<=3 .OR. M%IPS==7) THEN
   M%LSAVE = (ITRN+1)*JTRN*KTRN+7*ITRN+5*JTRN+6*KTRN+56
   M%LWORK = (ITRN+1)*JTRN*KTRN
ELSE
   N_LOOP: DO N=1,50
      IF ((JTRN+1)<=2**N) EXIT N_LOOP
   ENDDO N_LOOP
   M%LSAVE = KTRN*(6*N*(2**N)+2*N+19)+8*ITRN+7*JTRN+38
   M%LWORK = JTRN*(ITRN*(KTRN+1)+1)
ENDIF
 
ALLOCATE(M%SAVE(-3:M%LSAVE),STAT=IZERO)
CALL ChkMemErr('INIT','SAVE',IZERO)
ALLOCATE(M%SAVE2(-3:M%LSAVE),STAT=IZERO)
CALL ChkMemErr('INIT','SAVE2',IZERO)
ALLOCATE(M%WORK(M%LWORK),STAT=IZERO)
CALL ChkMemErr('INIT','WORK',IZERO)
ALLOCATE(M%PRHS(ITRN,JTRN,KTRN),STAT=IZERO)
CALL ChkMemErr('INIT','PRHS',IZERO)
   IF (KBAR>1) THEN
   IF (JBAR>1) ALLOCATE(M%BXS(JBP1,KBP1),STAT=IZERO)
   IF (JBAR==1) ALLOCATE(M%BXS(1,KBP1)   ,STAT=IZERO)
ELSE
               ALLOCATE(M%BXS(JBP1,1)   ,STAT=IZERO)
ENDIF
CALL ChkMemErr('INIT','BXS',IZERO)
IF (KBAR>1) THEN
   IF (JBAR>1) ALLOCATE(M%BXF(JBP1,KBP1),STAT=IZERO)
   IF (JBAR==1) ALLOCATE(M%BXF(1,KBP1)   ,STAT=IZERO)
ELSE
               ALLOCATE(M%BXF(JBP1,1)   ,STAT=IZERO)
ENDIF
CALL ChkMemErr('INIT','BXF',IZERO)
IF (KBAR>1) THEN
   ALLOCATE(M%BYS(IBP1,KBP1),STAT=IZERO)
ELSE
   ALLOCATE(M%BYS(IBP1,1),STAT=IZERO)
ENDIF
CALL ChkMemErr('INIT','BYS',IZERO)
IF (KBAR>1) THEN
   ALLOCATE(M%BYF(IBP1,KBP1),STAT=IZERO)
ELSE
   ALLOCATE(M%BYF(IBP1,1),STAT=IZERO)
ENDIF
CALL ChkMemErr('INIT','BYF',IZERO)
IF (JBAR>1) ALLOCATE(M%BZS(IBP1,JBP1),STAT=IZERO)
IF (JBAR==1) ALLOCATE(M%BZS(IBP1,1)   ,STAT=IZERO)
CALL ChkMemErr('INIT','BZS',IZERO)
IF (JBAR>1) ALLOCATE(M%BZF(IBP1,JBP1),STAT=IZERO)
IF (JBAR==1) ALLOCATE(M%BZF(IBP1,1)   ,STAT=IZERO)
CALL ChkMemErr('INIT','BZF',IZERO)
 
M%SAVE  = 0._EB
M%WORK  = 0._EB
M%PRHS  = 0._EB
M%BXS   = 0._EB
M%BXF   = 0._EB
M%BYS   = 0._EB
M%BYF   = 0._EB
M%BZS   = 0._EB
M%BZF   = 0._EB
 
! Initialize pressure solver   
 
XLM = 0._EB         ! No Helmholtz equation
XMU = 0._EB         ! No Helmholtz equation
LBC = 3
MBC = 3
NBC = 3
 
VENT_LOOP: DO N=1,M%N_VENT
   VT => M%VENTS(N)
   IF (VT%BOUNDARY_TYPE /= OPEN_BOUNDARY) CYCLE VENT_LOOP
   IF (VT%I1==0 .AND. VT%I2==0) THEN
      IF (LBC==3) LBC = 2
      IF (LBC==4) LBC = 1
   ENDIF
   IF (VT%I1==M%IBAR .AND. VT%I2==M%IBAR) THEN
      IF (LBC==3) LBC = 4
      IF (LBC==2) LBC = 1
   ENDIF
   IF (VT%J1==0 .AND. VT%J2==0) THEN
      IF (MBC==3) MBC = 2
      IF (MBC==4) MBC = 1
   ENDIF
   IF (VT%J1==M%JBAR .AND. VT%J2==M%JBAR) THEN
      IF (MBC==3) MBC = 4
      IF (MBC==2) MBC = 1
   ENDIF
   IF (VT%K1==0 .AND. VT%K2==0) THEN
      IF (NBC==3) NBC = 2
      IF (NBC==4) NBC = 1
   ENDIF
   IF (VT%K1==M%KBAR .AND. VT%K2==M%KBAR) THEN
      IF (NBC==3) NBC = 4
      IF (NBC==2) NBC = 1
   ENDIF
ENDDO VENT_LOOP
 
DO IW=1,M%NEWC
   IF (M%IJKW(9,IW)>0) THEN
      SELECT CASE(M%IJKW(4,IW))
         CASE( 1)
            IF (LBC==3) LBC = 2
            IF (LBC==4) LBC = 1
         CASE(-1)
            IF (LBC==3) LBC = 4
            IF (LBC==2) LBC = 1
         CASE( 2)
            IF (MBC==3) MBC = 2
            IF (MBC==4) MBC = 1
         CASE(-2)
            IF (MBC==3) MBC = 4
            IF (MBC==2) MBC = 1
         CASE( 3)
            IF (NBC==3) NBC = 2
            IF (NBC==4) NBC = 1
         CASE(-3)
            IF (NBC==3) NBC = 4
            IF (NBC==2) NBC = 1
      END SELECT
   ENDIF
ENDDO
 
! User over-rides of Poisson boundary conditions
 
IF (PBC(1,NM)==0 .AND. PBC(2,NM)==0) LBC = 1
IF (PBC(1,NM)==0 .AND. PBC(2,NM)==1) LBC = 2
IF (PBC(1,NM)==1 .AND. PBC(2,NM)==1) LBC = 3
IF (PBC(1,NM)==1 .AND. PBC(2,NM)==0) LBC = 4
IF (PBC(3,NM)==0 .AND. PBC(4,NM)==0) MBC = 1
IF (PBC(3,NM)==0 .AND. PBC(4,NM)==1) MBC = 2
IF (PBC(3,NM)==1 .AND. PBC(4,NM)==1) MBC = 3
IF (PBC(3,NM)==1 .AND. PBC(4,NM)==0) MBC = 4
IF (PBC(5,NM)==0 .AND. PBC(6,NM)==0) NBC = 1
IF (PBC(5,NM)==0 .AND. PBC(6,NM)==1) NBC = 2
IF (PBC(5,NM)==1 .AND. PBC(6,NM)==1) NBC = 3
IF (PBC(5,NM)==1 .AND. PBC(6,NM)==0) NBC = 4
 
! Poisson solver with stretching in the 1st coordinate
 
      SELECT_POISSON_SOLVER: SELECT CASE(M%IPS)
 
         CASE (0:1) SELECT_POISSON_SOLVER
            IF (.NOT.TWO_D) THEN
               CALL H3CZIS(XS,XF,IBAR,LBC,YS,YF,JBAR,MBC,ZS,ZF,KBAR,NBC,M%HX, XLM,ITRN,JTRN,IERR,M%SAVE)
               IF (PRESSURE_CORRECTION) CALL H3CZIS(XS,XF,IBAR,3,YS,YF,JBAR,3,ZS,ZF,KBAR,3,M%HX, XLM,ITRN,JTRN,IERR,M%SAVE2)
            ENDIF
            IF (TWO_D .AND. .NOT.CYLINDRICAL) THEN
               CALL H2CZIS(XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,M%HX,XLM, ITRN,IERR,M%SAVE)
               IF (PRESSURE_CORRECTION) CALL H2CZIS(XS,XF,IBAR,3,ZS,ZF,KBAR,3,M%HX,XLM,ITRN,IERR,M%SAVE2)
            ENDIF
            IF (TWO_D .AND. CYLINDRICAL) THEN
               IF (XS==0._EB .AND. LBC==1) LBC = 5
               IF (XS==0._EB .AND. LBC==2) LBC = 6
               IF (XS==0._EB .AND. LBC==3) LBC = 6
               IF (XS==0._EB .AND. LBC==4) LBC = 5
               CALL H2CYIS(XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,XLM,XMU,ITRN,IERR,M%SAVE)
            ENDIF
         CASE (2) SELECT_POISSON_SOLVER
            CALL H3CZIS(YS,YF,JBAR,MBC,XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,M%HY,XLM,ITRN,JTRN,IERR,M%SAVE)
         CASE (3) SELECT_POISSON_SOLVER
            IF (.NOT.TWO_D) CALL H3CZIS(ZS,ZF,KBAR,NBC,YS,YF,JBAR,MBC,XS,XF,IBAR,LBC,M%HZ,XLM,ITRN,JTRN,IERR,M%SAVE)
            IF (TWO_D) CALL H2CZIS(ZS,ZF,KBAR,NBC,XS,XF,IBAR,LBC,M%HZ,XLM,ITRN,IERR,M%SAVE)
         CASE (4) SELECT_POISSON_SOLVER
            CALL H3CSIS(XS,XF,IBAR,LBC,YS,YF,JBAR,MBC,ZS,ZF,KBAR,NBC,XLM,ITRN,JTRN,IERR,M%SAVE,M%WORK,M%HX,M%HY)
         CASE (5) SELECT_POISSON_SOLVER
            IF (.NOT.TWO_D) CALL H3CSIS(XS,XF,IBAR,LBC,ZS,ZF,KBAR,NBC,YS,YF,JBAR,MBC,XLM,ITRN,JTRN,IERR,M%SAVE,M%WORK,M%HX,M%HZ)
            IF (TWO_D) CALL H2CZIS(ZS,ZF,KBAR,NBC,XS,XF,IBAR,LBC,M%HZ,XLM,ITRN,IERR,M%SAVE)
         CASE (6) SELECT_POISSON_SOLVER
            CALL H3CSIS(ZS,ZF,KBAR,NBC,YS,YF,JBAR,MBC,XS,XF,IBAR,LBC,XLM,ITRN,JTRN,IERR,M%SAVE,M%WORK,M%HZ,M%HY)
         CASE (7) SELECT_POISSON_SOLVER
            CALL H2CZIS(XS,XF,IBAR,LBC,YS,YF,JBAR,MBC,M%HX,XLM,ITRN,IERR,M%SAVE)
 
      END SELECT SELECT_POISSON_SOLVER
 
! Check for errors with Poisson solver initialization
 
IF (IERR/=0) THEN
   WRITE(MESSAGE,'(A,I2,A,I3)') 'ERROR: Poisson initialization error, Number=',IERR, ', Mesh=',NM
   CALL SHUTDOWN(MESSAGE)
ENDIF
 
END SUBROUTINE INITIALIZE_POISSON_SOLVER
 
 
SUBROUTINE INITIAL_NOISE

! Generate random noise at the start of the simulation
 
REAL(EB) :: VFAC,RN
INTEGER  :: I,J,K
 
IF (EVACUATION_ONLY(NM)) RETURN
 
VFAC = 0.005_EB
 
DO K=0,M%KBAR
   DO J=0,M%JBAR
      LOOP_1: DO I=1,M%IBAR
         IF (M%SOLID(M%CELL_INDEX(I,J,K))   .OR. M%SOLID(M%CELL_INDEX(I,J,K+1)) .OR. &
             M%SOLID(M%CELL_INDEX(I,J+1,K)) .OR. M%SOLID(M%CELL_INDEX(I,J+1,K+1)))  CYCLE LOOP_1
         CALL RANDOM_NUMBER(RN)
         RN = VFAC*(-1._EB + 2._EB*RN)*M%DXMIN
         M%W(I,J,K)   = M%W(I,J,K)   - RN*M%RDY(J)
         M%W(I,J+1,K) = M%W(I,J+1,K) + RN*M%RDY(J+1)
         M%V(I,J,K)   = M%V(I,J,K)   + RN*M%RDZ(K)
         M%V(I,J,K+1) = M%V(I,J,K+1) - RN*M%RDZ(K+1)
      ENDDO LOOP_1
   ENDDO
ENDDO
DO K=0,M%KBAR
   DO J=1,M%JBAR
      LOOP_2: DO I=0,M%IBAR
         IF (M%SOLID(M%CELL_INDEX(I,J,K))   .OR. M%SOLID(M%CELL_INDEX(I,J,K+1)) .OR. &
             M%SOLID(M%CELL_INDEX(I+1,J,K)) .OR. M%SOLID(M%CELL_INDEX(I+1,J,K+1)))  CYCLE LOOP_2
         CALL RANDOM_NUMBER(RN)
         RN = VFAC*(-1._EB + 2._EB*RN)*M%DXMIN
         M%W(I,J,K)   = M%W(I,J,K)   - RN*M%RDX(I)*M%R(I)*M%RRN(I)
         M%W(I+1,J,K) = M%W(I+1,J,K) + RN*M%RDX(I+1)*M%R(I)*M%RRN(I+1)
         M%U(I,J,K)   = M%U(I,J,K)   + RN*M%RDZ(K)
         M%U(I,J,K+1) = M%U(I,J,K+1) - RN*M%RDZ(K+1)
      ENDDO LOOP_2
   ENDDO
ENDDO
DO K=1,M%KBAR
   DO J=0,M%JBAR
      LOOP_3: DO I=0,M%IBAR
         IF (M%SOLID(M%CELL_INDEX(I,J,K))   .OR. M%SOLID(M%CELL_INDEX(I,J+1,K)) .OR. &
             M%SOLID(M%CELL_INDEX(I+1,J,K)) .OR. M%SOLID(M%CELL_INDEX(I+1,J+1,K)))  CYCLE LOOP_3
         CALL RANDOM_NUMBER(RN)
         RN = VFAC*(-1._EB + 2._EB*RN)*M%DXMIN
         M%V(I,J,K)   = M%V(I,J,K)   - RN*M%RDX(I)
         M%V(I+1,J,K) = M%V(I+1,J,K) + RN*M%RDX(I+1)
         M%U(I,J,K)   = M%U(I,J,K)   + RN*M%RDY(J)
         M%U(I,J+1,K) = M%U(I,J+1,K) - RN*M%RDY(J+1)
      ENDDO LOOP_3
   ENDDO
ENDDO
 
END SUBROUTINE INITIAL_NOISE
 
 
SUBROUTINE INITIALIZE_DEVC

INTEGER :: III
TYPE (DEVICE_TYPE), POINTER :: DV
 
DEVICE_LOOP: DO N=1,N_DEVC
   DV => DEVICE(N)
   IF (NM/=DV%MESH .OR. DV%OUTPUT_INDEX>0)  CYCLE DEVICE_LOOP
   II  = GINV(DV%X-M%XS,1,NM)*M%RDXI   + 1._EB
   JJ  = GINV(DV%Y-M%YS,2,NM)*M%RDETA  + 1._EB
   KK  = GINV(DV%Z-M%ZS,3,NM)*M%RDZETA + 1._EB
   IOR = DV%IOR
   CALL GET_IW(II,JJ,KK,IOR,IW)
   IF (IW>0) THEN
      DV%IW = IW
      IF (DV%OUTPUT_INDEX==-6) THEN
         IBC = M%IJKW(5,IW)
         IF (SURFACE(IBC)%THERMAL_BC_INDEX /= THERMALLY_THICK) THEN
            WRITE(MESSAGE,'(A,I3,A)') 'DEViCe ',N, ' must be associated with a heat-conducting surface'
            CALL SHUTDOWN(MESSAGE)
         ENDIF
         DV%I_DEPTH = SURFACE(IBC)%N_CELLS
         DO III=SURFACE(IBC)%N_CELLS,1,-1
            IF (DV%DEPTH<=SURFACE(IBC)%X_S(III)) DV%I_DEPTH = III
         ENDDO
      ENDIF
   ELSE
      WRITE(MESSAGE,'(A,I4,A)') 'Reposition DEVC No.',DV%ORDINAL, '. FDS cannot determine which boundary cell to assign'
      CALL SHUTDOWN(MESSAGE)
   ENDIF
ENDDO DEVICE_LOOP
 
END SUBROUTINE INITIALIZE_DEVC
 
 
SUBROUTINE INITIALIZE_PROF

INTEGER :: NN
LOGICAL :: SUCCESS
TYPE (PROFILE_TYPE), POINTER :: PF

PROF_LOOP: DO N=1,N_PROF
   PF => PROFILE(N)
   IF (NM/=PF%MESH) CYCLE PROF_LOOP
   II  = GINV(PF%X-M%XS,1,NM)*M%RDXI   + 1._EB
   JJ  = GINV(PF%Y-M%YS,2,NM)*M%RDETA  + 1._EB
   KK  = GINV(PF%Z-M%ZS,3,NM)*M%RDZETA + 1._EB
   IOR = PF%IOR
   CALL GET_IW(II,JJ,KK,IOR,IW)
   IF (IW>0) THEN
      PF%IW = IW
      SF => SURFACE(IJKW(5,IW))