File: RELATIVE:/../../../mfix.git/model/des/set_bc_dem_mi.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: SET_BC_DEM_MI                                           !
4     !  Author: J.Musser                                   Date: 23-Nov-09  !
5     !                                                                      !
6     !  Purpose:                                                            !
7     !                                                                      !
8     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
9           SUBROUTINE SET_BC_DEM_MI
10     
11     !-----------------------------------------------
12     ! Modules
13     !-----------------------------------------------
14           USE bc
15           USE compar
16           USE constant
17           USE des_allocate
18           USE des_bc
19           USE discretelement
20           USE error_manager
21           USE funits
22           USE geometry
23           USE indices
24           USE param
25           USE param1
26           USE physprop
27           USE run
28           USE toleranc
29     
30           IMPLICIT NONE
31     
32           INTEGER :: BCV
33     
34     !-----------------------------------------------
35     ! Local variables
36     !-----------------------------------------------
37           INTEGER BCV_I      ! BC loop counter
38           INTEGER M, MM           ! Mass phase loop counter
39           INTEGER RANGE_TOP, RANGE_BOT ! Dummy values
40           INTEGER PHASE_CNT        ! Number of solid phases at bc
41           INTEGER PHASE_LIST(DIM_M) ! List of phases used in current bc
42     
43     ! the number of particles injected in a solids time step
44           DOUBLE PRECISION NPMpSEC(DES_MMAX) ! For solid phase m
45           DOUBLE PRECISION NPpSEC
46           DOUBLE PRECISION NPpDT        ! Total for BC
47           DOUBLE PRECISION SCALED_VAL
48           DOUBLE PRECISION MAX_DIA ! Max diameter of incoming particles at bc
49     
50           DOUBLE PRECISION :: EPs_ERR
51           DOUBLE PRECISION :: VOL_FLOW
52     
53           LOGICAL, parameter :: setDBG = .FALSE.
54           LOGICAL :: dFlag
55     
56           LOGICAL :: FATAL
57     
58     ! Temp inlet velocity for solids phase M
59           DOUBLE PRECISION VEL_TMP(DIM_M)
60           DOUBLE PRECISION EPs_TMP(DIM_M)
61     
62     ! Minimum/maximum solids velocity at inlet.  Also used in the iterative
63     ! steps as the starting and ending velocities
64           DOUBLE PRECISION  MAX_VEL
65     
66           DOUBLE PRECISION  MINIPV, MAXIPV
67     
68           INTEGER :: OCCUPANTS
69     ! jump_here
70     
71     !-----------------------------------------------
72     
73           CALL INIT_ERR_MSG("SET_BC_DEM_MI")
74     
75           dFlag = (DMP_LOG .AND. setDBG)
76           if(dFlag) write(*,"(2/,2x,'DEM inlet count: ',I4)") DEM_BCMI
77     
78     ! Allocate the MI data structures for NEW and RES2 cases. Allocation for
79     ! RES1 cases is done prior to reading the RES file.
80           IF(RUN_TYPE /= 'RESTART_1') CALL ALLOCATE_DEM_MI
81     
82     ! Loop over BCs that flagged for DEM mass inflow.
83           DO BCV_I = 1, DEM_BCMI
84     
85     ! Get the user defined BC ID.
86              BCV = DEM_BCMI_MAP(BCV_I)
87     
88              if(dFlag) write(*,"(2/,'Setting DEM_MI:',I3)") BCV_I
89     
90     ! The number of mass phases at this inlet.  While a system may be
91     ! polydisperse, the inlet could consist of a single mass phase
92              PHASE_CNT = 0
93     ! The mass phase indices of incoming particles at this inlet
94              PHASE_LIST(:) = -1
95     ! The max diameter of incoming particles at this inlet
96              MAX_DIA = ZERO
97     
98     ! Determine if the inlet is mono or polydisperse
99              DO M=1, SMAX + DES_MMAX
100                 IF(SOLIDS_MODEL(M) /= 'DEM') CYCLE
101                 IF(BC_ROP_s(BCV,M) == UNDEFINED) CYCLE
102                 IF(COMPARE(BC_ROP_s(BCV,M),ZERO)) CYCLE
103                 PHASE_CNT = PHASE_CNT + 1
104                 PHASE_LIST(PHASE_CNT) = M-SMAX
105                 MAX_DIA = MAX(MAX_DIA,DES_D_P0(M-SMAX))
106              ENDDO
107     ! Set the polydispersity flag.
108              DEM_MI(BCV_I)%POLYDISPERSE = (PHASE_CNT > 1)
109     
110     ! Layout the feed pattern.
111              CALL LAYOUT_MI_DEM(BCV, BCV_I, MAX_DIA)
112     
113     ! Initialize
114              MAX_VEL = ZERO
115              NPMpSEC(:) = ZERO
116     
117     ! Calculate the individual velocities for each solid phase
118              DO MM = 1, PHASE_CNT
119                 M = PHASE_LIST(MM)
120     
121     ! Pull off the BC velocity normal to the flow plane.
122                 SELECT CASE(BC_PLANE(BCV))
123                 CASE('N','S'); VEL_TMP(M) = abs(BC_V_s(BCV,M+SMAX))
124                 CASE('E','W'); VEL_TMP(M) = abs(BC_U_s(BCV,M+SMAX))
125                 CASE('T','B'); VEL_TMP(M) = abs(BC_W_s(BCV,M+SMAX))
126                 END SELECT
127     
128     ! Check for min/max inlet velocity
129                 MAX_VEL = MAX(ABS(VEL_TMP(M)), MAX_VEL)
130     ! Calculate volumetric flow rate to convert to particle count. BC_AREA
131     ! was already corrected for cut cells and velocity was recalculated
132     ! to ensure user-specified mass or volumetric flow rates.
133                 VOL_FLOW = VEL_TMP(M) * BC_AREA(BCV) * BC_EP_S(BCV,M+SMAX)
134     ! Calculate the number of particles of mass phase M are injected per
135     ! second for each solid phase present at the boundary
136                 NPMpSEC(M) = VOL_FLOW / (PI/6.d0*DES_D_P0(M)**3)
137     ! Write some debugging information if needed.
138                 if(dFlag) write(*,1100) M, VEL_TMP(M), NPMpSEC(M)
139              ENDDO
140     
141     ! Total number of particles at BCV injected per second
142              NPpSEC = sum(NPMpSEC)
143     
144      1100 FORMAT(/2x,'Conversion Info: Phase ',I2,/4x,'Velocity: ',g12.5,/ &
145              4X,'NPMpSEC = ',F11.1)
146     
147              if(dFlag) write(*,"(/2x,'Max Velocity:',3x,g12.5)") MAX_VEL
148              if(dFlag) write(*,"( 2x,'NPpSEC:',3x,F11.1)") NPpSEC
149     
150     ! The number of total particles per solid time step DTSOLID
151              NPpDT = NPpSEC * DTSOLID
152     
153     ! Inject one particle every PI_FACTOR solids time steps.
154              IF(NPpDT .LT. 1)THEN
155                 PI_COUNT(BCV_I) = 1
156                 PI_FACTOR(BCV_I) = FLOOR(real(1.d0/NPpDT))
157     ! Inject PI_COUNT particles every soilds time step.
158              ELSE
159                 PI_COUNT(BCV_I) = CEILING(real(NPpDT))
160                 PI_FACTOR(BCV_I) = 1
161              ENDIF
162     
163              OCCUPANTS = DEM_MI(BCV_I)%OCCUPANTS
164     
165     ! Calculate the minimum inlet velocity. The cutoff is associated with
166     ! square packing of disks on a plane.
167              MINIPV = MAX_DIA/( DTSOLID*dble(PI_FACTOR(BCV_I)) * dble(     &
168                 FLOOR( real(OCCUPANTS)/real(PI_COUNT(BCV_I)))))
169     ! Calculate the velocity needed to ensure that half the inlet is free.
170     ! Inlets with velocities greater than this value can be randomly seeded,
171     ! otherwise, particles are seeded in according to the grid.
172              MAXIPV = MAX_DIA/( DTSOLID*dble(PI_FACTOR(BCV_I)) * dble(     &
173                 FLOOR(CEILING(real(OCCUPANTS)/2.0)/real(PI_COUNT(BCV_I)))))
174     
175              if(dFlag) write(*,"(/2x,'MaxIPV:',3x,g12.5)") MAXIPV
176              if(dFlag) write(*,"( 2x,'MinIPV:',3x,g12.5)") MINIPV
177     
178              IF(MAX_VEL < MINIPV) THEN
179                 WRITE(ERR_MSG,1110) BCV, MAX_VEL, MINIPV
180                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
181              ENDIF
182     
183      1110 FORMAT('Error 1110: Solids velocity for BC ',I3,' is too low ', &
184              'to satisfy DEM',/'inlet packing restrictions. Potential ',  &
185              'solutions:',//,' > If the BC velocities (BC_U_s, BC_V_s, ', &
186              'BC_W_s) are defined, specify',/3x,'a larger value for the ',&
187              'velocity normal to the flow plane.',//' > If MASSFLOW or ', &
188              'VOLFLOW are defined, decrease the solids volume',/3x,       &
189              'fraction to increase solids velocity.',//2x,'Max user-',    &
190              'specified BC velocity:   ',g12.5,/2x,'Minimum required ',   &
191              'solids Velocity: ',g12.5)
192     
193     ! Set all BC solids velocities to the largest velocity and recalculate
194     ! BC_EP_s to determine the magnitude of the change.
195              EPs_ERR = ZERO
196              EPs_TMP = ZERO
197              DO MM = 1, PHASE_CNT
198                 M = PHASE_LIST(MM)
199                 EPs_TMP(M) = BC_EP_s(BCV,M+SMAX) * (VEL_TMP(M) / MAX_VEL)
200                 EPs_ERR = EPs_ERR + (BC_EP_s(BCV,M+SMAX) - EPs_TMP(M))
201     
202     ! Over-write the current BC value.
203                 SELECT CASE(BC_PLANE(BCV))
204                 CASE('N'); BC_V_s(BCV,M) =  abs(MAX_VEL)
205                 CASE('S'); BC_V_s(BCV,M) = -abs(MAX_VEL)
206                 CASE('E'); BC_U_s(BCV,M) =  abs(MAX_VEL)
207                 CASE('W'); BC_U_s(BCV,M) = -abs(MAX_VEL)
208                 CASE('T'); BC_W_s(BCV,M) =  abs(MAX_VEL)
209                 CASE('B'); BC_W_s(BCV,M) = -abs(MAX_VEL)
210                 END SELECT
211     
212              ENDDO
213     
214     ! If the net change in solids volume fraction is greatere than 0.01,
215     ! flag this as an error and exit. >> Let the user fix the input.
216              IF(.NOT.COMPARE(EPs_ERR,SMALL_NUMBER)) THEN
217                 IF(EPs_ERR > 0.01) THEN
218                    WRITE(ERR_MSG,1200) BCV, MAX_VEL, EPs_ERR
219                    CALL FLUSH_ERR_MSG(FOOTER=.FALSE.)
220                    FATAL = .TRUE.
221     
222     ! Report the amount of changes imposed on the BC in setting a
223     ! uniform inlet velocity.
224                 ELSE
225                    WRITE(ERR_MSG,1205) BCV, MAX_VEL, EPs_ERR
226                    CALL FLUSH_ERR_MSG(FOOTER=.FALSE.)
227                    FATAL = .FALSE.
228                 ENDIF
229     
230                 WRITE(ERR_MSG, 1210)
231                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
232                 DO MM = 1, PHASE_CNT
233                 M = PHASE_LIST(MM)
234                    WRITE(ERR_MSG,1211) M, BC_EP_s(BCV,M), EPs_TMP(M), &
235                       (BC_EP_s(BCV,M)-EPs_TMP(M))
236                    CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
237                 ENDDO
238                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., ABORT=FATAL)
239              ENDIF
240     
241     
242      1200 FORMAT('Error 1200: Unable to impose a uniform solids velocity ',&
243              'on BC',I3,'.',/'Setting all solids to the highest velocity ',&
244              'results in a large error',/'in the solids volume fraction. ',&
245              'Please correct the mfix.dat file.',2/5X,'Max Inlet Velocity',&
246              ':',1X,ES11.4,/5x,'Total BC_EP_s Error:',ES11.4)
247     
248     
249      1205 FORMAT('Warning 1201: Uniform solids velocity imposed on BC',I3, &
250              '.',2/,5X,'Uniform Inlet Velocity:',1X,ES11.4,/5X,'Total ',   &
251              'BC_EP_s Error:',4X,ES11.4,/' ')
252     
253      1210 FORMAT(/,5X,'|',11('-'),'|',3(14('-'),'|'),/5X,'|',3X,'Phase',3X,&
254              '|',4X,'BC_EP_s',3X,'|',2X,'Calculated',2X,'|',3X,'ABS ',   &
255              'Error',2X,'|',/5X,'|',11('-'),'|',3(14('-'),'|'))
256     
257      1211 FORMAT(5X,'|',4X,I2,5X,'| ',1X,ES11.4,1X,'|',1X,ES11.4,2X,'|',   &
258              1X,ES11.4,1X,' |',/5X,'|',11('-'),'|',3(14('-'),'|'))
259     
260     
261     ! For polydisperse inlets, construct the DES_POLY_LAYOUT array
262              IF(DEM_MI(BCV_I)%POLYDISPERSE) THEN
263                 RANGE_BOT = 1
264                 DO MM=1,PHASE_CNT - 1
265                    M = PHASE_LIST(MM)
266                    SCALED_VAL = dble(NUMFRAC_LIMIT)*(NPMpSEC(M)/NPpSEC)
267                    RANGE_TOP = FLOOR(SCALED_VAL) + (RANGE_BOT-1)
268                    DEM_BC_POLY_LAYOUT(BCV_I,RANGE_BOT:RANGE_TOP) = M
269                    RANGE_BOT = RANGE_TOP+1
270                 ENDDO
271     
272                 M = PHASE_LIST(PHASE_CNT)
273                 DEM_BC_POLY_LAYOUT(BCV_I,RANGE_BOT:NUMFRAC_LIMIT) = M
274     ! For monodisperse inlets, store the single mass phase used
275              ELSE
276                 DEM_BC_POLY_LAYOUT(BCV_I,:) = PHASE_LIST(1)
277              ENDIF
278     
279     
280     ! Calculate des mass inlet time; time between injection.  If the run
281     ! type is RESTART_1, DES_MI_TIME will be picked up from the restart file
282     ! with an updated value.
283              IF(RUN_TYPE /= 'RESTART_1') &
284                 DEM_MI_TIME(BCV_I) = TIME + dble(PI_FACTOR(BCV_I)) * DTSOLID
285     
286     
287              WRITE(ERR_MSG,1000) BCV, NPpDT, int(NPpDT/DTSOLID), &
288                 PI_FACTOR(BCV_I), PI_COUNT(BCV_I), DEM_MI_TIME(BCV_I)
289              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
290     
291      1000 FORMAT(2/,2X,'For mass inlet BC: ', I3,/,&
292              4X,'No. particles injected per solids time step = ', ES15.8,/,&
293              4X,'No. particles injected per second = ', I10,/,&
294              4X,'PI_FACTOR = ', I10,' PI_COUNT = ', I5,/,&
295              4X,'start DES_MI_TIME = ', ES15.8,/'  ')
296     
297           ENDDO
298     
299     !
300           CALL SET_DEM_BCMI_IJK
301     
302     
303           CALL FINL_ERR_MSG
304     
305     
306           RETURN
307           END SUBROUTINE SET_BC_DEM_MI
308     
309     
310     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
311     !                                                                      !
312     !  Subroutine: SET_BC_DEM                                              !
313     !  Author: J.Musser                                   Date: 13-Jul-09  !
314     !                                                                      !
315     !  Purpose: Check the data provided for the des mass inflow boundary   !
316     !  condition and flag errors if the data is improper.  This module is  !
317     !  also used to convert the proveded information into the format       !
318     !  necessary for the dependent subrountines to function properly.      !
319     !                                                                      !
320     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
321           SUBROUTINE SET_DEM_BCMI_IJK
322     
323           use bc, only: BC_PLANE
324           use bc, only: BC_X_w, BC_X_e, BC_Y_s, BC_Y_n, BC_Z_b, BC_Z_t
325     
326           use des_bc, only: DEM_BCMI, DEM_BCMI_MAP, DEM_BCMI_IJK
327           use des_bc, only: DEM_BCMI_IJKSTART, DEM_BCMI_IJKEND
328     
329           use funits, only: DMP_LOG
330     
331           use mpi_utility
332           use error_manager
333           use desgrid, only: dg_IMIN1, dg_IMAX1
334           use desgrid, only: dg_JMIN1, dg_JMAX1
335           use desgrid, only: dg_KMIN1, dg_KMAX1
336     
337           use desgrid, only: DG_FUNIJK
338           use desgrid, only: IofPOS, JofPOS, KofPOS
339           use desgrid, only: dg_is_ON_myPE_plus1layers
340     
341           IMPLICIT NONE
342     
343           INTEGER, ALLOCATABLE :: LOC_DEM_BCMI_IJK(:)
344     
345           INTEGER :: BCV, BCV_I
346     
347           INTEGER :: LC
348     
349           INTEGER :: MAX_CELLS
350     
351           INTEGER :: BND1, BND2
352     
353           LOGICAL, parameter :: setDBG = .FALSE.
354           LOGICAL :: dFlag
355     
356           INTEGER :: I,J,K,IJK
357           INTEGER :: I_w, I_e, J_s, J_n, K_b, K_t
358     
359           CALL INIT_ERR_MSG("SET_DEM_BCMI_IJK")
360     
361           dFlag = (DMP_LOG .AND. setDBG)
362     
363           if(dFlag) write(*,"(2/,2x,'From: SET_DEM_BCMI_IJK')")
364     
365     
366     
367     ! Loop over all inflow BCs to get an approximate count of the number
368     ! of fluid cells that are adjacent to them.
369           MAX_CELLS = 0
370           DO BCV_I=1, DEM_BCMI
371              BCV = DEM_BCMI_MAP(BCV_I)
372     
373     ! Set the search area a little bigger than the inlet area.
374              if(dFlag) WRITE(*,"(/2x,'Adding cells for BCV: ',I3)") BCV
375              SELECT CASE (BC_PLANE(BCV))
376              CASE('N','S')
377                 BND1 = min(IofPOS(BC_X_e(BCV))+1,dg_IMAX1) - &
378                    max(IofPOS(BC_X_w(BCV))-1,dg_IMIN1)
379                 BND2 = min(KofPOS(BC_Z_t(BCV))+1,dg_KMAX1) - &
380                    max(KofPOS(BC_Z_b(BCV))-1,dg_KMIN1)
381     
382              CASE('E','W')
383                 BND1 = min(JofPOS(BC_Y_n(BCV))+1,dg_JMAX1) - &
384                    max(JofPOS(BC_Y_s(BCV))-1,dg_JMIN1)
385                 BND2 = min(KofPOS(BC_Z_t(BCV))+1,dg_KMAX1) - &
386                    max(KofPOS(BC_Z_b(BCV))-1,dg_KMIN1)
387     
388              CASE('T','B')
389                 BND1 = min(IofPOS(BC_X_e(BCV))+1,dg_IMAX1) - &
390                    max(IofPOS(BC_X_w(BCV))-1,dg_IMIN1)
391                 BND2 = min(JofPOS(BC_Y_n(BCV))+1,dg_JMAX1) - &
392                    max(JofPOS(BC_Y_s(BCV))-1,dg_JMIN1)
393              END SELECT
394     
395              MAX_CELLS = MAX_CELLS + (BND1 + 1)*(BND2 + 1)
396              if(dFlag) WRITE(*,"(4x,'Plane:   ',A)") BC_PLANE(BCV)
397              if(dFlag) WRITE(*,"(4x,'Cells: ', I8)") (BND1+1)*(BND2+1)
398           ENDDO
399     
400           if(dFlag) write(*,"(2x,'Max Cells: ',I8)") MAX_CELLS
401     
402     ! Allocate an array to hold the IJK values. This array should be
403     ! more than enough to store all the IJKs.
404           allocate( LOC_DEM_BCMI_IJK(MAX_CELLS) )
405     
406     
407     ! Loop over the IJKs for each BC and store only the IJKs that you
408     ! own as well as the start/end array positions for each BC.
409           LC = 1
410           DO BCV_I = 1, DEM_BCMI
411     
412              DEM_BCMI_IJKSTART(BCV_I) = LC
413              BCV = DEM_BCMI_MAP(BCV_I)
414     
415              if(dFlag) write(*,"(/2x,'Searching for fluid cells:',I3)") BCV
416     
417              I_w = max(IofPOS(BC_X_w(BCV))-1,dg_IMIN1)
418              I_e = min(IofPOS(BC_X_e(BCV))+1,dg_IMAX1)
419              J_s = max(JofPOS(BC_Y_s(BCV))-1,dg_JMIN1)
420              J_n = min(JofPOS(BC_Y_n(BCV))+1,dg_JMAX1)
421              K_b = max(KofPOS(BC_Z_b(BCV))-1,dg_KMIN1)
422              K_t = min(KofPOS(BC_Z_t(BCV))+1,dg_KMAX1)
423     
424     ! Depending on the flow plane, the 'common' index needs set to reference
425     ! the fluid cell.
426              SELECT CASE (BC_PLANE(BCV))
427              CASE('N')
428                 J_s = JofPOS(BC_Y_s(BCV))+1
429                 J_n = JofPOS(BC_Y_n(BCV))+1
430              CASE('S')
431                 J_s = JofPOS(BC_Y_s(BCV))-1
432                 J_n = JofPOS(BC_Y_n(BCV))-1
433              CASE('E')
434                 I_w = IofPOS(BC_X_w(BCV))+1
435                 I_e = IofPOS(BC_X_e(BCV))+1
436              CASE('W')
437                 I_w = IofPOS(BC_X_w(BCV))-1
438                 I_e = IofPOS(BC_X_e(BCV))-1
439              CASE('T')
440                 K_b = KofPOS(BC_Z_b(BCV))+1
441                 K_t = KofPOS(BC_Z_t(BCV))+1
442              CASE('B')
443                 K_b = KofPOS(BC_Z_b(BCV))-1
444                 K_t = KofPOS(BC_Z_t(BCV))-1
445              END SELECT
446     
447              if(dFlag) then
448                 write(*,"(4x,'Search bounds: ')")
449                 write(*,"(6x,'I_w/I_e:',2(2x,I6))") I_w, I_e
450                 write(*,"(6x,'J_s/J_n:',2(2x,I6))") J_s, J_n
451                 write(*,"(6x,'K_b/K_t:',2(2x,I6))") K_b, K_t
452              endif
453     
454     ! Store the IJKs.
455              DO K = K_b, K_t
456              DO J = J_s, J_n
457              DO I = I_w, I_e
458     ! Skip cells that this rank does not own or are considered dead.
459                 IF(.NOT.dg_is_ON_myPE_plus1layers(I,J,K))CYCLE
460     
461                 IJK = DG_FUNIJK(I,J,K)
462                 LOC_DEM_BCMI_IJK(LC) = IJK
463                 LC = LC+1
464              ENDDO
465              ENDDO
466              ENDDO
467     
468              DEM_BCMI_IJKEND(BCV_I) = LC-1
469     
470              IF(dFLAG) write(*,1111) BCV, BCV_I,                           &
471                 DEM_BCMI_IJKSTART(BCV_I), DEM_BCMI_IJKEND(BCV_I)
472     
473           ENDDO
474     
475      1111 FORMAT(/2x,'DEM Mass Inflow:',/4x,'BC:',I4,3x,'MAP:',I4,         &
476              /4x,'IJKSTART:',I6,/4x,'IJKEND:  ',I6)
477     
478     
479     ! Allocate the global store arrary array. This changes across MPI ranks.
480           IF(LC > 1) THEN
481              allocate( DEM_BCMI_IJK(LC-1) )
482              DEM_BCMI_IJK(1:LC-1) = LOC_DEM_BCMI_IJK(1:LC-1)
483           ELSE
484              allocate( DEM_BCMI_IJK(1) )
485              DEM_BCMI_IJK(1) = LOC_DEM_BCMI_IJK(1)
486           ENDIF
487     
488           deallocate(LOC_DEM_BCMI_IJK)
489     
490     
491           CALL FINL_ERR_MSG
492     
493     
494           RETURN
495           END SUBROUTINE SET_DEM_BCMI_IJK
496