Merge pull request #14206 from gforney/smoke

FDS source: output 3D smoke using density

Author: gforney

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -10070,7 +10070,6 @@ \section{Plot3D Static Data Dumps}
 
 \section{SMOKE3D: Realistic Smoke and Fire}
 \label{info:SMOKE3D}
-\label{info:smoke16}
 
 For any simulation involving combustion, FDS automatically creates three output files that are rendered by Smokeview as realistic looking smoke, fire, and/or hot gas. By default, the output quantities are the \ct{'DENSITY'} of \ct{'SOOT'}, the \ct{'HRRPUV'} (Heat Release Rate Per Unit Volume), and the \ct{TEMPERATURE} of the gases. These quantities are rendered as semi-transparent slices that give you the illusion of looking through smoke or fire. The greater the value of the quantity, the darker the slices will appear en masse. For \ct{HRRPUV}, the lowest value rendered is:
 \be
@@ -10085,7 +10084,7 @@ \section{SMOKE3D: Realistic Smoke and Fire}
 \end{lstlisting}
 The \ct{MASS_EXTINCTION_COEFFICIENT} is passed to Smokeview to be used for visualization.
 
-FDS outputs 3D smoke quantities as 8 bit integers compressed using run length encoding.  Soot density, HRRPUV or temperatures are first scaled to 8 bit integers (soot density is converted to an opacity first). Repeated integers are replaced by $n$I where $n$ is the number of repeats and I is the value repeated.  When higher resolution is desired, set \ct{SMOKE3D_16=T} on the \ct{MISC} line to output these quantities as 16 bit integers.  For each time frame, FDS then outputs the time, the minimum value, the maximum value, and 16 bit integers scaled between the minimum and maximum value.  This data is then output to a file with extension \ct{.s16} rather than \ct{.s3d} . These 16 bit integer values are not compressed.
+FDS outputs 3D smoke quantities as 8 bit integers compressed using run length encoding.  Soot density, HRRPUV or temperatures are first scaled to 8 bit integers (soot density is converted to an opacity first). Repeated integers are replaced by $n$I where $n$ is the number of repeats and I is the value repeated.  
 
 \newpage
 
@@ -12655,7 +12654,6 @@ \section{\texorpdfstring{{\tt MISC}}{MISC} (Miscellaneous Parameters)}
 \ct{RND_SEED}                                     & Integer      & Section~\ref{info:restart}              &            &                   \\ \hline
 \ct{SC}                                            & Real         & Section~\ref{info:NOISE}                &            & 0.5               \\ \hline
 \ct{SIMULATION_MODE}                              & Character    & Section~\ref{Sim_Mode}                  &            & \ct{'VLES'}      \\ \hline
-\ct{SMOKE3D_16}                                   & Logical      & Section~\ref{info:smoke16}              &            & \ct{F}           \\ \hline
 \ct{SMOKE_ALBEDO}                                 & Real         & Reference~\cite{Smokeview_Users_Guide}  &            & 0.3               \\ \hline
 \ct{SOLID_PHASE_ONLY}                            & Logical      & Section~\ref{solid_phase_verification}  &            & \ct{F}           \\ \hline
 \ct{SOOT_DENSITY}                                 & Real         & Section~\ref{info:deposition}           &            & 1800              \\ \hline

Source/cons.f90

@@ -267,7 +267,6 @@ MODULE GLOBAL_CONSTANTS
 LOGICAL :: NO_PRESSURE_ZONES=.FALSE.                !< Flag to suppress pressure zones
 LOGICAL :: CTRL_DIRECT_FORCE=.FALSE.                !< Allow adjustable direct force via CTRL logic
 LOGICAL :: REACTING_THIN_OBSTRUCTIONS=.FALSE.       !< Thin obstructions that off-gas are present
-LOGICAL :: SMOKE3D_16=.FALSE.                       !< Output 3D smoke values using 16 bit integers
 LOGICAL :: CHECK_BOUNDARY_ONE_D_ARRAYS=.FALSE.      !< Flag that indicates that ONE_D array dimensions need to be checked
 LOGICAL :: TENSOR_DIFFUSIVITY=.FALSE.               !< If true, use experimental tensor diffusivity model for spec and tmp
 LOGICAL :: OXPYRO_MODEL=.FALSE.                     !< Flag to use oxidative pyrolysis mass transfer model
