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