File: /nfs/home/0/users/jenkins/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
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(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(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: undefined, 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
380           use toleranc, only: tmin
381     
382           use indices, only: im1, jm1, km1
383           use functions, only: is_on_mype_plus2layers
384           use boundfunijk, only: bound_funijk
385           use compar, only: dead_cell_at
386           IMPLICIT NONE
387     
388     ! Dummy arguments
389     !--------------------------------------------------------------------//
390     ! index for boundary condition
391           INTEGER, INTENT(IN) :: BCV
392     
393     ! Local variables
394     !--------------------------------------------------------------------//
395     ! indices
396           INTEGER :: I, J, K, IJK, M
397     ! ijk index for setting normal component of velocity
398           INTEGER :: FIJK
399     ! number densities for use in GHD theory only
400           DOUBLE PRECISION :: nM, nTOT
401     ! calculation for normal component of mixture velocity
402           DOUBLE PRECISION :: lvel_s
403     !--------------------------------------------------------------------//
404     
405           DO K = BC_K_B(BCV), BC_K_T(BCV)
406           DO J = BC_J_S(BCV), BC_J_N(BCV)
407           DO I = BC_I_W(BCV), BC_I_E(BCV)
408              IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
409              IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
410              IJK = BOUND_FUNIJK(I,J,K)
411     
412              P_STAR(IJK) = ZERO
413              P_G(IJK) = SCALE(BC_P_G(BCV))
414              EP_G(IJK) = BC_EP_G(BCV)
415              T_G(IJK) = BC_T_G(BCV)
416     
417              IF (NMAX(0) > 0) &
418                X_G(IJK,:NMAX(0)) = BC_X_G(BCV,:NMAX(0))
419     
420              IF (NScalar > 0) &
421                Scalar(IJK,:NScalar) = BC_Scalar(BCV,:NScalar)
422     
423              IF (K_Epsilon) THEN
424                K_Turb_G(IJK) = BC_K_Turb_G(BCV)
425                E_Turb_G(IJK) = BC_E_Turb_G(BCV)
426              ENDIF
427     
428              DO M = 1, SMAX
429                ROP_S(IJK,M) = BC_ROP_S(BCV,M)
430                T_S(IJK,M) = BC_T_S(BCV,M)
431                THETA_M(IJK,M) = BC_THETA_M(BCV,M)
432                IF (NMAX(M) > 0) &
433                   X_S(IJK,M,:NMAX(M)) = BC_X_S(BCV,M,:NMAX(M))
434              ENDDO
435     
436              U_G(IJK) = BC_U_G(BCV)
437              V_G(IJK) = BC_V_G(BCV)
438              W_G(IJK) = BC_W_G(BCV)
439              IF (SMAX > 0) THEN
440                 U_S(IJK,:SMAX) = BC_U_S(BCV,:SMAX)
441                 V_S(IJK,:SMAX) = BC_V_S(BCV,:SMAX)
442                 W_S(IJK,:SMAX) = BC_W_S(BCV,:SMAX)
443              ENDIF
444     
445     ! When the boundary plane is located on the E, N, T side of the domain
446     ! (fluid cell is located w, s, b), set the component of velocity normal
447     ! to the boundary plane of the adjacent fluid cell
448              SELECT CASE (TRIM(BC_PLANE(BCV)))
449                 CASE ('W')
450                    FIJK = BOUND_FUNIJK(IM1(I),J,K)
451                    U_G(FIJK) = BC_U_G(BCV)
452                    IF (SMAX >0) U_S(FIJK,:SMAX) = BC_U_S(BCV,:SMAX)
453                 CASE ('S')
454                    FIJK = BOUND_FUNIJK(I,JM1(J),K)
455                    V_G(FIJK) = BC_V_G(BCV)
456                    IF(SMAX>0) V_S(FIJK,:SMAX) = BC_V_S(BCV,:SMAX)
457                 CASE ('B')
458                    FIJK = BOUND_FUNIJK(I,J,KM1(K))
459                    W_G(FIJK) = BC_W_G(BCV)
460                    IF (SMAX>0) W_S(FIJK,:SMAX) = BC_W_S(BCV,:SMAX)
461              END SELECT
462     
463     ! for GHD theory to compute mixture BC of velocity and density
464              IF(KT_TYPE_ENUM == GHD_2007) THEN
465                 ROP_S(IJK,MMAX) = ZERO
466                 nTOT = ZERO
467                 THETA_M(IJK,MMAX) = ZERO
468                 U_S(IJK,MMAX) =  ZERO
469                 V_S(IJK,MMAX) =  ZERO
470                 W_S(IJK,MMAX) =  ZERO
471                 lvel_s = zero
472                 DO M = 1, SMAX
473                    ROP_S(IJK,MMAX) = ROP_S(IJK,MMAX) + &
474                       BC_ROP_S(BCV,M)
475                    nM = BC_ROP_S(BCV,M)*6d0/ &
476                       (PI*D_p(IJK,M)**3*RO_S(IJK,M))
477                    nTOT = nTOT + nM
478                    THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) + &
479                       nM*BC_THETA_M(BCV,M)
480                    U_S(IJK,MMAX) = U_S(IJK,MMAX) + &
481                       BC_ROP_S(BCV,M)*BC_U_S(BCV,M)
482                    V_S(IJK,MMAX) = V_S(IJK,MMAX) + &
483                       BC_ROP_S(BCV,M)*BC_V_S(BCV,M)
484                    W_S(IJK,MMAX) = W_S(IJK,MMAX) + &
485                       BC_ROP_S(BCV,M)*BC_W_S(BCV,M)
486     ! set velocity component normal to plane in adjacent fluid cell
487                    SELECT CASE (TRIM(BC_PLANE(BCV)))
488                       CASE ('W')
489                          lvel_s = lvel_s + BC_ROP_S(BCV,M)*BC_U_S(BCV,M)
490                       CASE ('S')
491                          lvel_s = lvel_s + BC_ROP_S(BCV,M)*BC_V_S(BCV,M)
492                       CASE ('B')
493                          lvel_s = lvel_s + BC_ROP_S(BCV,M)*BC_W_S(BCV,M)
494                    END SELECT
495                 ENDDO
496                 IF(ROP_S(IJK,MMAX) > ZERO) THEN
497                    THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) / nTOT
498                    U_S(IJK,MMAX) = U_S(IJK,MMAX) / ROP_S(IJK,MMAX)
499                    V_S(IJK,MMAX) = V_S(IJK,MMAX) / ROP_S(IJK,MMAX)
500                    W_S(IJK,MMAX) = W_S(IJK,MMAX) / ROP_S(IJK,MMAX)
501                    lvel_s = lvel_s/ROP_S(IJK,MMAX)
502                 ENDIF
503                 SELECT CASE (TRIM(BC_PLANE(BCV)))
504                    CASE ('W')
505                      FIJK = BOUND_FUNIJK(IM1(I),J,K)
506                      U_S(FIJK,MMAX) = lvel_s
507                    CASE ('S')
508                      FIJK = BOUND_FUNIJK(I,JM1(J),K)
509                      V_S(FIJK,MMAX) = lvel_s
510                    CASE ('B')
511                      FIJK = BOUND_FUNIJK(I,J,KM1(K))
512                      W_S(FIJK,MMAX) = lvel_s
513                 END SELECT
514     
515              ENDIF  ! end if (trim(kt_type) ='ghd')
516     
517     ! Set MMS BCs when MI boundary condition is used.
518              IF (USE_MMS) THEN
519                 P_G(IJK) = SCALE(MMS_P_G(IJK))
520                 EP_G(IJK) = MMS_EP_G(IJK)
521                 T_G(IJK) = MMS_T_G(IJK)
522     
523                 DO M = 1, SMAX
524                    ROP_S(IJK,M) = MMS_ROP_S(IJK)
525                    T_S(IJK,M) = MMS_T_S(IJK)
526                    THETA_M(IJK,M) = MMS_THETA_M(IJK)
527                 ENDDO
528     
529                 U_G(IJK) = MMS_U_G(IJK)
530                 V_G(IJK) = MMS_V_G(IJK)
531                 W_G(IJK) = MMS_W_G(IJK)
532                 IF (SMAX > 0) THEN
533                    U_S(IJK,:SMAX) = MMS_U_S(IJK)
534                    V_S(IJK,:SMAX) = MMS_V_S(IJK)
535                    W_S(IJK,:SMAX) = MMS_W_S(IJK)
536                 ENDIF
537     ! When the boundary plane is W, S, or B the normal component of velocity
538     ! needs to be set for both sides of the boundary cell.
539                 SELECT CASE (TRIM(BC_PLANE(BCV)))
540                    CASE ('W')
541                      FIJK = BOUND_FUNIJK(IM1(I),J,K)
542                      U_G(FIJK) = MMS_U_G(FIJK)
543                      IF(SMAX>0) U_S(FIJK,:SMAX) = MMS_U_S(FIJK)
544                    CASE ('S')
545                      FIJK = BOUND_FUNIJK(I,JM1(J),K)
546                      V_G(FIJK) = MMS_V_G(FIJK)
547                      IF(SMAX>0) V_S(FIJK,:SMAX) = MMS_V_S(FIJK)
548                    CASE ('B')
549                      FIJK = BOUND_FUNIJK(I,J,KM1(K))
550                      W_G(FIJK) = MMS_W_G(FIJK)
551                      IF(SMAX>0) W_S(FIJK,:SMAX) = MMS_W_S(FIJK)
552                 END SELECT
553              ENDIF ! end if(USE_MMS)
554     
555           ENDDO   ! do i
556           ENDDO   ! do j
557           ENDDO   ! do k
558     
559           RETURN
560           END SUBROUTINE SET_BC0_INFLOW
561     
562     
563     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
564     !                                                                      C
565     !  Subroutine: set_bc0_init_jet                                        C
566     !  Purpose: initializing time dependent jet conditions                 C
567     !                                                                      C
568     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
569           SUBROUTINE SET_BC0_INIT_JET(BCV)
570     
571     ! Modules
572     !--------------------------------------------------------------------//
573           use bc, only: bc_plane
574           use bc, only: bc_jet_g, bc_jet_g0
575           use bc, only: bc_dt_0, bc_time
576           use bc, only: bc_u_g, bc_v_g, bc_w_g
577           use param1, only: undefined
578           use run, only: time
579           IMPLICIT NONE
580     
581     ! Dummy arguments
582     !--------------------------------------------------------------------//
583     ! index of boundary
584           INTEGER, INTENT(IN) :: BCV
585     !--------------------------------------------------------------------//
586     
587           BC_JET_G(BCV) = UNDEFINED
588           IF (BC_DT_0(BCV) /= UNDEFINED) THEN
589              BC_TIME(BCV) = TIME + BC_DT_0(BCV)
590              BC_JET_G(BCV) = BC_JET_G0(BCV)
591              IF (BC_JET_G(BCV) /= UNDEFINED) THEN
592                 SELECT CASE (TRIM(BC_PLANE(BCV)))
593                 CASE ('W', 'E')
594                    BC_U_G(BCV) = BC_JET_G(BCV)
595                 CASE ('S', 'N')
596                    BC_V_G(BCV) = BC_JET_G(BCV)
597                 CASE ('B', 'T')
598                    BC_W_G(BCV) = BC_JET_G(BCV)
599                 END SELECT
600              ENDIF
601           ELSE
602              BC_TIME(BCV) = UNDEFINED
603           ENDIF
604           RETURN
605           END SUBROUTINE SET_BC0_INIT_JET
606     
607     
608     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
609     !                                                                      C
610     !  Subroutine: set_bc0_init_bcdt_calcs                                 C
611     !  Purpose: initializing time dependent outflow calculations for       C
612     !  modifying outflow conditions (MO type) or simple reporting          C
613     !  outflow conditions (PO or O types)                                  C
614     !                                                                      C
615     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
616           SUBROUTINE SET_BC0_INIT_BCDT_CALCS(BCV)
617     
618     ! Modules
619     !--------------------------------------------------------------------//
620           use bc, only: bc_dt_0, bc_time
621           use bc, only: bc_out_n
622           use bc, only: bc_mout_g, bc_mout_s
623           use bc, only: bc_vout_g, bc_vout_s
624           use physprop, only: smax
625           use param1, only: undefined, zero
626           use run, only: time
627           IMPLICIT NONE
628     ! Dummy arguments
629     !--------------------------------------------------------------------//
630     ! index of boundary
631           INTEGER, INTENT(IN) :: BCV
632     !--------------------------------------------------------------------//
633     
634     ! initializing for time dependent outflow reporting calculation
635           IF (BC_DT_0(BCV) /= UNDEFINED) THEN
636              BC_TIME(BCV) = TIME + BC_DT_0(BCV)
637              BC_OUT_N(BCV) = 0
638              BC_MOUT_G(BCV) = ZERO
639              BC_VOUT_G(BCV) = ZERO
640              IF (SMAX > 0) THEN
641                 BC_MOUT_S(BCV,:SMAX) = ZERO
642                 BC_VOUT_S(BCV,:SMAX) = ZERO
643              ENDIF
644           ELSE
645              BC_TIME(BCV) = UNDEFINED
646           ENDIF
647           RETURN
648           END SUBROUTINE SET_BC0_INIT_BCDT_CALCS
649     
650     
651     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
652     !                                                                      C
653     !  Subroutine: set_bc0_init_scalars                                    C
654     !  Purpose: Provide a mechanism to initialize certain field variables  C
655     !  in the boundary cell to the value of the neighboring fluid cell.    C
656     !  It is not clear to what extent initialization is necessary as the   C
657     !  governing matrix equation for each of the field variables is set    C
658     !  such that the bc cells are not actually solved but are set to the   C
659     !  value of the adjacent fluid cell(see bc_phi). However, there may be C
660     !  some instance in the code where the bc cells are referenced before  C
661     !  the code has called the governing equation.                         C
662     !  This subroutine illustrates how these particular variables should   C
663     !  not need (or use) bc values defined by the user!                    C
664     !                                                                      C
665     !  Comments: this call replaces some initializations that were         C
666     !  'inadvertently' occurring in set_bc1                                C
667     !                                                                      C
668     !  WARNING: this routine only serves to initialize field variables     C
669     !  whose governing solver routine invokes bc_phi! do not insert any    C
670     !  other initializations here as it is likely not appropriate!!        C
671     !                                                                      C
672     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
673           SUBROUTINE set_bc0_init_bcscalars(BCV)
674     
675     ! Modules
676     !--------------------------------------------------------------------//
677           use bc, only: bc_plane
678           use bc, only: bc_k_b, bc_k_t
679           use bc, only: bc_j_s, bc_j_n
680           use bc, only: bc_i_w, bc_i_e
681     
682           use fldvar, only: t_g, t_s, theta_m
683           use fldvar, only: x_g, x_s, scalar
684           use fldvar, only: k_turb_g, e_turb_g
685           use physprop, only: smax, mmax
686           use physprop, only: nmax
687           use run, only: k_epsilon, kt_type_enum, ghd_2007
688           use scalars, only: nscalar
689     
690           use indices, only: im1, ip1, jm1, jp1, km1, kp1
691           use functions, only: is_on_mype_plus2layers
692           use boundfunijk, only: bound_funijk
693           use compar, only: dead_cell_at
694           IMPLICIT NONE
695     
696     ! Dummy arguments
697     !--------------------------------------------------------------------//
698     ! index of boundary
699           INTEGER, INTENT(IN) :: BCV
700     
701     ! Local variables
702     !--------------------------------------------------------------------//
703     ! indices of current boundary cell
704           INTEGER :: IJK, I, J, K
705     ! index of neighboring fluid cell
706           INTEGER :: FIJK
707     ! local indices
708           INTEGER :: N, M
709     !--------------------------------------------------------------------//
710     
711           DO K = BC_K_B(BCV), BC_K_T(BCV)
712           DO J = BC_J_S(BCV), BC_J_N(BCV)
713           DO I = BC_I_W(BCV), BC_I_E(BCV)
714              IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
715              IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
716              IJK = BOUND_FUNIJK(I,J,K)
717     
718              SELECT CASE (TRIM(BC_PLANE(BCV)))
719                 CASE ('W')   ! fluid cell at west
720                    FIJK = BOUND_FUNIJK(IM1(I),J,K)
721                 CASE ('E')
722                    FIJK = BOUND_FUNIJK(IP1(I),J,K)
723                 CASE ('S')   ! fluid cell at south
724                    FIJK = BOUND_FUNIJK(I,JM1(J),K)
725                 CASE ('N')
726                    FIJK = BOUND_FUNIJK(I,JP1(J),K)
727                 CASE ('B')   ! fluid cell at bottom
728                    FIJK = BOUND_FUNIJK(I,J,KM1(K))
729                 CASE ('T')
730                    FIJK = BOUND_FUNIJK(I,J,KP1(K))
731              END SELECT
732     
733              T_G(IJK) = T_G(FIJK)
734              IF (NMAX(0) > 0) &
735                 X_G(IJK,:NMAX(0)) = X_G(FIJK,:NMAX(0))
736     
737     ! setting scalar quantities
738              IF (NScalar >0) &
739                 Scalar(IJK, :NScalar) = Scalar(FIJK, :NScalar)
740     
741     ! setting turbulence quantities
742              IF(K_Epsilon) THEN
743                 K_Turb_G(IJK) = K_Turb_G(FIJK)
744                 E_Turb_G(IJK) = E_Turb_G(FIJK)
745              ENDIF
746     
747     ! this should only loop of tfm phases
748              DO M = 1, SMAX
749                 T_S(IJK,M) = T_S(FIJK,M)
750                 THETA_M(IJK,M) =  THETA_M(FIJK,M)
751                 IF (NMAX(M) > 0) &
752                    X_S(IJK,M,:NMAX(M)) = X_S(FIJK,M,:NMAX(M))
753              ENDDO
754     
755              IF(KT_TYPE_ENUM == GHD_2007) &
756                  THETA_M(IJK,MMAX) =  THETA_M(FIJK,MMAX)
757     
758           ENDDO
759           ENDDO
760           ENDDO
761           RETURN
762           END SUBROUTINE set_bc0_init_bcscalars
763     
764     
765     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
766     !                                                                      !
767     !  Subroutine: SET_IJK_P_G                                             !
768     !  Purpose: Pick an appropriate control volume to specify Ppg.         !
769     !                                                                      !
770     !  Author: J. Musser                                  Date: 07-Nov-13  !
771     !  Reviewer:                                          Date:            !
772     !                                                                      !
773     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
774           SUBROUTINE SET_IJK_P_G
775     
776     ! IJK location where Ppg is fixed.
777           use bc, only: IJK_P_G
778     ! Specified constant gas density.
779           use physprop, only: RO_G0
780     
781           use geometry, only: CYCLIC
782           use geometry, only: CYCLIC_X, CYCLIC_X_PD, CYCLIC_X_MF
783           use geometry, only: CYCLIC_Y, CYCLIC_Y_PD, CYCLIC_Y_MF
784           use geometry, only: CYCLIC_Z, CYCLIC_Z_PD, CYCLIC_Z_MF
785           use geometry, only: iMAX1, iMin1
786           use geometry, only: jMAX1, jMin1
787           use geometry, only: kMAX1, kMin1
788     
789           use geometry, only: do_K
790     
791           use funits, only: DMP_LOG
792     
793           use bc, only: BC_I_w, BC_I_e
794           use bc, only: BC_J_s, BC_J_n
795           use bc, only: BC_K_b, BC_K_t
796           use bc, only: BC_DEFINED
797           use bc, only: BC_TYPE
798           use bc, only: BC_PLANE
799     
800     ! MFIX Runtime parameters:
801           use param, only: DIMENSION_BC
802           use param1, only: UNDEFINED
803           use param1, only: UNDEFINED_I
804     
805           use mpi_utility
806     
807           implicit none
808     !--------------------------------------------------------------------//
809           INTEGER :: BCV
810     
811           CHARACTER(len=7) :: Map
812           CHARACTER(len=128) :: lMsg
813     
814           INTEGER :: l3
815           INTEGER :: l2, u2
816           INTEGER :: l1, u1
817     
818           LOGICAL, parameter :: setDBG = .FALSE.
819           LOGICAL :: dFlag
820           INTEGER :: iErr
821     
822           EXTERNAL JKI_MAP, IKJ_MAP, KIJ_MAP
823     !--------------------------------------------------------------------//
824     
825           dFlag = (DMP_LOG .AND. setDBG)
826     
827     ! Initialize.
828           iErr = 0
829           IJK_P_G = UNDEFINED_I
830     
831     ! This is not needed for compressible cases.
832           IF(RO_G0 == UNDEFINED) THEN
833              IF(dFlag) write(*,"(3x,A)")                                   &
834                 'Compressible: IJK_P_g remaining undefined.'
835              return
836           ENDIF
837     
838     ! If there are no cyclic boundaries, look for a pressure outflow.
839           lpBCV: DO BCV = 1, DIMENSION_BC
840              IF (.NOT.BC_DEFINED(BCV)) cycle lpBCV
841              IF (BC_TYPE(BCV) == 'P_OUTFLOW' .OR. &
842                  BC_TYPE(BCV) == 'P_INFLOW') THEN
843                 IF(dFlag) write(*,"(3x,A)")                                &
844                    'Outflow PC defined: IJK_P_g remaining undefined.'
845                 RETURN
846              ENDIF
847           ENDDO lpBCV
848     
849     ! Initialize.
850              l3 = UNDEFINED_I
851              l2 = UNDEFINED_I; u2=l2
852              l1 = UNDEFINED_I; u1=l1
853     
854     ! If there are cyclic boundaries, flag a cell along the positive
855     ! domain extreme in the cyclic direction (e.g., JMAX1).
856           IF(CYCLIC_Y .OR. CYCLIC_Y_PD .OR. CYCLIC_Y_MF) THEN
857     
858              Map = 'JKI_MAP'
859              l3 = JMAX1
860              l2 = KMIN1;  u2 = KMAX1
861              l1 = IMIN1;  u1 = IMAX1
862              lMsg='Cyclic in Y'
863     
864           ELSEIF(CYCLIC_X .OR. CYCLIC_X_PD .OR. CYCLIC_X_MF) THEN
865     
866              Map = 'IKJ_MAP'
867              l3 = IMAX1
868              l2 = KMIN1;  u2 = KMAX1
869              l1 = JMIN1;  u1 = JMAX1
870              lMsg='Cyclic in X'
871     
872           ELSEIF(CYCLIC_Z .OR. CYCLIC_Z_PD .OR. CYCLIC_Z_MF) THEN
873     
874              Map = 'KIJ_MAP'
875              l3 = KMAX1
876              l2 = IMIN1;  u2 = IMAX1
877              l1 = JMIN1;  u1 = JMAX1
878              lMsg='Cyclic in Z'
879     
880           ENDIF
881     
882     ! No cyclic boundaries or pressure outflows. The IJ plane is used in
883     ! this case to maximize search region for 2D problems.
884           IF(l3 == UNDEFINED_I) THEN
885              Map = 'KIJ_MAP'
886              l3 = merge((KMAX1 - KMIN1)/2 + 1, KMIN1, do_K)
887              l2 = IMIN1;  u2 = IMAX1
888              l1 = JMIN1;  u1 = JMAX1
889              lMsg='Center of domain'
890           ENDIF
891     
892     ! Debugging messages.
893           IF(dFlag) THEN
894              write(*,"(3/,3x,'Map: ',A)") Map
895              write(*,"(/5x,'l3:',2x,I4)") l3
896              write(*,"( 5x,'l2:',2(2x,I4))") l2, u2
897              write(*,"( 5x,'l1:',2(2x,I4))") l1, u1
898              write(*,"( 5x,'Msg: ',A)") trim(lMsg)
899           ENDIF
900     
901     ! Invoke the search routine.
902           SELECT CASE (Map)
903           CASE ('JKI_MAP')
904              CALL IJK_Pg_SEARCH(l3, l2, u2, l1, u1, JKI_MAP, dFlag, iErr)
905           CASE ('IKJ_MAP')
906              CALL IJK_Pg_SEARCH(l3, l2, u2, l1, u1, IKJ_MAP, dFlag, iErr)
907           CASE ('KIJ_MAP')
908              CALL IJK_Pg_SEARCH(l3, l2, u2, l1, u1, KIJ_MAP, dFlag, iErr)
909           CASE DEFAULT
910              iErr = 1001
911           END SELECT
912     
913           IF(iErr == 0) RETURN
914     
915     ! Error management.
916           IF(DMP_LOG) THEN
917              SELECT CASE (iErr)
918              CASE ( 1001);  WRITE(UNIT_LOG, 1001); WRITE(*,1001)
919              CASE ( 2000);  WRITE(UNIT_LOG, 2000); WRITE(*,2000)
920              CASE ( 2001);  WRITE(UNIT_LOG, 2001); WRITE(*,2001)
921              CASE ( 2002);  WRITE(UNIT_LOG, 2002); WRITE(*,2002)
922              CASE DEFAULT
923                 WRITE(UNIT_LOG, 1000) iErr
924                 WRITE(*,1000) iErr
925              END SELECT
926     
927              WRITE(UNIT_LOG, 9000) MAP(1:1), l3, MAP(2:2),                 &
928                 l2, u2, MAP(3:3), l1, u1
929              WRITE(*, 9000) MAP(1:1), l3, MAP(2:2),                        &
930                 l2, u2, MAP(3:3), l1, u1
931     
932              WRITE(*, 9999)
933              WRITE(UNIT_LOG, 9999)
934     
935           ENDIF
936     
937     
938           CALL MFIX_EXIT(myPE)
939     
940     
941      1000 FORMAT(//1X,70('*')/' From: SET_IJK_Pg',/,                       &
942              ' Error 1000: Unknown error reported. x', I4.4)
943     
944      1001 FORMAT(//1X,70('*')/' From: SET_IJK_Pg',/,                       &
945              ' Error 1001: Invalid mapping function.')
946     
947      2000 FORMAT(//1X,70('*')/' From: SET_IJK_Pg > IJK_Pg_SEARCH',/,       &
948              ' Error 2000: Unknown error reported from IJK_Pg_SEARCH.')
949     
950      2001 FORMAT(//1X,70('*')/' From: SET_IJK_Pg > IJK_Pg_SEARCH',/,       &
951              ' Error 2001: Unable to locate fluid cell in search region.')
952     
953      2002 FORMAT(//1X,70('*')/' From: SET_IJK_Pg > IJK_Pg_SEARCH',/,       &
954              ' Error 2002: Unable to locate fluid cell owner.')
955     
956      9000 FORMAT(/' Search plane information:',/,3x,A1,': ',I8,            &
957               2(/3x,A1,': ',I8,' x ',I8))
958     
959      9999 FORMAT(/' Fatal Error --> Invoking MFIX_EXIT',/1x,70('*'),2/)
960     
961     
962           END SUBROUTINE SET_IJK_P_G
963     
964     
965     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
966     !                                                                      !
967     !  Author: J. Musser                                  Date: 07-Nov-13  !
968     !  Reviewer:                                          Date:            !
969     !                                                                      !
970     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
971           SUBROUTINE IJK_Pg_SEARCH(ll3, ll2, lu2, ll1, lu1, lMAP,          &
972              ldFlag, iErr)
973     
974     ! Modules
975     !--------------------------------------------------------------------//
976     ! IJK location where Ppg is fixed.
977           use bc, only: IJK_P_g
978           use indices
979           use mpi_utility
980           use functions
981           implicit none
982     
983     ! Dummy arguments
984     !--------------------------------------------------------------------//
985           INTEGER, INTENT(IN)  :: ll3
986           INTEGER, INTENT(IN)  :: ll2, lu2
987           INTEGER, INTENT(IN)  :: ll1, lu1
988           LOGICAL, intent(in) :: ldFlag
989           INTEGER, INTENT(OUT)  :: iErr
990     
991     ! Local variables
992     !--------------------------------------------------------------------//
993           EXTERNAL lMAP
994     
995           INTEGER :: lc2, lS2, lE2
996           INTEGER :: lc1, lS1, lE1
997           INTEGER :: I, J, K, IJK
998           LOGICAL :: recheck
999           INTEGER :: IJK_Pg_Owner, proc
1000           INTEGER :: gIJK(0:numPEs-1)
1001           INTEGER :: I_J_K_Pg(3)
1002           INTEGER :: lpCnt
1003     
1004           CHARACTER(len=32) :: cInt
1005     
1006     !--------------------------------------------------------------------//
1007     
1008     ! Initialize Error Flag
1009           iErr = 2000
1010     
1011     ! Initialize the Owner ID
1012           IJK_Pg_Owner = UNDEFINED_I
1013     
1014     ! Set the initial search region, a single cell.
1015           lS1 = ll1 + (lu1-ll1)/2 + 1; lE1 = lS1
1016           lS2 = ll2 + (lu2-ll2)/2 + 1; lE2 = lS2
1017     
1018           lpCnt = 1
1019           recheck = .TRUE.
1020           do while(recheck)
1021     
1022     ! Initialize the global IJK array to zero. Resetting this array inside
1023     ! this do-loop is most likely overkill. This loop should only cycle
1024     ! if gIJK is zero.
1025              gIJK = 0
1026     
1027     ! Report debugging information for the search region.
1028              if(ldFlag) then
1029                 write(*,"(/5x,'Pass: ',I4)") lpCnt
1030                 write(*,"( 5x,'lp2 bounds:',2(2x,I4))")lS2, lE2
1031                 write(*,"( 5x,'lp1 bounds:',2(2x,I4))")lS1, lE1
1032              endif
1033     
1034              lp2: do lc2 = lS2, lE2
1035              lp1: do lc1 = lS1, lE1
1036     ! Map the loop counters to I/J/K indices.
1037                 CALL lMAP(lc1, lc2, ll3, I, J, K)
1038     
1039     ! Only the rank that owns this I/J/K proceeds.
1040                 if(.NOT.IS_ON_myPE_owns(I,J,K)) cycle
1041     ! Calculate the triple loop index.
1042                 IJK = funijk(I,J,K)
1043     ! If there is fluid at this location, store the IJK and exit loops.
1044                 if(fluid_at(IJK)) then
1045                    gIJK(myPE) = IJK
1046                    exit lp2
1047                 endif
1048              enddo lp1
1049              enddo lp2
1050     
1051     ! Sync gIJK across all processes. Select the lowest ranked process that
1052     ! has gIJK defined. This choice is arbitray and doesn't really matter.
1053     ! It just needs to be consistent.
1054              CALL global_all_sum(gIJK)
1055              proc_lp: do proc=0, numPEs-1
1056                 if(gIJK(proc) /= 0) then
1057                    IJK_P_g = gIJK(proc)
1058                    IJK_Pg_Owner = proc
1059                    recheck = .FALSE.
1060                    exit proc_lp
1061                 endif
1062              enddo proc_lp
1063     
1064     ! If the proceeding section did not find a fluid cell, expand the search
1065     ! area and try again.
1066              if(recheck) then
1067                 if(lS1 > ll1 .OR. lE1 < lu1 .OR.                           &
1068                    lS2 > ll2 .OR. lE2 < lu2) then
1069     ! Expand the 1-axis
1070                    lS1 = max((lS1-1), ll1)
1071                    lE1 = min((lE1+1), lu1)
1072     ! Expand the 2-axis
1073                    lS2 = max((lS2-1), ll2)
1074                    lE2 = min((lE2+1), lu2)
1075     ! The entire seach plane was checked with no fluid cell identified.
1076     ! Force IJK_P_g to undefined for later error checking.
1077                 else
1078                    recheck = .FALSE.
1079                    IJK_P_g = UNDEFINED_I
1080                 endif
1081              endif
1082           enddo
1083     
1084     ! Verify that one fluid cell was detected. Otherwise flag the possible
1085     ! errors and return.
1086           if(IJK_P_G == UNDEFINED_I) then
1087              iErr = 2001
1088              return
1089           elseif(IJK_Pg_Owner == UNDEFINED_I) then
1090              iErr = 2002
1091              return
1092           endif
1093     
1094     
1095     ! The owner if the IJK_Pg gets the global I/J/K values and sends
1096     ! them to all ranks.
1097           if(myPE == IJK_Pg_Owner) then
1098              I_J_K_Pg(1) = I_OF(IJK_P_G)
1099              I_J_K_Pg(2) = J_OF(IJK_P_G)
1100              I_J_K_Pg(3) = K_OF(IJK_P_G)
1101           endif
1102           CALL BCAST(I_J_K_Pg, IJK_Pg_Owner)
1103     
1104           I = I_J_K_Pg(1)
1105           J = I_J_K_Pg(2)
1106           K = I_J_K_Pg(3)
1107     
1108     ! If debugging, have PE_IO report some information before the
1109     ! data is overwritten.
1110           if(ldFlag) then
1111              write(*,"(/3x,'IJK_P_g successfully identified!')")
1112              cInt=''; write(cInt,*) IJK_Pg_Owner
1113              write(*,"( 5x,'Owner Rank: ',A)")trim(adjustl(cInt))
1114              cInt=''; write(cInt,*) IJK_P_G
1115              write(*,"(5x, 'IJK: ',A)", advance='no') trim(adjustl(cInt))
1116              write(*,"(' :: ')", advance='no')
1117              cInt=''; write(cInt,*) I
1118              write(*,"('(',A)",advance='no') trim(adjustl(cInt))
1119              cInt=''; write(cInt,*) J
1120              write(*,"(',',A)",advance='no') trim(adjustl(cInt))
1121              cInt=''; write(cInt,*) K
1122              write(*,"(',',A,')',2/)") trim(adjustl(cInt))
1123           endif
1124     
1125     ! Ranks that 'see' IJK_P_g store their local IJK value. Everyone else
1126     ! resets IJK_P_g to UNDEFINED_I. This removes the need for getting
1127     ! I/J/K values later on in source_PPg.
1128     !      IJK_P_g = merge(funijk(I,J,K), UNDEFINED_I,                      &
1129     !         IS_ON_myPE_plus2layers(I,J,K))
1130     
1131           IF(IS_ON_myPE_plus2layers(I,J,K)) THEN
1132              IJK_P_g = funijk(I,J,K)
1133           ELSE
1134              IJK_P_g = UNDEFINED_I
1135           ENDIF
1136     
1137           IERR = 0
1138           RETURN
1139           END SUBROUTINE IJK_Pg_SEARCH
1140     
1141     
1142     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
1143     !                                                                      !
1144     !  Author: J. Musser                                  Date: 09-Oct-13  !
1145     !  Reviewer:                                          Date:            !
1146     !                                                                      !
1147     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
1148           SUBROUTINE JKI_MAP(in1, in2, in3, lI, lJ, lK)
1149           implicit none
1150     
1151     ! Dummy arguments
1152     !--------------------------------------------------------------------//
1153           INTEGER, intent(in) :: in1, in2, in3
1154           INTEGER, intent(out) :: lI, lJ, lK
1155     
1156           lI=in1; lJ=in3; lK=in2; return
1157           return
1158           END SUBROUTINE JKI_MAP
1159     
1160     
1161           SUBROUTINE IKJ_MAP(in1, in2, in3, lI, lJ, lK)
1162           implicit none
1163     ! Dummy arguments
1164     !--------------------------------------------------------------------//
1165           INTEGER, intent(in) :: in1, in2, in3
1166           INTEGER, intent(out) :: lI, lJ, lK
1167     
1168           lI=in3; lJ=in1; lK=in2; return
1169           return
1170           END SUBROUTINE IKJ_MAP
1171     
1172     
1173           SUBROUTINE KIJ_MAP(in1, in2, in3, lI, lJ, lK)
1174           implicit none
1175     ! Dummy arguments
1176     !--------------------------------------------------------------------//
1177           INTEGER, intent(in) :: in1, in2, in3
1178           INTEGER, intent(out) :: lI, lJ, lK
1179     
1180           lI=in2; lJ=in1; lK=in3; return
1181           return
1182           END SUBROUTINE KIJ_MAP
1183