@@ -736,17 +735,17 @@ MODULE GLOBAL_CONSTANTS
 REAL(EB), ALLOCATABLE, DIMENSION(:) :: VIEW_ANGLE_FACTOR        !< View angle area ORIENTATION_VECTOR
 INTEGER :: N_ORIENTATION_VECTOR                               !< Number of ORIENTATION_VECTORs
 
-INTEGER :: TGA_MESH_INDEX=HUGE(INTEGER_ONE)  !< Mesh for the special TGA calculation
-INTEGER :: TGA_SURF_INDEX=-100               !< Surface properties to use for special TGA calculation
-INTEGER :: TGA_WALL_INDEX=-100               !< Wall index to use for special TGA calculation
-INTEGER :: TGA_PARTICLE_INDEX=-100           !< Particle index to use for special TGA calculation
-REAL(EB) :: TGA_DT=0.01_EB                   !< Time step (s) to use for special TGA calculation
-REAL(EB) :: TGA_DUMP=1._EB                   !< Temperature output interval (K), starting at TMPA, to use for special TGA calculation
-REAL(EB) :: TGA_HEATING_RATE=5._EB           !< Heat rate (K/min) to use for special TGA calculation
-REAL(EB) :: TGA_FINAL_TEMPERATURE=800._EB    !< Final Temperature (C) to use for special TGA calculation
-REAL(EB) :: TGA_CONVERSION_FACTOR=1._EB      !< Conversion factor for TGA output
-REAL(EB) :: MCC_CONVERSION_FACTOR=1._EB      !< Conversion factor for MCC output
-REAL(EB) :: DSC_CONVERSION_FACTOR=1._EB      !< Conversion factor for DSC output
+INTEGER :: TGA_MESH_INDEX=HUGE(INTEGER_ONE) !< Mesh for the special TGA calculation
+INTEGER :: TGA_SURF_INDEX=-100              !< Surface properties to use for special TGA calculation
+INTEGER :: TGA_WALL_INDEX=-100              !< Wall index to use for special TGA calculation
+INTEGER :: TGA_PARTICLE_INDEX=-100          !< Particle index to use for special TGA calculation
+REAL(EB) :: TGA_DT=0.01_EB                  !< Time step (s) to use for special TGA calculation
+REAL(EB) :: TGA_DUMP=1._EB                  !< Temperature output interval (K), starting at TMPA, to use for special TGA calculation
+REAL(EB) :: TGA_HEATING_RATE=5._EB          !< Heat rate (K/min) to use for special TGA calculation
+REAL(EB) :: TGA_FINAL_TEMPERATURE=800._EB   !< Final Temperature (C) to use for special TGA calculation
+REAL(EB) :: TGA_CONVERSION_FACTOR=1._EB     !< Conversion factor for TGA output
+REAL(EB) :: MCC_CONVERSION_FACTOR=1._EB     !< Conversion factor for MCC output
+REAL(EB) :: DSC_CONVERSION_FACTOR=1._EB     !< Conversion factor for DSC output
 
 LOGICAL :: IBLANK_SMV=.TRUE.  !< Parameter passed to smokeview (in .smv file) to control generation of blockages
 

Source/dump.f90

@@ -393,13 +393,8 @@ SUBROUTINE ASSIGN_FILE_NAMES
    ALLOCATE(FN_TERRAIN(NMESHES))
    ALLOCATE(LU_TERRAIN(NMESHES))
 ENDIF
-IF (.NOT.SMOKE3D_16) THEN
-   ALLOCATE(FN_SMOKE3D(N_SMOKE3D*2,NMESHES))
-   ALLOCATE(LU_SMOKE3D(N_SMOKE3D*2,NMESHES))
-ELSE
-   ALLOCATE(FN_SMOKE3D(N_SMOKE3D*3,NMESHES)) ! also allocate unit numbers and file names for the size files
-   ALLOCATE(LU_SMOKE3D(N_SMOKE3D*3,NMESHES)) ! and the SMOKE3D_16 files (16 bit integers)
-ENDIF
+ALLOCATE(FN_SMOKE3D(N_SMOKE3D*3,NMESHES)) ! also allocate unit numbers and file names for the size files
+ALLOCATE(LU_SMOKE3D(N_SMOKE3D*3,NMESHES)) ! and the SMOKE3D_DENSITY files
 ALLOCATE(FN_PART(2*NMESHES))
 ALLOCATE(LU_PART(2*NMESHES))
 ALLOCATE(FN_CORE(NMESHES))
