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