File: RELATIVE:/../../../mfix.git/model/set_bc0.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: set_bc0                                                 C
4     !  Purpose: This subroutine does the initial setting of all boundary   C
5     !  conditions. The user specifications of the boundary conditions are  C
6     !  checked for veracity in various check_data/ routines:               C
7     !  (e.g., check_boundary_conditions).                                  C
8     !                                                                      C
9     !  Author: M. Syamlal                                 Date: 29-JAN-92  C
10     !                                                                      C
11     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
12           SUBROUTINE SET_BC0
13     
14     ! Modules
15     !--------------------------------------------------------------------//
16           use bc, only: ijk_p_g
17           use bc, only: bc_type, bc_defined
18           use fldvar, only: x_g, t_g, p_g
19           use mms, only: calculate_mms, calculate_mms_source, use_mms
20     
21           use param, only: dimension_bc
22           use param1, only: undefined_i
23     
24           use mpi_utility
25           use sendrecv
26           IMPLICIT NONE
27     
28     ! Local variables
29     !--------------------------------------------------------------------//
30     ! Local index for boundary condition
31           INTEGER ::  L
32     !--------------------------------------------------------------------//
33     
34     ! Incompressible cases require that Ppg specified for one cell.
35     ! The following attempts to pick an appropriate cell.
36           CALL SET_IJK_P_G
37     
38           IF(USE_MMS) THEN
39     ! IJK_P_G is set as UNDEFINED for MMS since current IJK_P_G setting
40     ! is not a second order operation.
41              IJK_P_G = UNDEFINED_I
42     ! Calculate MMS variables. Better place might be inside interate before
43     ! every time-step. (?)
44              CALL CALCULATE_MMS
45              CALL CALCULATE_MMS_SOURCE
46           END IF
47     
48           DO L = 1, DIMENSION_BC
49              IF (BC_DEFINED(L)) THEN
50     
51                 SELECT CASE (TRIM(BC_TYPE(L)))
52     
53                 CASE ('FREE_SLIP_WALL')
54                    CALL SET_BC0_WALLS(L)
55                 CASE ('NO_SLIP_WALL')
56                    CALL SET_BC0_WALLS(L)
57                 CASE ('PAR_SLIP_WALL')
58                    CALL SET_BC0_WALLS(L)
59                 CASE ('P_OUTFLOW')
60                    CALL SET_BC0_INIT_BCDT_CALCS(L)
61                    CALL SET_BC0_OUTFLOW(L)
62                    CALL SET_BC0_INIT_BCSCALARS(L)
63                 CASE ('MASS_OUTFLOW')
64                    CALL SET_BC0_INIT_BCDT_CALCS(L)
65                    CALL SET_BC0_INFLOW(L)
66                    CALL SET_BC0_INIT_BCSCALARS(L)
67                 CASE ('OUTFLOW')
68                    CALL SET_BC0_INIT_BCDT_CALCS(L)
69                    CALL SET_BC0_OUTFLOW(L)
70                    CALL SET_BC0_INIT_BCSCALARS(L)
71                 CASE ('MASS_INFLOW')
72                    CALL SET_BC0_INIT_JET(L)
73                    CALL SET_BC0_INFLOW(L)
74                 CASE ('P_INFLOW')
75                    CALL SET_BC0_INFLOW(L)
76                 END SELECT
77              ENDIF
78           ENDDO
79     
80     ! Make T_g nonzero in k=0,1 ghost layers when k-decomposition employed
81           call send_recv(T_G,2)
82           call send_recv(P_G,2)
83           call send_recv(X_G,2)
84     
85           RETURN
86           END SUBROUTINE SET_BC0
87     
88     
89     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
90     !                                                                      C
91     !  Subroutine: set_bc0_walls                                           C
92     !  Purpose: Define certain field variables at the wall boundaries      C
93     !  according to the user specifications. These are not the real        C
94     !  values in the wall cells, only initial guesses.                     C
95     !                                                                      C
96     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
97           SUBROUTINE SET_BC0_WALLS(BCV)
98     
99     ! Modules
100     !--------------------------------------------------------------------//
101           use bc, only: bc_k_b, bc_k_t
102           use bc, only: bc_j_s, bc_j_n
103           use bc, only: bc_i_w, bc_i_e
104           use bc, only: bc_tw_g, bc_tw_s, bc_thetaw_m
105           use bc, only: bc_xw_g, bc_xw_s, bc_scalarw
106     
107           use fldvar, only: x_g, x_s, scalar
108           use fldvar, only: t_g, t_s, theta_m
109     
110           use param1, only: undefined
111           use physprop, only: smax, mmax, nmax
112           use scalars, only: nscalar
113     
114           use functions, only: is_on_mype_plus2layers
115           use boundfunijk, only: bound_funijk
116           use compar, only: dead_cell_at
117           IMPLICIT NONE
118     
119     ! Dummy arguments
120     !--------------------------------------------------------------------//
121     ! index of boundary
122           INTEGER, INTENT(IN) :: BCV
123     
124     ! Local variables
125     !--------------------------------------------------------------------//
126     ! indices
127           INTEGER :: I, J, K, IJK, M
128     !--------------------------------------------------------------------//
129     
130           DO K = BC_K_B(BCV), BC_K_T(BCV)
131           DO J = BC_J_S(BCV), BC_J_N(BCV)
132           DO I = BC_I_W(BCV), BC_I_E(BCV)
133     
134              IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
135              IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
136              IJK = BOUND_FUNIJK(I,J,K)
137     
138              IF(BC_Tw_g(BCV) /= UNDEFINED) T_g(IJK) = BC_Tw_g(BCV)
139     
140              IF(NMAX(0) > 0) &
141                 WHERE (BC_Xw_G(BCV,:NMAX(0)) /= UNDEFINED) &
142                 X_G(IJK,:NMAX(0)) = BC_Xw_G(BCV,:NMAX(0))
143     
144              IF(SMAX > 0) &
145                 WHERE(BC_Tw_s(BCV,:SMAX) /= UNDEFINED) &
146                 T_s(IJK,:SMAX) = BC_Tw_s(BCV,:SMAX)
147     
148              IF(MMAX > 0) &
149                 WHERE(BC_Thetaw_m(BCV,:MMAX) /= UNDEFINED) &
150                 Theta_m(IJK,:MMAX) = BC_Thetaw_m(BCV,:MMAX)
151     
152              DO M = 1, SMAX
153                 IF(NMAX(M) > 0) &
154                    WHERE (BC_Xw_s(BCV,M,:NMAX(M)) /= UNDEFINED) &
155                    X_s(IJK,M,:NMAX(M)) = BC_Xw_s(BCV,M,:NMAX(M))
156              ENDDO
157     
158              IF(NScalar > 0) &
159                 WHERE (BC_ScalarW(BCV,:NScalar) /= UNDEFINED) &
160                 Scalar(IJK,:NScalar) = BC_ScalarW(BCV,:NScalar)
161     
162           ENDDO   ! end do (i=bc_i_w(l),bc_i_e(l))
163           ENDDO   ! end do (j=bc_j_s(l),bc_j_n(l))
164           ENDDO   ! end do (k=bc_k_b(l),bc_k_t(l))
165     
166           RETURN
167           END SUBROUTINE SET_BC0_WALLS
168     
169     
170     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
171     !                                                                      C
172     !  Subroutine: set_bc0_outflow                                         C
173     !  Purpose: Set the initial settings of the boundary conditions (if    C
174     !  they are defined) for pressure outflow (PO) or outflow (O)          C
175     !  boundary types.                                                     C
176     !                                                                      C
177     !  Comments: For a new run the field variables are undefined in the    C
178     !  boundary cell locations, while for a restart run the field variable C
179     !  may have an existing value based on the preceding simulation.       C
180     !  Regardless, a user defined BC value will supercede any existing     C
181     !  value.                                                              C
182     !                                                                      C
183     !  Note: Scalar field quantities (i.e., T_G, T_s, x_g, x_s, Theta_m,   C
184     !  scalar, k_turb_g, e_turb_g) do not (should not) need to be defined  C
185     !  in any type of outflow boundary as the boundary values are unused   C
186     !  by the corresponding matrix equation (bc_phi). The only reason      C
187     !  this is potentially necessary is during a calculation another cell  C
188     !  references the value of this boundary cell before the corresponding C
189     !  governing equation solver is called.                                C
190     !                                                                      C
191     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
192           SUBROUTINE SET_BC0_OUTFLOW(BCV)
193     
194     ! Modules
195     !--------------------------------------------------------------------//
196           use bc, only: bc_k_b, bc_k_t
197           use bc, only: bc_j_s, bc_j_n
198           use bc, only: bc_i_w, bc_i_e
199           use bc, only: bc_p_g
200           use bc, only: bc_t_g, bc_t_s, bc_theta_m
201           use bc, only: bc_x_g, bc_x_s, bc_scalar
202           use bc, only: bc_ep_g, bc_rop_s
203           use bc, only: bc_k_turb_g, bc_e_turb_g
204     
205           use constant, only: pi
206     
207           use fldvar, only: x_g, x_s, scalar
208           use fldvar, only: t_g, t_s, theta_m
209           use fldvar, only: p_g, p_star
210           use fldvar, only: k_turb_g, e_turb_g
211           use fldvar, only: ep_g, rop_s
212           use fldvar, only: d_p, ro_s
213     
214           use mms, only: use_mms
215           use mms, only: mms_p_g, mms_t_g, mms_ep_g
216           use mms, only: mms_theta_m, mms_t_s, mms_rop_s
217     
218           use param1, only: undefined, zero
219           use physprop, only: smax, mmax, nmax
220           use run, only: k_epsilon, kt_type_enum, ghd_2007
221           use scalars, only: nscalar
222           use scales, only: scale_pressure
223           use toleranc, only: tmin
224     
225           use functions, only: is_on_mype_plus2layers
226           use compar, only: dead_cell_at
227           use boundfunijk, only: bound_funijk
228           IMPLICIT NONE
229     
230     ! Dummy arguments
231     !--------------------------------------------------------------------//
232     ! index of boundary condition
233           INTEGER, INTENT(IN) :: BCV
234     
235     ! Local variables
236     !--------------------------------------------------------------------//
237     ! indices
238           INTEGER :: I, J, K, IJK, M
239     ! number densities for use in GHD theory only
240           DOUBLE PRECISION :: nM, nTOT
241     !--------------------------------------------------------------------//
242     
243     
244           DO K = BC_K_B(BCV), BC_K_T(BCV)
245           DO J = BC_J_S(BCV), BC_J_N(BCV)
246           DO I = BC_I_W(BCV), BC_I_E(BCV)
247              IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
248              IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
249              IJK = BOUND_FUNIJK(I,J,K)
250     
251              P_STAR(IJK) = ZERO
252              P_G(IJK) = SCALE_PRESSURE(BC_P_G(BCV))
253              IF (BC_EP_G(BCV) /= UNDEFINED) EP_G(IJK) = BC_EP_G(BCV)
254     
255              T_G(IJK)= merge(BC_T_G(BCV), TMIN,&
256                 BC_T_G(BCV) /= UNDEFINED)
257     
258              IF (NMAX(0) > 0) &
259                 WHERE (BC_X_G(BCV,:NMAX(0)) /= UNDEFINED) &
260                 X_G(IJK,:NMAX(0)) = BC_X_G(BCV,:NMAX(0))
261     
262              IF (NScalar > 0) &
263                 WHERE (BC_Scalar(BCV,:NScalar) /= UNDEFINED) &
264                 Scalar(IJK,:NScalar) = BC_Scalar(BCV,:NScalar)
265     
266              IF (K_Epsilon) THEN
267                 IF (BC_K_Turb_G(BCV) /= UNDEFINED) &
268                    K_Turb_G(IJK) = BC_K_Turb_G(BCV)
269                 IF (BC_E_Turb_G(BCV) /= UNDEFINED) &
270                     E_Turb_G(IJK) = BC_E_Turb_G(BCV)
271              ENDIF
272     
273              DO M = 1, SMAX
274                 IF (BC_ROP_S(BCV,M) /= UNDEFINED) &
275                    ROP_S(IJK,M) = BC_ROP_S(BCV,M)
276                 IF(BC_T_S(BCV,M) /= UNDEFINED) &
277                    T_S(IJK,M)=BC_T_S(BCV,M)
278                 IF (BC_THETA_M(BCV,M) /= UNDEFINED) &
279                    THETA_M(IJK,M) = BC_THETA_M(BCV,M)
280     
281                 IF (NMAX(M) > 0) &
282                    WHERE (BC_X_S(BCV,M,:NMAX(M)) /= UNDEFINED) &
283                    X_S(IJK,M,:NMAX(M)) = BC_X_S(BCV,M,:NMAX(M))
284              ENDDO
285     
286     ! for GHD theory to compute mixture BC of velocity and density
287              IF(KT_TYPE_ENUM == GHD_2007) THEN
288                 ROP_S(IJK,MMAX) = ZERO
289                 nTOT = ZERO
290                 THETA_M(IJK,MMAX) = ZERO
291                 DO M = 1, SMAX
292                    IF (BC_ROP_S(BCV,M) /= UNDEFINED) THEN
293                       ROP_S(IJK,MMAX) = ROP_S(IJK,MMAX) + &
294                          BC_ROP_S(BCV,M)
295                       nM = BC_ROP_S(BCV,M)*6d0 / &
296                          (PI*D_p(IJK,M)**3*RO_S(IJK,M))
297                       nTOT = nTOT + nM
298                       IF (BC_THETA_M(BCV,M) /= UNDEFINED) &
299                          THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) + &
300                          nM*BC_THETA_M(BCV,M)
301                    ENDIF
302                 ENDDO
303                 IF(ROP_S(IJK,MMAX) > ZERO) &
304                    THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) / nTOT
305              ENDIF   ! end if(kt_type_enum== ghd_2007)
306     
307     ! Set MMS BCs when PO boundary condition is used.
308              IF (USE_MMS) THEN
309                 P_G(IJK) = SCALE_PRESSURE(MMS_P_G(IJK))
310                 EP_G(IJK) = MMS_EP_G(IJK)
311                 T_G(IJK) = MMS_T_G(IJK)
312     
313                 M = 1 ! Single solid phase for MMS cases
314                 ROP_S(IJK,M) = MMS_ROP_S(IJK)
315                 T_S(IJK,M) = MMS_T_S(IJK)
316                    THETA_M(IJK,M) = MMS_THETA_M(IJK)
317              ENDIF ! end if(USE_MMS)
318     
319           ENDDO   ! do i
320           ENDDO   ! do j
321           ENDDO   ! do k
322     
323           RETURN
324           END SUBROUTINE SET_BC0_OUTFLOW
325     
326     
327     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
328     !                                                                      C
329     !  Subroutine: set_bc0_inflow                                          C
330     !  Purpose: Set the initial settings of the boundary conditions        C
331     !  pressure inflow (PI) and. mass inflow (MI) boundary types. Also do  C
332     !  mass outflow (MO) boundary types... due to velocity...              C
333     !                                                                      C
334     !  Comments: Unlike the treament of PO or O boundary types no checks   C
335     !  are made for these boundary types to determine whether the given    C
336     !  BC value is defined before it is assigned to the field variable.    C
337     !  However, the corresponding check routines generally ensure such BC  C
338     !  quantities are defined for MI or PI boundaries if they are needed   C
339     !  for the simulation.                                                 C
340     !                                                                      C
341     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
342           SUBROUTINE SET_BC0_INFLOW(BCV)
343     
344     ! Modules
345     !--------------------------------------------------------------------//
346           use bc, only: bc_plane
347           use bc, only: bc_k_b, bc_k_t
348           use bc, only: bc_j_s, bc_j_n
349           use bc, only: bc_i_w, bc_i_e
350           use bc, only: bc_u_g, bc_v_g, bc_w_g
351           use bc, only: bc_u_s, bc_v_s, bc_w_s
352           use bc, only: bc_p_g
353           use bc, only: bc_t_g, bc_t_s, bc_theta_m
354           use bc, only: bc_x_g, bc_x_s, bc_scalar
355           use bc, only: bc_ep_g, bc_rop_s
356           use bc, only: bc_k_turb_g, bc_e_turb_g
357     
358           use constant, only: pi
359     
360           use fldvar, only: u_g, v_g, w_g
361           use fldvar, only: u_s, v_s, w_s
362           use fldvar, only: x_g, x_s, scalar
363           use fldvar, only: t_g, t_s, theta_m
364           use fldvar, only: p_g, p_star
365           use fldvar, only: k_turb_g, e_turb_g
366           use fldvar, only: ep_g, rop_s
367           use fldvar, only: d_p, ro_s
368     
369           use mms, only: use_mms
370           use mms, only: mms_u_g, mms_v_g, mms_w_g
371           use mms, only: mms_u_s, mms_v_s, mms_w_s
372           use mms, only: mms_p_g, mms_t_g, mms_ep_g
373           use mms, only: mms_theta_m, mms_t_s, mms_rop_s
374     
375           use param1, only: zero
376           use physprop, only: smax, mmax, nmax
377           use run, only: k_epsilon, kt_type_enum, ghd_2007
378           use scalars, only: nscalar
379           use scales, only: scale_pressure
380     
381           use indices, only: im1, jm1, km1
382           use functions, only: is_on_mype_plus2layers
383           use boundfunijk, only: bound_funijk
384           use compar, only: dead_cell_at
385           IMPLICIT NONE
386     
387     ! Dummy arguments
388     !--------------------------------------------------------------------//
389     ! index for boundary condition
390           INTEGER, INTENT(IN) :: BCV
391     
392     ! Local variables
393     !--------------------------------------------------------------------//
394     ! indices
395           INTEGER :: I, J, K, IJK, M
396     ! ijk index for setting normal component of velocity
397           INTEGER :: FIJK
398     ! number densities for use in GHD theory only
399           DOUBLE PRECISION :: nM, nTOT
400     ! calculation for normal component of mixture velocity
401           DOUBLE PRECISION :: lvel_s
402     !--------------------------------------------------------------------//
403     
404           DO K = BC_K_B(BCV), BC_K_T(BCV)
405           DO J = BC_J_S(BCV), BC_J_N(BCV)
406           DO I = BC_I_W(BCV), BC_I_E(BCV)
407              IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
408              IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
409              IJK = BOUND_FUNIJK(I,J,K)
410     
411              P_STAR(IJK) = ZERO
412              P_G(IJK) = SCALE_PRESSURE(BC_P_G(BCV))
413              EP_G(IJK) = BC_EP_G(BCV)
414              T_G(IJK) = BC_T_G(BCV)
415     
416              IF (NMAX(0) > 0) &
417                X_G(IJK,:NMAX(0)) = BC_X_G(BCV,:NMAX(0))
418     
419              IF (NScalar > 0) &
420                Scalar(IJK,:NScalar) = BC_Scalar(BCV,:NScalar)
421     
422              IF (K_Epsilon) THEN
423                K_Turb_G(IJK) = BC_K_Turb_G(BCV)
424                E_Turb_G(IJK) = BC_E_Turb_G(BCV)
425              ENDIF
426     
427              DO M = 1, SMAX
428                ROP_S(IJK,M) = BC_ROP_S(BCV,M)
429                T_S(IJK,M) = BC_T_S(BCV,M)
430                THETA_M(IJK,M) = BC_THETA_M(BCV,M)
431                IF (NMAX(M) > 0) &
432                   X_S(IJK,M,:NMAX(M)) = BC_X_S(BCV,M,:NMAX(M))
433              ENDDO
434     
435              U_G(IJK) = BC_U_G(BCV)
436              V_G(IJK) = BC_V_G(BCV)
437              W_G(IJK) = BC_W_G(BCV)
438              IF (SMAX > 0) THEN
439                 U_S(IJK,:SMAX) = BC_U_S(BCV,:SMAX)
440                 V_S(IJK,:SMAX) = BC_V_S(BCV,:SMAX)
441                 W_S(IJK,:SMAX) = BC_W_S(BCV,:SMAX)
442              ENDIF
443     
444     ! When the boundary plane is located on the E, N, T side of the domain
445     ! (fluid cell is located w, s, b), set the component of velocity normal
446     ! to the boundary plane of the adjacent fluid cell
447              SELECT CASE (TRIM(BC_PLANE(BCV)))
448                 CASE ('W')
449                    FIJK = BOUND_FUNIJK(IM1(I),J,K)
450                    U_G(FIJK) = BC_U_G(BCV)
451                    IF (SMAX >0) U_S(FIJK,:SMAX) = BC_U_S(BCV,:SMAX)
452                 CASE ('S')
453                    FIJK = BOUND_FUNIJK(I,JM1(J),K)
454                    V_G(FIJK) = BC_V_G(BCV)
455                    IF(SMAX>0) V_S(FIJK,:SMAX) = BC_V_S(BCV,:SMAX)
456                 CASE ('B')
457                    FIJK = BOUND_FUNIJK(I,J,KM1(K))
458                    W_G(FIJK) = BC_W_G(BCV)
459                    IF (SMAX>0) W_S(FIJK,:SMAX) = BC_W_S(BCV,:SMAX)
460              END SELECT
461     
462     ! for GHD theory to compute mixture BC of velocity and density
463              IF(KT_TYPE_ENUM == GHD_2007) THEN
464                 ROP_S(IJK,MMAX) = ZERO
465                 nTOT = ZERO
466                 THETA_M(IJK,MMAX) = ZERO
467                 U_S(IJK,MMAX) =  ZERO
468                 V_S(IJK,MMAX) =  ZERO
469                 W_S(IJK,MMAX) =  ZERO
470                 lvel_s = zero
471                 DO M = 1, SMAX
472                    ROP_S(IJK,MMAX) = ROP_S(IJK,MMAX) + &
473                       BC_ROP_S(BCV,M)
474                    nM = BC_ROP_S(BCV,M)*6d0/ &
475                       (PI*D_p(IJK,M)**3*RO_S(IJK,M))
476                    nTOT = nTOT + nM
477                    THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) + &
478                       nM*BC_THETA_M(BCV,M)
479                    U_S(IJK,MMAX) = U_S(IJK,MMAX) + &
480                       BC_ROP_S(BCV,M)*BC_U_S(BCV,M)
481                    V_S(IJK,MMAX) = V_S(IJK,MMAX) + &
482                       BC_ROP_S(BCV,M)*BC_V_S(BCV,M)
483                    W_S(IJK,MMAX) = W_S(IJK,MMAX) + &
484                       BC_ROP_S(BCV,M)*BC_W_S(BCV,M)
485     ! set velocity component normal to plane in adjacent fluid cell
486                    SELECT CASE (TRIM(BC_PLANE(BCV)))
487                       CASE ('W')
488                          lvel_s = lvel_s + BC_ROP_S(BCV,M)*BC_U_S(BCV,M)
489                       CASE ('S')
490                          lvel_s = lvel_s + BC_ROP_S(BCV,M)*BC_V_S(BCV,M)
491                       CASE ('B')
492                          lvel_s = lvel_s + BC_ROP_S(BCV,M)*BC_W_S(BCV,M)
493                    END SELECT
494                 ENDDO
495                 IF(ROP_S(IJK,MMAX) > ZERO) THEN
496                    THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) / nTOT
497                    U_S(IJK,MMAX) = U_S(IJK,MMAX) / ROP_S(IJK,MMAX)
498                    V_S(IJK,MMAX) = V_S(IJK,MMAX) / ROP_S(IJK,MMAX)
499                    W_S(IJK,MMAX) = W_S(IJK,MMAX) / ROP_S(IJK,MMAX)
500                    lvel_s = lvel_s/ROP_S(IJK,MMAX)
501                 ENDIF
502                 SELECT CASE (TRIM(BC_PLANE(BCV)))
503                    CASE ('W')
504                      FIJK = BOUND_FUNIJK(IM1(I),J,K)
505                      U_S(FIJK,MMAX) = lvel_s
506                    CASE ('S')
507                      FIJK = BOUND_FUNIJK(I,JM1(J),K)
508                      V_S(FIJK,MMAX) = lvel_s
509                    CASE ('B')
510                      FIJK = BOUND_FUNIJK(I,J,KM1(K))
511                      W_S(FIJK,MMAX) = lvel_s
512                 END SELECT
513     
514              ENDIF  ! end if (trim(kt_type) ='ghd')
515     
516     ! Set MMS BCs when MI boundary condition is used.
517              IF (USE_MMS) THEN
518                 P_G(IJK) = SCALE_PRESSURE(MMS_P_G(IJK))
519                 EP_G(IJK) = MMS_EP_G(IJK)
520                 T_G(IJK) = MMS_T_G(IJK)
521     
522                 DO M = 1, SMAX
523                    ROP_S(IJK,M) = MMS_ROP_S(IJK)
524                    T_S(IJK,M) = MMS_T_S(IJK)
525                    THETA_M(IJK,M) = MMS_THETA_M(IJK)
526                 ENDDO
527     
528                 U_G(IJK) = MMS_U_G(IJK)
529                 V_G(IJK) = MMS_V_G(IJK)
530                 W_G(IJK) = MMS_W_G(IJK)
531                 IF (SMAX > 0) THEN
532                    U_S(IJK,:SMAX) = MMS_U_S(IJK)
533                    V_S(IJK,:SMAX) = MMS_V_S(IJK)
534                    W_S(IJK,:SMAX) = MMS_W_S(IJK)
535                 ENDIF
536     ! When the boundary plane is W, S, or B the normal component of velocity
537     ! needs to be set for both sides of the boundary cell.
538                 SELECT CASE (TRIM(BC_PLANE(BCV)))
539                    CASE ('W')
540                      FIJK = BOUND_FUNIJK(IM1(I),J,K)
541                      U_G(FIJK) = MMS_U_G(FIJK)
542                      IF(SMAX>0) U_S(FIJK,:SMAX) = MMS_U_S(FIJK)
543                    CASE ('S')
544                      FIJK = BOUND_FUNIJK(I,JM1(J),K)
545                      V_G(FIJK) = MMS_V_G(FIJK)
546                      IF(SMAX>0) V_S(FIJK,:SMAX) = MMS_V_S(FIJK)
547                    CASE ('B')
548                      FIJK = BOUND_FUNIJK(I,J,KM1(K))
549                      W_G(FIJK) = MMS_W_G(FIJK)
550                      IF(SMAX>0) W_S(FIJK,:SMAX) = MMS_W_S(FIJK)
551                 END SELECT
552              ENDIF ! end if(USE_MMS)
553     
554           ENDDO   ! do i
555           ENDDO   ! do j
556           ENDDO   ! do k
557     
558           RETURN
559           END SUBROUTINE SET_BC0_INFLOW
560     
561     
562     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
563     !                                                                      C
564     !  Subroutine: set_bc0_init_jet                                        C
565     !  Purpose: initializing time dependent jet conditions                 C
566     !                                                                      C
567     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
568           SUBROUTINE SET_BC0_INIT_JET(BCV)
569     
570     ! Modules
571     !--------------------------------------------------------------------//
572           use bc, only: bc_plane
573           use bc, only: bc_jet_g, bc_jet_g0
574           use bc, only: bc_dt_0, bc_time
575           use bc, only: bc_u_g, bc_v_g, bc_w_g
576           use param1, only: undefined
577           use run, only: time
578           IMPLICIT NONE
579     
580     ! Dummy arguments
581     !--------------------------------------------------------------------//
582     ! index of boundary
583           INTEGER, INTENT(IN) :: BCV
584     !--------------------------------------------------------------------//
585     
586           BC_JET_G(BCV) = UNDEFINED
587           IF (BC_DT_0(BCV) /= UNDEFINED) THEN
588              BC_TIME(BCV) = TIME + BC_DT_0(BCV)
589              BC_JET_G(BCV) = BC_JET_G0(BCV)
590              IF (BC_JET_G(BCV) /= UNDEFINED) THEN
591                 SELECT CASE (TRIM(BC_PLANE(BCV)))
592                 CASE ('W', 'E')
593                    BC_U_G(BCV) = BC_JET_G(BCV)
594                 CASE ('S', 'N')
595                    BC_V_G(BCV) = BC_JET_G(BCV)
596                 CASE ('B', 'T')
597                    BC_W_G(BCV) = BC_JET_G(BCV)
598                 END SELECT
599              ENDIF
600           ELSE
601              BC_TIME(BCV) = UNDEFINED
602           ENDIF
603           RETURN
604           END SUBROUTINE SET_BC0_INIT_JET
605     
606     
607     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
608     !                                                                      C
609     !  Subroutine: set_bc0_init_bcdt_calcs                                 C
610     !  Purpose: initializing time dependent outflow calculations for       C
611     !  modifying outflow conditions (MO type) or simple reporting          C
612     !  outflow conditions (PO or O types)                                  C
613     !                                                                      C
614     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
615           SUBROUTINE SET_BC0_INIT_BCDT_CALCS(BCV)
616     
617     ! Modules
618     !--------------------------------------------------------------------//
619           use bc, only: bc_dt_0, bc_time
620           use bc, only: bc_out_n
621           use bc, only: bc_mout_g, bc_mout_s
622           use bc, only: bc_vout_g, bc_vout_s
623           use physprop, only: smax
624           use param1, only: undefined, zero
625           use run, only: time
626           IMPLICIT NONE
627     ! Dummy arguments
628     !--------------------------------------------------------------------//
629     ! index of boundary
630           INTEGER, INTENT(IN) :: BCV
631     !--------------------------------------------------------------------//
632     
633     ! initializing for time dependent outflow reporting calculation
634           IF (BC_DT_0(BCV) /= UNDEFINED) THEN
635              BC_TIME(BCV) = TIME + BC_DT_0(BCV)
636              BC_OUT_N(BCV) = 0
637              BC_MOUT_G(BCV) = ZERO
638              BC_VOUT_G(BCV) = ZERO
639              IF (SMAX > 0) THEN
640                 BC_MOUT_S(BCV,:SMAX) = ZERO
641                 BC_VOUT_S(BCV,:SMAX) = ZERO
642              ENDIF
643           ELSE
644              BC_TIME(BCV) = UNDEFINED
645           ENDIF
646           RETURN
647           END SUBROUTINE SET_BC0_INIT_BCDT_CALCS
648     
649     
650     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
651     !                                                                      C
652     !  Subroutine: set_bc0_init_scalars                                    C
653     !  Purpose: Provide a mechanism to initialize certain field variables  C
654     !  in the boundary cell to the value of the neighboring fluid cell.    C
655     !  It is not clear to what extent initialization is necessary as the   C
656     !  governing matrix equation for each of the field variables is set    C
657     !  such that the bc cells are not actually solved but are set to the   C
658     !  value of the adjacent fluid cell(see bc_phi). However, there may be C
659     !  some instance in the code where the bc cells are referenced before  C
660     !  the code has called the governing equation.                         C
661     !  This subroutine illustrates how these particular variables should   C
662     !  not need (or use) bc values defined by the user!                    C
663     !                                                                      C
664     !  Comments: this call replaces some initializations that were         C
665     !  'inadvertently' occurring in set_bc1                                C
666     !                                                                      C
667     !  WARNING: this routine only serves to initialize field variables     C
668     !  whose governing solver routine invokes bc_phi! do not insert any    C
669     !  other initializations here as it is likely not appropriate!!        C
670     !                                                                      C
671     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
672           SUBROUTINE set_bc0_init_bcscalars(BCV)
673     
674     ! Modules
675     !--------------------------------------------------------------------//
676           use bc, only: bc_plane
677           use bc, only: bc_k_b, bc_k_t
678           use bc, only: bc_j_s, bc_j_n
679           use bc, only: bc_i_w, bc_i_e
680     
681           use fldvar, only: t_g, t_s, theta_m
682           use fldvar, only: x_g, x_s, scalar
683           use fldvar, only: k_turb_g, e_turb_g
684           use physprop, only: smax, mmax
685           use physprop, only: nmax
686           use run, only: k_epsilon, kt_type_enum, ghd_2007
687           use scalars, only: nscalar
688     
689           use indices, only: im1, ip1, jm1, jp1, km1, kp1
690           use functions, only: is_on_mype_plus2layers
691           use boundfunijk, only: bound_funijk
692           use compar, only: dead_cell_at
693           IMPLICIT NONE
694     
695     ! Dummy arguments
696     !--------------------------------------------------------------------//
697     ! index of boundary
698           INTEGER, INTENT(IN) :: BCV
699     
700     ! Local variables
701     !--------------------------------------------------------------------//
702     ! indices of current boundary cell
703           INTEGER :: IJK, I, J, K
704     ! index of neighboring fluid cell
705           INTEGER :: FIJK
706     ! local indices
707           INTEGER :: M
708     !--------------------------------------------------------------------//
709     
710           DO K = BC_K_B(BCV), BC_K_T(BCV)
711           DO J = BC_J_S(BCV), BC_J_N(BCV)
712           DO I = BC_I_W(BCV), BC_I_E(BCV)
713              IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
714              IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
715              IJK = BOUND_FUNIJK(I,J,K)
716     
717              SELECT CASE (TRIM(BC_PLANE(BCV)))
718                 CASE ('W')   ! fluid cell at west
719                    FIJK = BOUND_FUNIJK(IM1(I),J,K)
720                 CASE ('E')
721                    FIJK = BOUND_FUNIJK(IP1(I),J,K)
722                 CASE ('S')   ! fluid cell at south
723                    FIJK = BOUND_FUNIJK(I,JM1(J),K)
724                 CASE ('N')
725                    FIJK = BOUND_FUNIJK(I,JP1(J),K)
726                 CASE ('B')   ! fluid cell at bottom
727                    FIJK = BOUND_FUNIJK(I,J,KM1(K))
728                 CASE ('T')
729                    FIJK = BOUND_FUNIJK(I,J,KP1(K))
730              END SELECT
731     
732              T_G(IJK) = T_G(FIJK)
733              IF (NMAX(0) > 0) &
734                 X_G(IJK,:NMAX(0)) = X_G(FIJK,:NMAX(0))
735     
736     ! setting scalar quantities
737              IF (NScalar >0) &
738                 Scalar(IJK, :NScalar) = Scalar(FIJK, :NScalar)
739     
740     ! setting turbulence quantities
741              IF(K_Epsilon) THEN
742                 K_Turb_G(IJK) = K_Turb_G(FIJK)
743                 E_Turb_G(IJK) = E_Turb_G(FIJK)
744              ENDIF
745     
746     ! this should only loop of tfm phases
747              DO M = 1, SMAX
748                 T_S(IJK,M) = T_S(FIJK,M)
749                 THETA_M(IJK,M) =  THETA_M(FIJK,M)
750                 IF (NMAX(M) > 0) &
751                    X_S(IJK,M,:NMAX(M)) = X_S(FIJK,M,:NMAX(M))
752              ENDDO
753     
754              IF(KT_TYPE_ENUM == GHD_2007) &
755                  THETA_M(IJK,MMAX) =  THETA_M(FIJK,MMAX)
756     
757           ENDDO
758           ENDDO
759           ENDDO
760           RETURN
761           END SUBROUTINE set_bc0_init_bcscalars
762     
763     
764     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
765     !                                                                      !
766     !  Subroutine: SET_IJK_P_G                                             !
767     !  Purpose: Pick an appropriate control volume to specify Ppg.         !
768     !                                                                      !
769     !  Author: J. Musser                                  Date: 07-Nov-13  !
770     !  Reviewer:                                          Date:            !
771     !                                                                      !
772     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
773           SUBROUTINE SET_IJK_P_G
774     
775     ! IJK location where Ppg is fixed.
776           use bc, only: IJK_P_G
777     ! Specified constant gas density.
778           use physprop, only: RO_G0
779     
780           use geometry, only: CYCLIC_X, CYCLIC_X_PD, CYCLIC_X_MF
781           use geometry, only: CYCLIC_Y, CYCLIC_Y_PD, CYCLIC_Y_MF
782           use geometry, only: CYCLIC_Z, CYCLIC_Z_PD, CYCLIC_Z_MF
783           use geometry, only: iMAX1, iMin1
784           use geometry, only: jMAX1, jMin1
785           use geometry, only: kMAX1, kMin1
786     
787           use geometry, only: do_K
788     
789           use funits, only: DMP_LOG
790     
791           use bc, only: BC_DEFINED
792           use bc, only: BC_TYPE
793     
794     ! MFIX Runtime parameters:
795           use param, only: DIMENSION_BC
796           use param1, only: UNDEFINED
797           use param1, only: UNDEFINED_I
798     
799           use mpi_utility
800     
801           implicit none
802     !--------------------------------------------------------------------//
803           INTEGER :: BCV
804     
805           CHARACTER(len=7) :: Map
806           CHARACTER(len=128) :: lMsg
807     
808           INTEGER :: l3
809           INTEGER :: l2, u2
810           INTEGER :: l1, u1
811     
812           LOGICAL, parameter :: setDBG = .FALSE.
813           LOGICAL :: dFlag
814           INTEGER :: iErr
815     
816     !--------------------------------------------------------------------//
817     
818           dFlag = (DMP_LOG .AND. setDBG)
819     
820     ! Initialize.
821           iErr = 0
822           IJK_P_G = UNDEFINED_I
823     
824     ! This is not needed for compressible cases.
825           IF(RO_G0 == UNDEFINED) THEN
826              IF(dFlag) write(*,"(3x,A)")                                   &
827                 'Compressible: IJK_P_g remaining undefined.'
828              return
829           ELSEIF(RO_G0 == 0.0d0) THEN
830              IF(dFlag) write(*,"(3x,A)")                                   &
831                 'No gas phase: IJK_P_g remaining undefined.'
832              return
833           ENDIF
834     
835     ! If there are no cyclic boundaries, look for a pressure outflow.
836           lpBCV: DO BCV = 1, DIMENSION_BC
837              IF (.NOT.BC_DEFINED(BCV)) cycle lpBCV
838              IF (BC_TYPE(BCV) == 'P_OUTFLOW' .OR. &
839                  BC_TYPE(BCV) == 'P_INFLOW') THEN
840                 IF(dFlag) write(*,"(3x,A)")                                &
841                    'Outflow PC defined: IJK_P_g remaining undefined.'
842                 RETURN
843              ENDIF
844           ENDDO lpBCV
845     
846     ! Initialize.
847              l3 = UNDEFINED_I
848              l2 = UNDEFINED_I; u2=l2
849              l1 = UNDEFINED_I; u1=l1
850     
851     ! If there are cyclic boundaries, flag a cell along the positive
852     ! domain extreme in the cyclic direction (e.g., JMAX1).
853           IF(CYCLIC_Y .OR. CYCLIC_Y_PD .OR. CYCLIC_Y_MF) THEN
854     
855              Map = 'JKI_MAP'
856              l3 = JMAX1
857              l2 = KMIN1;  u2 = KMAX1
858              l1 = IMIN1;  u1 = IMAX1
859              lMsg='Cyclic in Y'
860     
861           ELSEIF(CYCLIC_X .OR. CYCLIC_X_PD .OR. CYCLIC_X_MF) THEN
862     
863              Map = 'IKJ_MAP'
864              l3 = IMAX1
865              l2 = KMIN1;  u2 = KMAX1
866              l1 = JMIN1;  u1 = JMAX1
867              lMsg='Cyclic in X'
868     
869           ELSEIF(CYCLIC_Z .OR. CYCLIC_Z_PD .OR. CYCLIC_Z_MF) THEN
870     
871              Map = 'KIJ_MAP'
872              l3 = KMAX1
873              l2 = IMIN1;  u2 = IMAX1
874              l1 = JMIN1;  u1 = JMAX1
875              lMsg='Cyclic in Z'
876     
877           ENDIF
878     
879     ! No cyclic boundaries or pressure outflows. The IJ plane is used in
880     ! this case to maximize search region for 2D problems.
881           IF(l3 == UNDEFINED_I) THEN
882              Map = 'KIJ_MAP'
883              l3 = merge(max((KMAX1-KMIN1)/2+1,2), KMIN1, do_K)
884              l2 = IMIN1;  u2 = IMAX1
885              l1 = JMIN1;  u1 = JMAX1
886              lMsg='Center of domain'
887           ENDIF
888     
889     ! Debugging messages.
890           IF(dFlag) THEN
891              write(*,"(3/,3x,'Map: ',A)") Map
892              write(*,"(/5x,'l3:',2x,I4)") l3
893              write(*,"( 5x,'l2:',2(2x,I4))") l2, u2
894              write(*,"( 5x,'l1:',2(2x,I4))") l1, u1
895              write(*,"( 5x,'Msg: ',A)") trim(lMsg)
896           ENDIF
897     
898     ! Invoke the search routine.
899           CALL IJK_Pg_SEARCH(l3, l2, u2, l1, u1, MAP, dFlag, iErr)
900     
901           IF(iErr == 0) RETURN
902     
903     ! Error management.
904           IF(DMP_LOG) THEN
905              SELECT CASE (iErr)
906              CASE ( 1001);  WRITE(UNIT_LOG, 1001); WRITE(*,1001)
907              CASE ( 2000);  WRITE(UNIT_LOG, 2000); WRITE(*,2000)
908              CASE ( 2001);  WRITE(UNIT_LOG, 2001); WRITE(*,2001)
909              CASE ( 2002);  WRITE(UNIT_LOG, 2002); WRITE(*,2002)
910              CASE DEFAULT
911                 WRITE(UNIT_LOG, 1000) iErr
912                 WRITE(*,1000) iErr
913              END SELECT
914     
915              WRITE(UNIT_LOG, 9000) MAP(1:1), l3, MAP(2:2),                 &
916                 l2, u2, MAP(3:3), l1, u1
917              WRITE(*, 9000) MAP(1:1), l3, MAP(2:2),                        &
918                 l2, u2, MAP(3:3), l1, u1
919     
920              WRITE(*, 9999)
921              WRITE(UNIT_LOG, 9999)
922     
923           ENDIF
924     
925     
926           CALL MFIX_EXIT(myPE)
927     
928     
929      1000 FORMAT(//1X,70('*')/' From: SET_IJK_Pg',/,                       &
930              ' Error 1000: Unknown error reported. x', I4.4)
931     
932      1001 FORMAT(//1X,70('*')/' From: SET_IJK_Pg',/,                       &
933              ' Error 1001: Invalid mapping function.')
934     
935      2000 FORMAT(//1X,70('*')/' From: SET_IJK_Pg > IJK_Pg_SEARCH',/,       &
936              ' Error 2000: Unknown error reported from IJK_Pg_SEARCH.')
937     
938      2001 FORMAT(//1X,70('*')/' From: SET_IJK_Pg > IJK_Pg_SEARCH',/,       &
939              ' Error 2001: Unable to locate fluid cell in search region.')
940     
941      2002 FORMAT(//1X,70('*')/' From: SET_IJK_Pg > IJK_Pg_SEARCH',/,       &
942              ' Error 2002: Unable to locate fluid cell owner.')
943     
944      9000 FORMAT(/' Search plane information:',/,3x,A1,': ',I8,            &
945               2(/3x,A1,': ',I8,' x ',I8))
946     
947      9999 FORMAT(/' Fatal Error --> Invoking MFIX_EXIT',/1x,70('*'),2/)
948     
949           END SUBROUTINE SET_IJK_P_G
950     
951     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
952     !                                                                      !
953     !  Author: J. Musser                                  Date: 07-Nov-13  !
954     !  Reviewer:                                          Date:            !
955     !                                                                      !
956     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
957           SUBROUTINE IJK_Pg_SEARCH(ll3, ll2, lu2, ll1, lu1, lMAP,          &
958              ldFlag, iErr)
959     
960     ! Modules
961     !--------------------------------------------------------------------//
962     ! IJK location where Ppg is fixed.
963           use bc, only: IJK_P_g
964           use indices
965           use param1, only: undefined_i
966           use mpi_utility
967           use functions
968           implicit none
969     
970     ! Dummy arguments
971     !--------------------------------------------------------------------//
972           INTEGER, INTENT(IN)  :: ll3
973           INTEGER, INTENT(IN)  :: ll2, lu2
974           INTEGER, INTENT(IN)  :: ll1, lu1
975           LOGICAL, INTENT(IN) :: ldFlag
976           INTEGER, INTENT(OUT)  :: iErr
977           CHARACTER(len=*), INTENT(IN) :: lMAP
978     
979     ! Local variables
980     !--------------------------------------------------------------------//
981           INTEGER :: lc2, lS2, lE2
982           INTEGER :: lc1, lS1, lE1
983           INTEGER :: I, J, K, IJK
984           LOGICAL :: recheck
985           INTEGER :: IJK_Pg_Owner, proc
986           INTEGER :: gIJK(0:numPEs-1)
987           INTEGER :: I_J_K_Pg(3)
988           INTEGER :: lpCnt
989     
990           CHARACTER(len=32) :: cInt
991     
992     !--------------------------------------------------------------------//
993     
994     ! Initialize Error Flag
995           iErr = 2000
996     
997     ! Initialize the Owner ID
998           IJK_Pg_Owner = UNDEFINED_I
999     
1000     ! Set the initial search region, a single cell.
1001           lS1 = ll1 + (lu1-ll1)/2 + 1; lE1 = lS1
1002           lS2 = ll2 + (lu2-ll2)/2 + 1; lE2 = lS2
1003     
1004           lpCnt = 1
1005           recheck = .TRUE.
1006           do while(recheck)
1007     
1008     ! Initialize the global IJK array to zero. Resetting this array inside
1009     ! this do-loop is most likely overkill. This loop should only cycle
1010     ! if gIJK is zero.
1011              gIJK = 0
1012     
1013     ! Report debugging information for the search region.
1014              if(ldFlag) then
1015                 write(*,"(/5x,'Pass: ',I4)") lpCnt
1016                 write(*,"( 5x,'lp2 bounds:',2(2x,I4))")lS2, lE2
1017                 write(*,"( 5x,'lp1 bounds:',2(2x,I4))")lS1, lE1
1018              endif
1019     
1020              lp2: do lc2 = lS2, lE2
1021              lp1: do lc1 = lS1, lE1
1022     ! Map the loop counters to I/J/K indices.
1023                 SELECT CASE (lMap)
1024                 CASE ('JKI_MAP')
1025                    I=lc1; J=ll3; K=lc2
1026                 CASE ('IKJ_MAP')
1027                    I=ll3; J=lc1; K=lc2
1028                 CASE ('KIJ_MAP')
1029                    I=lc2; J=lc1; K=ll3
1030                 CASE DEFAULT
1031                    iErr = 1001
1032                 END SELECT
1033     
1034     ! Only the rank that owns this I/J/K proceeds.
1035                 if(.NOT.IS_ON_myPE_owns(I,J,K)) cycle
1036     ! Calculate the triple loop index.
1037                 IJK = funijk(I,J,K)
1038     ! If there is fluid at this location, store the IJK and exit loops.
1039                 if(fluid_at(IJK)) then
1040                    gIJK(myPE) = IJK
1041                    exit lp2
1042                 endif
1043              enddo lp1
1044              enddo lp2
1045     
1046     ! Sync gIJK across all processes. Select the lowest ranked process that
1047     ! has gIJK defined. This choice is arbitray and doesn't really matter.
1048     ! It just needs to be consistent.
1049              CALL global_all_sum(gIJK)
1050              proc_lp: do proc=0, numPEs-1
1051                 if(gIJK(proc) /= 0) then
1052                    IJK_P_g = gIJK(proc)
1053                    IJK_Pg_Owner = proc
1054                    recheck = .FALSE.
1055                    exit proc_lp
1056                 endif
1057              enddo proc_lp
1058     
1059     ! If the proceeding section did not find a fluid cell, expand the search
1060     ! area and try again.
1061              if(recheck) then
1062                 if(lS1 > ll1 .OR. lE1 < lu1 .OR.                           &
1063                    lS2 > ll2 .OR. lE2 < lu2) then
1064     ! Expand the 1-axis
1065                    lS1 = max((lS1-1), ll1)
1066                    lE1 = min((lE1+1), lu1)
1067     ! Expand the 2-axis
1068                    lS2 = max((lS2-1), ll2)
1069                    lE2 = min((lE2+1), lu2)
1070     ! The entire seach plane was checked with no fluid cell identified.
1071     ! Force IJK_P_g to undefined for later error checking.
1072                 else
1073                    recheck = .FALSE.
1074                    IJK_P_g = UNDEFINED_I
1075                 endif
1076              endif
1077           enddo
1078     
1079     ! Verify that one fluid cell was detected. Otherwise flag the possible
1080     ! errors and return.
1081           if(IJK_P_G == UNDEFINED_I) then
1082              iErr = 2001
1083              return
1084           elseif(IJK_Pg_Owner == UNDEFINED_I) then
1085              iErr = 2002
1086              return
1087           endif
1088     
1089     
1090     ! The owner if the IJK_Pg gets the global I/J/K values and sends
1091     ! them to all ranks.
1092           if(myPE == IJK_Pg_Owner) then
1093              I_J_K_Pg(1) = I_OF(IJK_P_G)
1094              I_J_K_Pg(2) = J_OF(IJK_P_G)
1095              I_J_K_Pg(3) = K_OF(IJK_P_G)
1096           endif
1097           CALL BCAST(I_J_K_Pg, IJK_Pg_Owner)
1098     
1099           I = I_J_K_Pg(1)
1100           J = I_J_K_Pg(2)
1101           K = I_J_K_Pg(3)
1102     
1103     ! If debugging, have PE_IO report some information before the
1104     ! data is overwritten.
1105           if(ldFlag) then
1106              write(*,"(/3x,'IJK_P_g successfully identified!')")
1107              cInt=''; write(cInt,*) IJK_Pg_Owner
1108              write(*,"( 5x,'Owner Rank: ',A)")trim(adjustl(cInt))
1109              cInt=''; write(cInt,*) IJK_P_G
1110              write(*,"(5x, 'IJK: ',A)", advance='no') trim(adjustl(cInt))
1111              write(*,"(' :: ')", advance='no')
1112              cInt=''; write(cInt,*) I
1113              write(*,"('(',A)",advance='no') trim(adjustl(cInt))
1114              cInt=''; write(cInt,*) J
1115              write(*,"(',',A)",advance='no') trim(adjustl(cInt))
1116              cInt=''; write(cInt,*) K
1117              write(*,"(',',A,')',2/)") trim(adjustl(cInt))
1118           endif
1119     
1120     ! Ranks that 'see' IJK_P_g store their local IJK value. Everyone else
1121     ! resets IJK_P_g to UNDEFINED_I. This removes the need for getting
1122     ! I/J/K values later on in source_PPg.
1123     !      IJK_P_g = merge(funijk(I,J,K), UNDEFINED_I,                      &
1124     !         IS_ON_myPE_plus2layers(I,J,K))
1125     
1126           IF(IS_ON_myPE_plus2layers(I,J,K)) THEN
1127              IJK_P_g = funijk(I,J,K)
1128           ELSE
1129              IJK_P_g = UNDEFINED_I
1130           ENDIF
1131     
1132           IERR = 0
1133           RETURN
1134           END SUBROUTINE IJK_Pg_SEARCH
1135