@@ -445,10 +440,8 @@ SUBROUTINE ASSIGN_FILE_NAMES
       WRITE(FN_SMOKE3D(N,NM),  '(A,A,I0,A,I0,A)') TRIM(RESULTS_DIR)//TRIM(CHID),'_',NM,'_',N,'.s3d'
       LU_SMOKE3D(N+N_SMOKE3D,NM) = GET_FILE_NUMBER()
       WRITE(FN_SMOKE3D(N+N_SMOKE3D,NM),'(A,A,I0,A,I0,A)') TRIM(RESULTS_DIR)//TRIM(CHID),'_',NM,'_',N,'.s3d.sz'
-      IF (SMOKE3D_16) THEN
-         LU_SMOKE3D(N+2*N_SMOKE3D,NM) = GET_FILE_NUMBER()
-         WRITE(FN_SMOKE3D(N+2*N_SMOKE3D,NM),  '(A,A,I0,A,I0,A)') TRIM(RESULTS_DIR)//TRIM(CHID),'_',NM,'_',N,'.s16'
-      ENDIF
+      LU_SMOKE3D(N+2*N_SMOKE3D,NM) = GET_FILE_NUMBER()
+      WRITE(FN_SMOKE3D(N+2*N_SMOKE3D,NM),  '(A,A,I0,A,I0,A)') TRIM(RESULTS_DIR)//TRIM(CHID),'_',NM,'_',N,'.s3dd'
    ENDDO
 
    ! Slice Files
@@ -906,7 +899,6 @@ SUBROUTINE INITIALIZE_MESH_DUMPS(NM)
 TYPE(PATCH_TYPE), POINTER :: PA
 TYPE(MESH_TYPE), POINTER :: M4
 INTEGER :: CC_VAL,NSTEPS
-INTEGER :: SMOKE3D_16_COMPRESS, SMOKE3D_16_VERSION
 LOGICAL :: OVERLAPPING_X,OVERLAPPING_Y,OVERLAPPING_Z
 TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC
 TYPE (RAD_FILE_TYPE), POINTER :: RF
@@ -1075,11 +1067,9 @@ SUBROUTINE INITIALIZE_MESH_DUMPS(NM)
          OPEN(LU_SMOKE3D(N+N_SMOKE3D,NM),FILE=FN_SMOKE3D(N+N_SMOKE3D,NM),FORM='FORMATTED',STATUS='REPLACE')
          WRITE(LU_SMOKE3D(N+N_SMOKE3D,NM),*) INTEGER_ZERO
          CLOSE(LU_SMOKE3D(N+N_SMOKE3D,NM))
-         IF (SMOKE3D_16) THEN
-            SMOKE3D_16_COMPRESS = 0
-            SMOKE3D_16_VERSION = 1
+         IF (SMOKE3D_FILE(N)%DISPLAY_TYPE=='GAS') THEN
             OPEN(LU_SMOKE3D(N+2*N_SMOKE3D,NM), FILE=FN_SMOKE3D(N+2*N_SMOKE3D,NM),FORM='UNFORMATTED',STATUS='REPLACE')
-            WRITE(LU_SMOKE3D(N+2*N_SMOKE3D,NM)) INTEGER_ONE,SMOKE3D_16_VERSION,SMOKE3D_16_COMPRESS,M%IBP1,M%JBP1,M%KBP1
+            WRITE(LU_SMOKE3D(N+2*N_SMOKE3D,NM)) INTEGER_ONE,INTEGER_ZERO,0,M%IBAR,0,M%JBAR,0,M%KBAR
             CLOSE(LU_SMOKE3D(N+2*N_SMOKE3D,NM))
          ENDIF
          IF (M%N_STRINGS+5>M%N_STRINGS_MAX) CALL RE_ALLOCATE_STRINGS(NM)
@@ -4451,7 +4441,7 @@ SUBROUTINE DUMP_SMOKE3D(T,DT,NM)
 
    ALLOCATE(QQ_PACK(IBP1*JBP1*KBP1))
    QQ_PACK = PACK(QQ(0:IBAR,0:JBAR,0:KBAR,1),MASK=.TRUE.)
-   CALL SMOKE3D_TO_FILE(NM,STIME,DXX,N,QQ_PACK,IBP1,JBP1,KBP1,SMOKE3D_16)
+   CALL SMOKE3D_TO_FILE(NM,STIME,DXX,N,QQ_PACK,IBP1,JBP1,KBP1)
    DEALLOCATE(QQ_PACK)
 
 ENDDO DATA_FILE_LOOP

Source/read.f90

@@ -1743,7 +1743,7 @@ SUBROUTINE READ_MISC
                 P_INF,RAMP_GX,RAMP_GY,RAMP_GZ,&
                 RAMP_UX,RAMP_UY,RAMP_UZ,RAMP_VX,RAMP_VY,RAMP_VZ,RAMP_WX,RAMP_WY,RAMP_WZ,&
                 RESTART,RESTART_CHID,SC,&
-                RND_SEED,SIMULATION_MODE,SMOKE3D_16,SMOKE_ALBEDO,SOLID_PHASE_ONLY,SOOT_DENSITY,SOOT_OXIDATION,&
+                RND_SEED,SIMULATION_MODE,SMOKE_ALBEDO,SOLID_PHASE_ONLY,SOOT_DENSITY,SOOT_OXIDATION,&
                 TAU_DEFAULT,TENSOR_DIFFUSIVITY,TERRAIN_IMAGE,TEST_NEW_CHAR_MODEL,TEXTURE_ORIGIN,&
                 THERMOPHORETIC_DEPOSITION,THERMOPHORETIC_SETTLING,THICKEN_OBSTRUCTIONS,&
                 TMPA,TURBULENCE_MODEL,TURBULENT_DEPOSITION,UVW_FILE,&
@@ -1800,7 +1800,6 @@ SUBROUTINE READ_MISC
 MAX_LEAK_PATHS       = 200
 FLUX_LIMITER         = 'null'
 LES_FILTER_TYPE      = 'null'
-SMOKE3D_16           = .FALSE.
 
 ! Initial read of the MISC line
 

Source/smvv.f90

@@ -921,22 +921,17 @@ END SUBROUTINE SLICE_TO_RLEFILE
 
 !> \brief Routine to generate SMOKE3D compressed data to send to Smokeview
 
-SUBROUTINE SMOKE3D_TO_FILE(NM,TIME,DX,SMOKE3D_INDEX,VALS,NX,NY,NZ,SMOKE3D_16_FLAG)
+SUBROUTINE SMOKE3D_TO_FILE(NM,TIME,DX,SMOKE3D_INDEX,VALS,NX,NY,NZ)
 
 USE OUTPUT_DATA, ONLY: SMOKE3D_TYPE,SMOKE3D_FILE,N_SMOKE3D
 USE GLOBAL_CONSTANTS, ONLY: TMPA,TMPM,LU_SMOKE3D,FN_SMOKE3D,TEMP_MAX_SMV,HRRPUV_MAX_SMV
 INTEGER, INTENT(IN) :: NX,NY,NZ,NM,SMOKE3D_INDEX
-LOGICAL, INTENT(IN) :: SMOKE3D_16_FLAG
 REAL(FB), INTENT(IN) :: TIME, DX
 REAL(FB), INTENT(IN), DIMENSION(NX*NY*NZ) :: VALS
-CHARACTER(LEN=1), DIMENSION(:), ALLOCATABLE :: BUFFER_IN, BUFFER_OUT
-CHARACTER(LEN=1), DIMENSION(:), ALLOCATABLE, TARGET :: BUFFER16_IN
-REAL(FB) :: FACTOR,TEMP_MIN,VAL_FDS,VAL_SMV,MAX_VAL
-INTEGER :: I,NCHARS_OUT,NVALS,NCHARS_IN
-INTEGER :: VAL16, FACTOR16, NVALS16_OUT
-INTEGER :: VAL_LOW, VAL_HIGH
+CHARACTER(LEN=1), DIMENSION(:), ALLOCATABLE :: BUFFER_IN, BUFFER_OUT, BUFFER_DENSITY_IN, BUFFER_DENSITY_OUT
+REAL(FB) :: FACTOR,TEMP_MIN,VAL_FDS,VAL_SMV,MAX_VAL,MAX_DENSITY_VAL
+INTEGER :: I,NCHARS_OUT,NVALS,NCHARS_IN,NCHARS_DENSITY_OUT
 TYPE(SMOKE3D_TYPE), POINTER :: S3
-REAL(FB) :: SMOKE3D_16_VALMIN, SMOKE3D_16_VALMAX
 
 S3 => SMOKE3D_FILE(SMOKE3D_INDEX)
 
@@ -950,13 +945,30 @@ SUBROUTINE SMOKE3D_TO_FILE(NM,TIME,DX,SMOKE3D_INDEX,VALS,NX,NY,NZ,SMOKE3D_16_FLA
 
 IF (S3%DISPLAY_TYPE=='GAS') THEN
 
+   MAX_DENSITY_VAL=0.0_FB
    FACTOR=-REAL(S3%MASS_EXTINCTION_COEFFICIENT,FB)*DX
    DO I=1,NVALS
       VAL_FDS = MAX(0.0_FB,VALS(I))
       VAL_SMV = 254*(1.0_FB-EXP(FACTOR*VAL_FDS))
       BUFFER_IN(I) = CHAR(NINT(VAL_SMV))
       MAX_VAL = MAX(VAL_SMV,MAX_VAL)  ! If MAX_VAL=0, soot in mesh is completely transparent
+      MAX_DENSITY_VAL = MAX(VALS(I),MAX_DENSITY_VAL)
    ENDDO
+   ALLOCATE(BUFFER_DENSITY_IN(NVALS))
+   ALLOCATE(BUFFER_DENSITY_OUT(NVALS))
+   IF (MAX_DENSITY_VAL .LE. 0.0_FB) THEN
+      MAX_DENSITY_VAL = 0.0_FB
+      DO I=1,NVALS
+         BUFFER_DENSITY_IN(I) = CHAR(0)
+      ENDDO
+   ELSE
+      DO I=1,NVALS
+         VAL_FDS = 254.0_FB*(VALS(I)/MAX_DENSITY_VAL)
+         VAL_FDS = MAX(  0.0_FB, VAL_FDS)
+         VAL_FDS = MIN(254.0_FB, VAL_FDS)
+         BUFFER_DENSITY_IN(I) = CHAR(NINT(VAL_FDS))
+      ENDDO
+   ENDIF
 
 ELSEIF (S3%DISPLAY_TYPE=='FIRE') THEN
 
@@ -988,7 +1000,7 @@ SUBROUTINE SMOKE3D_TO_FILE(NM,TIME,DX,SMOKE3D_INDEX,VALS,NX,NY,NZ,SMOKE3D_16_FLA
 
 ! Write size information to a text file
 
-IF (.NOT.SMOKE3D_16_FLAG) THEN ! write out an extra column below if SMOKE3D_16_FLAG is true
+IF (S3%DISPLAY_TYPE.NE.'GAS') THEN ! write out an extra column below if SMOKE3D_DENSITY_FLAG is true
    OPEN(LU_SMOKE3D(SMOKE3D_INDEX+N_SMOKE3D,NM),FILE=FN_SMOKE3D(SMOKE3D_INDEX+N_SMOKE3D,NM),&
         FORM='FORMATTED',STATUS='OLD',POSITION='APPEND')
    WRITE(LU_SMOKE3D(SMOKE3D_INDEX+N_SMOKE3D,NM),*) TIME,NCHARS_IN,NCHARS_OUT,MAX_VAL
@@ -1002,43 +1014,24 @@ SUBROUTINE SMOKE3D_TO_FILE(NM,TIME,DX,SMOKE3D_INDEX,VALS,NX,NY,NZ,SMOKE3D_16_FLA
 WRITE(LU_SMOKE3D(SMOKE3D_INDEX,NM)) NCHARS_IN,NCHARS_OUT
 WRITE(LU_SMOKE3D(SMOKE3D_INDEX,NM)) (BUFFER_OUT(I),I=1,NCHARS_OUT)
 CLOSE(LU_SMOKE3D(SMOKE3D_INDEX,NM))
-! write out data as 2 byte integers
-
-IF (SMOKE3D_16_FLAG) THEN
-   NVALS16_OUT = NVALS   ! after compression is added, NVALS16_OUT will be smaller than NVALS
-   ALLOCATE(BUFFER16_IN(2*NVALS))
-   SMOKE3D_16_VALMIN = VALS(1)
-   SMOKE3D_16_VALMAX = SMOKE3D_16_VALMIN
-   DO I=2,NVALS
-      SMOKE3D_16_VALMIN = MIN(SMOKE3D_16_VALMIN,VALS(I))
-      SMOKE3D_16_VALMAX = MAX(SMOKE3D_16_VALMAX,VALS(I))
-   ENDDO
-   IF (SMOKE3D_16_VALMIN == SMOKE3D_16_VALMAX) THEN
-      SMOKE3D_16_VALMAX = SMOKE3D_16_VALMIN + 1.0_FB
-   ENDIF
-   FACTOR16 = 2**16 - 1
-   DO I=1,NVALS
-      VAL16 = INT(FACTOR16*(VALS(I)-SMOKE3D_16_VALMIN)/(SMOKE3D_16_VALMAX-SMOKE3D_16_VALMIN))
-      VAL16 = MIN(MAX(0, VAL16), FACTOR16)
-      VAL_LOW = MOD(VAL16, 256)
-      VAL_HIGH = VAL16/256
-      BUFFER16_IN(2*I-1) = CHAR(VAL_HIGH)  
-      BUFFER16_IN(2*i)   = CHAR(VAL_LOW)
-   ENDDO
 
+IF ( S3%DISPLAY_TYPE=='GAS') THEN
+! Pack the values into a compressed buffer
+
+   CALL RLE(BUFFER_DENSITY_IN,NCHARS_IN,BUFFER_DENSITY_OUT,NCHARS_DENSITY_OUT)
    OPEN(LU_SMOKE3D(SMOKE3D_INDEX+N_SMOKE3D,NM),FILE=FN_SMOKE3D(SMOKE3D_INDEX+N_SMOKE3D,NM),&
         FORM='FORMATTED',STATUS='OLD',POSITION='APPEND')
-   WRITE(LU_SMOKE3D(SMOKE3D_INDEX+N_SMOKE3D,NM),*) TIME,NCHARS_IN,NCHARS_OUT,MAX_VAL,NVALS16_OUT
+   WRITE(LU_SMOKE3D(SMOKE3D_INDEX+N_SMOKE3D,NM),*) TIME,NCHARS_IN,NCHARS_OUT,MAX_VAL,NCHARS_DENSITY_OUT,MAX_DENSITY_VAL
    CLOSE(LU_SMOKE3D(SMOKE3D_INDEX+N_SMOKE3D,NM))
 
    OPEN(LU_SMOKE3D(SMOKE3D_INDEX+2*N_SMOKE3D,NM),FILE=FN_SMOKE3D(SMOKE3D_INDEX+2*N_SMOKE3D,NM),FORM='UNFORMATTED',&
         STATUS='OLD',POSITION='APPEND')
-   WRITE(LU_SMOKE3D(SMOKE3D_INDEX + 2*N_SMOKE3D,NM)) TIME
-   WRITE(LU_SMOKE3D(SMOKE3D_INDEX + 2*N_SMOKE3D,NM)) NVALS,NVALS16_OUT,SMOKE3D_16_VALMIN,SMOKE3D_16_VALMAX
-   ! output (2**16 - 1)(val-valmin)/(valmax-valmin)
-   WRITE(LU_SMOKE3D(SMOKE3D_INDEX + 2*N_SMOKE3D,NM)) (BUFFER16_IN(I),I=1,2*NVALS)
-   CLOSE(LU_SMOKE3D(SMOKE3D_INDEX + 2*N_SMOKE3D,NM))
-   DEALLOCATE(BUFFER16_IN)
+   WRITE(LU_SMOKE3D(SMOKE3D_INDEX+2*N_SMOKE3D,NM)) TIME
+   WRITE(LU_SMOKE3D(SMOKE3D_INDEX+2*N_SMOKE3D,NM)) NCHARS_IN,NCHARS_DENSITY_OUT
+   WRITE(LU_SMOKE3D(SMOKE3D_INDEX+2*N_SMOKE3D,NM)) (BUFFER_DENSITY_OUT(I),I=1,NCHARS_DENSITY_OUT)
+   CLOSE(LU_SMOKE3D(SMOKE3D_INDEX+2*N_SMOKE3D,NM))
+   DEALLOCATE(BUFFER_DENSITY_IN)
+   DEALLOCATE(BUFFER_DENSITY_OUT)
 ENDIF
 
 DEALLOCATE(BUFFER_IN)