File: N:\mfix\model\calc_trd_s.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_trD_s                                              C
4     !  Purpose: Calculate the trace of the solids phase rate of strain     C
5     !  tensor at i, j, k                                                   C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 19-DEC-96  C
8     !                                                                      C
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
10           SUBROUTINE CALC_TRD_S(lTRD_S)
11     
12     ! Modules
13     !-----------------------------------------------
14           USE param
15           USE param1
16           USE parallel
17           USE geometry
18           USE fldvar
19           USE indices
20           USE physprop
21           USE compar
22           USE sendrecv
23           USE functions
24           USE cutcell
25           IMPLICIT NONE
26     
27     ! Dummy arguments
28     !-----------------------------------------------
29     ! Strain rate tensor components for mth solids phase
30           DOUBLE PRECISION, INTENT(OUT) :: ltrD_s(DIMENSION_3, DIMENSION_M)
31     
32     ! Local variables
33     !-----------------------------------------------
34     ! Indices
35           INTEGER :: I, J, K
36           INTEGER :: IJK, IMJK, IJMK, IJKM
37           INTEGER :: IM, M
38     
39     !-----------------------------------------------
40     
41           DO M = 1, MMAX
42     !!$omp    parallel do private(ijk,i,j,k,im,imjk,ijmk,ijkm)
43              DO IJK = ijkstart3, ijkend3
44                 IF (.NOT.WALL_AT(IJK)) THEN
45                    I = I_OF(IJK)
46                    J = J_OF(IJK)
47                    K = K_OF(IJK)
48                    IM = IM1(I)
49                    IMJK = IM_OF(IJK)
50                    IJMK = JM_OF(IJK)
51                    IJKM = KM_OF(IJK)
52     
53                    IF(.NOT.CUT_CELL_AT(IJK)) THEN
54     ! at i, j, k:
55     ! du/dx + dv/dy + 1/x dw/dz + u/x =
56     ! 1/x d(xu)/dx + dv/dy + 1/x dw/dz =
57                       lTRD_S(IJK,M) = (X_E(I)*U_S(IJK,M)-&
58                          X_E(IM)*U_S(IMJK,M))*OX(I)*ODX(I) + &
59                          (V_S(IJK,M)-V_S(IJMK,M))*ODY(J) + &
60                          (W_S(IJK,M)-W_S(IJKM,M))*(OX(I)*ODZ(K))
61                    ELSE
62                       CALL CALC_CG_TRD_S(IJK,M,ltrD_s)
63                    ENDIF
64     
65                 ELSE
66                    ltrD_s(IJK,M) = ZERO
67                 ENDIF
68              ENDDO
69           ENDDO
70     
71           END SUBROUTINE CALC_TRD_S
72     
73     
74     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
75     !                                                                      C
76     !  Subroutine: CALC_CG_trD_s                                           C
77     !  Purpose: Calculate the trace of the solids phase rate of strain     C
78     !  tensor at i, j, k with cartesian grid modifications                 C
79     !                                                                      C
80     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
81     !                                                                      C
82     !                                                                      C
83     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
84           SUBROUTINE CALC_CG_TRD_S(IJK,M,ltrD_s)
85     
86     ! Modules
87     !-----------------------------------------------
88           USE param
89           USE param1
90           USE parallel
91           USE geometry
92           USE fldvar
93           USE indices
94           USE physprop
95           USE compar
96           USE sendrecv
97           USE functions
98           USE cutcell
99     
100           USE bc
101           USE cutcell
102           USE quadric
103           IMPLICIT NONE
104     
105     ! Dummy arguments
106     !-----------------------------------------------
107     ! index
108           INTEGER, INTENT(IN) :: IJK
109     ! phase index
110           INTEGER, INTENT(IN) :: M
111     ! Strain rate tensor components for mth solids phase
112           DOUBLE PRECISION, INTENT(INOUT) :: ltrD_s(DIMENSION_3, DIMENSION_M)
113     
114     ! Local variables
115     !-----------------------------------------------
116     ! Indices
117           INTEGER :: I, J, K
118           INTEGER :: IMJK, IJMK, IJKM
119           INTEGER :: IM
120     
121           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
122           DOUBLE PRECISION :: dudx,dvdy,dwdz
123           DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Vi,Wi,Sx,Sy,Sz
124           DOUBLE PRECISION :: UW_s,VW_s,WW_s
125     
126           LOGICAL :: U_NODE_AT_E, U_NODE_AT_W
127           LOGICAL :: V_NODE_AT_N, V_NODE_AT_S
128           LOGICAL :: W_NODE_AT_T, W_NODE_AT_B
129           INTEGER :: BCV
130           INTEGER :: BCT
131     !-----------------------------------------------
132     
133           I = I_OF(IJK)
134           J = J_OF(IJK)
135           K = K_OF(IJK)
136           IM = IM1(I)
137           IMJK = IM_OF(IJK)
138           IJMK = JM_OF(IJK)
139           IJKM = KM_OF(IJK)
140     
141           IF(FLOW_AT(IJK)) THEN
142               lTRD_S(IJK,M) = ZERO
143              RETURN
144           ENDIF
145     
146           BCV = BC_ID(IJK)
147     
148           IF(BCV > 0 ) THEN
149              BCT = BC_TYPE_ENUM(BCV)
150           ELSE
151              BCT = NONE
152           ENDIF
153           SELECT CASE (BCT)
154              CASE (CG_NSW)
155                 NOC_TRDS = .TRUE.
156                 UW_s = ZERO
157                 VW_s = ZERO
158                 WW_s = ZERO
159              CASE (CG_FSW)
160                 NOC_TRDS = .FALSE.
161                 UW_s = ZERO
162                 VW_s = ZERO
163                 WW_s = ZERO
164              CASE(CG_PSW)
165                 IF(BC_JJ_PS(BCV)==1) THEN   ! Johnson-Jackson partial slip bc
166                    NOC_TRDS = .FALSE.
167                    UW_s = BC_UW_S(BCV,M)
168                    VW_s = BC_VW_S(BCV,M)
169                    WW_s = BC_WW_S(BCV,M)
170                 ELSEIF(BC_HW_S(BCV,M)==UNDEFINED) THEN   ! same as NSW
171                    NOC_TRDS = .TRUE.
172                    UW_s = BC_UW_S(BCV,M)
173                    VW_s = BC_VW_S(BCV,M)
174                    WW_s = BC_WW_S(BCV,M)
175                 ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN    ! same as FSW
176                    NOC_TRDS = .FALSE.
177                 ELSE                                 ! partial slip
178                    NOC_TRDS = .FALSE.
179                    UW_s = ZERO
180                    VW_s = ZERO
181                    WW_s = ZERO
182                 ENDIF
183              CASE (CG_MI)
184                 lTRD_S(IJK,M) = ZERO
185                 RETURN
186              CASE (CG_PO)
187                 lTRD_S(IJK,M) = ZERO
188                 RETURN
189              CASE (NONE)
190                 lTRD_S(IJK,M) = ZERO
191                 RETURN
192           END SELECT
193     
194     
195     ! du/dx
196     !=======================================================================
197           U_NODE_AT_E = ((.NOT.BLOCKED_U_CELL_AT(IJK)) .AND.&
198                          (.NOT.WALL_U_AT(IJK)))
199           U_NODE_AT_W = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.&
200                          (.NOT.WALL_U_AT(IMJK)))
201           IF(U_NODE_AT_E.AND.U_NODE_AT_W) THEN
202              Ui = HALF * (U_S(IJK,M) + U_S(IMJK,M))
203              Xi = HALF * (X_U(IJK) + X_U(IMJK))
204              Yi = HALF * (Y_U(IJK) + Y_U(IMJK))
205              Zi = HALF * (Z_U(IJK) + Z_U(IMJK))
206              Sx = X_U(IJK) - X_U(IMJK)
207              Sy = Y_U(IJK) - Y_U(IMJK)
208              Sz = Z_U(IJK) - Z_U(IMJK)
209              CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
210              IF(Sx /= ZERO) THEN
211                 dudx =  (U_S(IJK,M) - U_S(IMJK,M))/Sx
212                 IF(NOC_TRDS) dudx = dudx - ((Ui-UW_s)/(Sx*DEL_H)*&
213                    (Sy*Ny+Sz*Nz))
214              ELSE
215                 dudx = ZERO
216              ENDIF
217     
218           ELSE
219              dudx = ZERO
220           ENDIF
221     
222     ! dv/dy
223     !=======================================================================
224           V_NODE_AT_N = ((.NOT.BLOCKED_V_CELL_AT(IJK)) .AND.&
225                          (.NOT.WALL_V_AT(IJK)))
226           V_NODE_AT_S = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.&
227                          (.NOT.WALL_V_AT(IJMK)))
228           IF(V_NODE_AT_N.AND.V_NODE_AT_S) THEN
229              Vi = HALF * (V_S(IJK,M) + V_S(IJMK,M))
230              Xi = HALF * (X_V(IJK) + X_V(IJMK))
231              Yi = HALF * (Y_V(IJK) + Y_V(IJMK))
232              Zi = HALF * (Z_V(IJK) + Z_V(IJMK))
233              Sx = X_V(IJK) - X_V(IJMK)
234              Sy = Y_V(IJK) - Y_V(IJMK)
235              Sz = Z_V(IJK) - Z_V(IJMK)
236              CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
237              IF(Sy /= ZERO) THEN
238                 dvdy =  (V_S(IJK,M) - V_S(IJMK,M))/Sy
239                 IF(NOC_TRDS) dvdy = dvdy - ((Vi-VW_s)/(Sy*DEL_H)*&
240                    (Sx*Nx+Sz*Nz))
241              ELSE
242                 dvdy =  ZERO
243              ENDIF
244           ELSE
245              dvdy = ZERO
246           ENDIF
247     
248     ! dw/dz
249     !=======================================================================
250           IF(NO_K) THEN
251              dwdz = ZERO
252           ELSE
253              W_NODE_AT_T = ((.NOT.BLOCKED_W_CELL_AT(IJK)) .AND.&
254                             (.NOT.WALL_W_AT(IJK)))
255              W_NODE_AT_B = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.&
256                             (.NOT.WALL_W_AT(IJKM)))
257              IF(W_NODE_AT_T.AND.W_NODE_AT_B) THEN
258                 Wi = HALF * (W_S(IJK,M) + W_S(IJKM,M))
259                 Xi = HALF * (X_W(IJK) + X_W(IJKM))
260                 Yi = HALF * (Y_W(IJK) + Y_W(IJKM))
261                 Zi = HALF * (Z_W(IJK) + Z_W(IJKM))
262                 Sx = X_W(IJK) - X_W(IJKM)
263                 Sy = Y_W(IJK) - Y_W(IJKM)
264                 Sz = Z_W(IJK) - Z_W(IJKM)
265                 CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
266                 IF(Sz /= ZERO) THEN
267                    dwdz =  (W_S(IJK,M) - W_S(IJKM,M))/Sz
268                    IF(NOC_TRDS) dwdz = dwdz - ((Wi-WW_s)/(Sz*DEL_H)*&
269                       (Sx*Nx+Sy*Ny))
270                 ELSE
271                    dwdz = ZERO
272                 ENDIF
273              ELSE
274                 dwdz = ZERO
275              ENDIF
276           ENDIF  ! NO_K
277     
278           lTRD_S(IJK,M) = dudx + dvdy + dwdz
279     
280           RETURN
281           END SUBROUTINE CALC_CG_TRD_S
282     
283     
284     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
285     !                                                                      C
286     !  Subroutine: GET_FACE_VEL_SOLIDS                                     C
287     !  Purpose: Evaluate the velocity components at each of the faces of   C
288     !  a scalar cell                                                       C
289     !                                                                      C
290     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
291           SUBROUTINE GET_FACE_VEL_SOLIDS(IJK, M, us_e, usw, usn, uss, ust, usb, uscc, &
292                                            vse, vsw, vsn, vss, vst, vsb, &
293                                            wse, wsw, wsn, wss, wst, wsb, wscc)
294     
295     ! Modules
296     !-----------------------------------------------------------------------
297           use fldvar, only: u_s, v_s, w_s
298           use functions, only: im_of, jm_of, km_of
299           use functions, only: ip_of, jp_of, kp_of
300           use functions, only: is_at_e, is_at_n, is_at_t
301           use functions, only: wall_at
302           use fun_avg, only: avg_x, avg_x_e
303           use fun_avg, only: avg_y, avg_y_n
304           use fun_avg, only: avg_z, avg_z_t
305           use geometry, only: cylindrical
306           use geometry, only: xlength, odx_e
307           use indices, only: i_of, j_of, k_of
308           use indices, only: im1, jm1, km1
309           USE is, only: any_is_defined
310           use param1, only: zero, one
311           use run, only: shear, V_SH
312           Use vshear, only: VSH
313           IMPLICIT NONE
314     
315     ! Dummy arguments
316     !-----------------------------------------------------------------------
317     ! index
318           INTEGER, INTENT(IN) :: IJK
319     ! solids index
320           INTEGER, INTENT(IN) :: M
321     ! U_s at the east (i+1/2, j, k) and west face (i-1/2, j, k)
322           DOUBLE PRECISION, INTENT(OUT) :: Us_E, UsW
323     ! U_s at the north (i, j+1/2, k) and south face (i, j-1/2, k)
324           DOUBLE PRECISION, INTENT(OUT) :: UsN, UsS
325     ! U_s at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
326           DOUBLE PRECISION, INTENT(OUT) :: UsT, UsB
327     ! U_s at the center of a scalar cell (i, j, k)
328     ! Calculated for Cylindrical coordinates only.
329           DOUBLE PRECISION, INTENT(OUT) :: UscC
330     
331     ! V_s at the east (i+1/2, j, k) and west face (i-1/2, j, k)
332           DOUBLE PRECISION, INTENT(OUT) :: VsE, VsW
333     ! V_s at the north (i, j+1/2, k) and south face (i, j-1/2, k)
334           DOUBLE PRECISION, INTENT(OUT) :: VsN, VsS
335     ! V_s at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
336           DOUBLE PRECISION, INTENT(OUT) :: VsT, VsB
337     
338     ! W_s at the east (i+1/2, j, k) and west face (i-1/2, j, k)
339           DOUBLE PRECISION, INTENT(OUT) :: WsE, WsW
340     ! W_s at the north (i, j+1/2, k) and south face (i, j-1/2, k)
341           DOUBLE PRECISION, INTENT(OUT) :: WsN, WsS
342     ! W_s at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
343           DOUBLE PRECISION, INTENT(OUT) :: WsT, WsB
344     ! W_s at the center of a scalar cell (i, j, k).
345     ! Calculated for Cylindrical coordinates only.
346           DOUBLE PRECISION, INTENT(OUT) :: WscC
347     
348     ! Local variables
349     !-----------------------------------------------------------------------
350     ! Cell indices
351           INTEGER :: I, J, K, IM, JM, KM
352           INTEGER :: IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
353           INTEGER :: IMJPK, IMJMK, IMJKP, IMJKM, IPJKM, IPJMK
354           INTEGER :: IJMKP, IJMKM, IJPKM
355     ! shear value calculation
356           DOUBLE PRECISION :: SRT
357     !-----------------------------------------------------------------------
358     
359           I = I_OF(IJK)
360           J = J_OF(IJK)
361           K = K_OF(IJK)
362           IM = IM1(I)
363           JM = JM1(J)
364           KM = KM1(K)
365           IMJK = IM_OF(IJK)
366           IPJK = IP_OF(IJK)
367           IJMK = JM_OF(IJK)
368           IJPK = JP_OF(IJK)
369           IJKM = KM_OF(IJK)
370           IJKP = KP_OF(IJK)
371           IMJPK = IM_OF(IJPK)
372           IMJMK = IM_OF(IJMK)
373           IMJKP = IM_OF(IJKP)
374           IMJKM = IM_OF(IJKM)
375           IPJKM = IP_OF(IJKM)
376           IPJMK = IP_OF(IJMK)
377           IJMKP = JM_OF(IJKP)
378           IJMKM = JM_OF(IJKM)
379           IJPKM = JP_OF(IJKM)
380     
381           UsN = AVG_Y(AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I),&
382                       AVG_X_E(U_s(IMJPK, M), U_s(IJPK, M), I), J)   !i, j+1/2, k
383           UsS = AVG_Y(AVG_X_E(U_s(IMJMK, M), U_s(IJMK, M), I),&
384                       AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I), JM)    !i, j-1/2, k
385           UsT = AVG_Z(AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I),&
386                       AVG_X_E(U_s(IMJKP, M), U_s(IJKP, M), I), K)   !i, j, k+1/2
387           UsB = AVG_Z(AVG_X_E(U_s(IMJKM, M), U_s(IJKM, M), I),&
388                       AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I), KM)    !i, j, k-1/2
389           Us_E = U_S(IJK,M)
390           UsW = U_S(IMJK,M)
391     
392           IF (SHEAR)  THEN
393              SRT=(2d0*V_sh/XLENGTH)
394              VsE = AVG_X(AVG_Y_N(V_s(IJMK, M), V_s(IJK, M)),&
395                          AVG_Y_N((V_s(IPJMK, M) - VSH(IPJMK) + &
396                                   VSH(IJMK) + SRT*ONE/oDX_E(I)), &
397                                  (V_s(IPJK, M) - VSH(IPJK) + &
398                                   VSH(IJK) + SRT*ONE/oDX_E(I))), I) !i+1/2, j, k
399              VsW = AVG_X(AVG_Y_N((V_s(IMJMK, M) - VSH(IMJMK) + &
400                                   VSH(IJMK) - SRT*ONE/oDX_E(IM)),&
401                                  (V_s(IMJK, M) - VSH(IMJK) + &
402                                   VSH(IJK) - SRT*ONE/oDX_E(IM)) ),&
403                          AVG_Y_N(V_s(IJMK, M), V_s(IJK, M)), IM)    !i-1/2, j, k
404           ELSE
405              VsE = AVG_X(AVG_Y_N(V_s(IJMK, M), V_s(IJK, M)),&
406                          AVG_Y_N(V_s(IPJMK, M), V_s(IPJK, M)), I )  !i+1/2, j, k
407              VsW = AVG_X(AVG_Y_N(V_s(IMJMK, M), V_s(IMJK, M)),&
408                          AVG_Y_N(V_s(IJMK, M), V_s(IJK, M)), IM )   !i-1/2, j, k
409           ENDIF
410           VsT = AVG_Z(AVG_Y_N(V_s(IJMK, M), V_s(IJK, M)),&
411                       AVG_Y_N(V_s(IJMKP, M), V_s(IJKP, M)), K )     !i, j, k+1/2
412           VsB = AVG_Z(AVG_Y_N(V_s(IJMKM, M), V_s(IJKM, M)),&
413                       AVG_Y_N(V_s(IJMK, M), V_s(IJK, M)), KM )      !i, j, k-1/2
414           VsN = V_S(IJK,M)
415           VsS = V_S(IJMK,M)
416     
417           WsN = AVG_Y(AVG_Z_T(W_s(IJKM, M), W_s(IJK, M)),&
418                       AVG_Z_T(W_s(IJPKM, M), W_s(IJPK, M)), J )     !i, j+1/2, k
419           WsS = AVG_Y(AVG_Z_T(W_s(IJMKM, M), W_s(IJMK, M)),&
420                       AVG_Z_T(W_s(IJKM, M), W_s(IJK, M)), JM )      !i, j-1/2, k
421           WsE = AVG_X(AVG_Z_T(W_s(IJKM, M), W_s(IJK, M)),&
422                       AVG_Z_T(W_s(IPJKM, M), W_s(IPJK, M)), I)      !i+1/2, j, k
423           WsW = AVG_X(AVG_Z_T(W_s(IMJKM, M), W_s(IMJK, M)),&
424                       AVG_Z_T(W_s(IJKM, M), W_s(IJK, M)), IM )      !i-1/2, j, k
425           WsT = W_S(IJK,M)
426           WsB = W_S(IJKM,M)
427     
428           IF(CYLINDRICAL) THEN
429              UscC = AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I)          !i, j, k
430              WscC = AVG_Z_T(W_s(IJKM, M), W_s(IJK, M))             !i, j, k
431           ELSE
432              UscC = ZERO
433              WscC = ZERO
434           ENDIF
435     
436     ! Check for IS surfaces and modify solids velocity-comp accordingly
437           IF(ANY_IS_DEFINED) THEN
438              IF(IS_AT_N(IJK)  .AND. .NOT.WALL_AT(IJPK)) &
439                 UsN = AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I)
440              IF(IS_AT_N(IJMK) .AND. .NOT.WALL_AT(IJMK)) &
441                 UsS = AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I)
442              IF(IS_AT_T(IJK)  .AND. .NOT.WALL_AT(IJKP)) &
443                 UsT = AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I)
444              IF(IS_AT_T(IJKM) .AND. .NOT.WALL_AT(IJKM)) &
445                 UsB = AVG_X_E(U_s(IMJK, M), U_s(IJK, M), I)
446     
447              IF(IS_AT_E(IJK)  .AND. .NOT.WALL_AT(IPJK)) &
448                 VsE = AVG_Y_N(V_s(IJMK, M), V_s(IJK, M))
449              IF(IS_AT_E(IMJK) .AND. .NOT.WALL_AT(IMJK)) &
450                 VsW = AVG_Y_N(V_s(IJMK, M), V_s(IJK, M))
451              IF(IS_AT_T(IJK)  .AND. .NOT.WALL_AT(IJKP)) &
452                 VsT = AVG_Y_N(V_s(IJMK, M), V_s(IJK, M))
453              IF(IS_AT_T(IJKM) .AND. .NOT.WALL_AT(IJKM)) &
454                 VsB = AVG_Y_N(V_s(IJMK, M), V_s(IJK, M))
455     
456              IF(IS_AT_N(IJK)  .AND. .NOT.WALL_AT(IJPK)) &
457                 WsN = AVG_Z_T(W_s(IJKM, M), W_s(IJK, M))
458              IF(IS_AT_N(IJMK) .AND. .NOT.WALL_AT(IJMK)) &
459                 WsS = AVG_Z_T(W_s(IJKM, M), W_s(IJK, M))
460              IF(IS_AT_E(IJK)  .AND. .NOT.WALL_AT(IPJK)) &
461                 WsE = AVG_Z_T(W_s(IJKM, M), W_s(IJK, M))
462              IF(IS_AT_E(IMJK) .AND. .NOT.WALL_AT(IMJK)) &
463                 WsW = AVG_Z_T(W_s(IJKM, M), W_s(IJK, M))
464           ENDIF
465     
466           RETURN
467           END SUBROUTINE GET_FACE_VEL_SOLIDS
468     
469     
470     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
471     !                                                                      C
472     !  Subroutine: CALC_DERIV_VEL_SOLIDS                                   C
473     !  Purpose: Calculate the gradient of the solids phase velocity and    C
474     !  the rate of strain tensor at i, j, k                                C
475     !                                                                      C
476     !         |  du/dx    du/dy   du/dz  |                                 C
477     ! DELV =  |  dv/dx    dv/dy   dv/dz  |  =  dVi/dxj                     C
478     !         |  dw/dx    dw/dy   dw/dz  |                                 C
479     !                                                                      C
480     ! 1/2(DELV+DELV^T) = 1/2(dVi/dxj+dVj/dxi)                              C
481     !                                                                      C
482     !                                                                      C
483     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
484           SUBROUTINE CALC_DERIV_VEL_SOLIDS(IJK, M, lVelGradS, lRateStrainS)
485     
486     ! Modules
487     !-----------------------------------------------------------------------
488           use cutcell, only: cut_cell_at
489           use geometry, only: odx, ody, odz
490           use geometry, only: ox
491           use indices, only: i_of, j_of, k_of
492           use param1, only: half
493     ! Dummy arguments
494     !-----------------------------------------------------------------------
495     ! index
496           INTEGER, INTENT(IN) :: IJK
497     ! solids phase index
498           INTEGER, INTENT(IN) :: M
499     ! gradient of velocity
500           DOUBLE PRECISION, INTENT(OUT) :: lVelGradS(3,3)    ! delV
501     ! rate of strain tensor
502           DOUBLE PRECISION, INTENT(OUT) :: lRateStrainS(3,3) ! D_s
503     
504     ! Local variables
505     !-----------------------------------------------------------------------
506     ! Cell indices
507           INTEGER :: I, J, K
508     ! face values of velocity
509           DOUBLE PRECISION :: us_e, usw, usn, uss, ust, usb, uscc
510           DOUBLE PRECISION :: vse, vsw, vsn, vss, vst, vsb
511           DOUBLE PRECISION :: wse, wsw, wsn, wss, wst, wsb, wscc
512     !-----------------------------------------------------------------------
513     
514           IF(.NOT.CUT_CELL_AT(IJK)) THEN
515              I = I_OF(IJK)
516              J = J_OF(IJK)
517              K = K_OF(IJK)
518     
519     ! Set the velocity components at the scalar faces
520              CALL GET_FACE_VEL_SOLIDS(IJK, M, &
521                 us_e, usw, usn, uss, ust, usb, uscc, &
522                  vse, vsw, vsn, vss, vst, vsb, &
523                  wse, wsw, wsn, wss, wst, wsb, wscc)
524     
525     ! Gradient of gas velocity  at cell center (i, j, k)
526              lVelGradS(1,1) = (Us_E-UsW)*ODX(I)
527              lVelGradS(1,2) = (UsN-UsS)*ODY(J)
528              lVelGradS(1,3) = (UsT-UsB)*(OX(I)*ODZ(K))-WscC*OX(I)
529     
530              lVelGradS(2,1) = (VsE-VsW)*ODX(I)
531              lVelGradS(2,2) = (VsN-VsS)*ODY(J)
532              lVelGradS(2,3) = (VsT-VsB)*(OX(I)*ODZ(K))
533     
534              lVelGradS(3,1) = (WsE-WsW)*ODX(I)
535              lVelGradS(3,2) = (WsN-WsS)*ODY(J)
536              lVelGradS(3,3) = (WsT-WsB)*(OX(I)*ODZ(K)) + UscC*OX(I)
537     
538     
539     ! Strain rate tensor, D_S, at cell center (i, j, k)
540     ! or calculated as 0.5*(DelS+DelS^T)
541              lRateStrainS(1,1) = (Us_E-UsW)*ODX(I)
542              lRateStrainS(1,2) = HALF*((UsN-UsS)*ODY(J)+&
543                                        (VsE-VsW)*ODX(I))
544              lRateStrainS(1,3) = HALF*((WsE-WsW)*ODX(I)+&
545                                        (UsT-UsB)*(OX(I)*ODZ(K))-&
546                                  WscC*OX(I))
547     
548              lRateStrainS(2,1) = lRateStrainS(1,2)
549              lRateStrainS(2,2) = (VsN-VsS)*ODY(J)
550              lRateStrainS(2,3) = HALF*((VsT-VsB)*(OX(I)*ODZ(K))+&
551                                        (WsN-WsS)*ODY(J))
552     
553              lRateStrainS(3,1) = lRateStrainS(1,3)
554              lRateStrainS(3,2) = lRateStrainS(2,3)
555              lRateStrainS(3,3) = (WsT-WsB)*(OX(I)*ODZ(K)) +&
556                                  UscC*OX(I)
557     
558           ELSE   ! cut cell
559              CALL CALC_CG_DERIV_VEL_SOLIDS(IJK,M,lVelGradS,lRateStrainS)
560           ENDIF
561     
562           RETURN
563           END SUBROUTINE CALC_DERIV_VEL_SOLIDS
564     
565     
566     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
567     !                                                                      C
568     !  Module name: CALC_CG_DERIV_VEL_SOLIDS                               C
569     !  Purpose: Calculate velocity derivatives in scalar cut-cell          C
570     !           Solids phase                                               C
571     !                                                                      C
572     !  Author: Jeff Dietiker                              Date: 25-JAN-96  C
573     !                                                                      C
574     !                                                                      C
575     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
576           SUBROUTINE CALC_CG_DERIV_VEL_SOLIDS(IJK, M, lVelGradS, lRateStrainS)
577     
578     ! Modules
579     !-----------------------------------------------
580           USE param
581           USE param1
582           USE parallel
583           USE geometry
584           USE fldvar
585           USE indices
586           USE compar
587           USE sendrecv
588           USE bc
589           USE cutcell
590           USE quadric
591           USE functions
592           IMPLICIT NONE
593     
594     ! Dummy arguments
595     !-----------------------------------------------
596     ! index
597           INTEGER, INTENT(IN) :: IJK
598     ! solids phase index
599           INTEGER, INTENT(IN) :: M
600     ! velocity gradient of solids phase M
601     !         |  du/dx    du/dy   du/dz  |
602     ! DELV =  |  dv/dx    dv/dy   dv/dz  |  =  dUi/dxj
603     !         |  dw/dx    dw/dy   dw/dz  |
604           DOUBLE PRECISION, INTENT(OUT) :: lVelGradS(3,3)
605     ! rate of strain tensor solids phase
606           DOUBLE PRECISION, INTENT(OUT) :: lRateStrainS(3,3)
607     
608     ! Local variables
609     !-----------------------------------------------
610     ! Indices
611           INTEGER :: I, J, K, IMJK, IJMK, IJKM
612     ! loop counters
613           INTEGER :: P, Q
614     !
615           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
616           DOUBLE PRECISION :: dudx,dudy,dudz
617           DOUBLE PRECISION :: dvdx,dvdy,dvdz
618           DOUBLE PRECISION :: dwdx,dwdy,dwdz
619           DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Vi,Wi,Sx,Sy,Sz
620           DOUBLE PRECISION :: UW_s,VW_s,WW_s
621     
622           LOGICAL :: U_NODE_AT_E, U_NODE_AT_W
623           LOGICAL :: V_NODE_AT_N, V_NODE_AT_S
624           LOGICAL :: W_NODE_AT_T, W_NODE_AT_B
625           INTEGER :: BCV
626           INTEGER :: BCT
627     !-----------------------------------------------
628     
629     !!$omp  parallel do private(IJK, I, J, K, IMJK, IJMK, IJKM) &
630     !!!!$omp& schedule(dynamic,chunk_size)
631     
632     ! initialize
633           lVelGradS = ZERO
634           lRateStrainS = ZERO
635     
636           I = I_OF(IJK)
637           J = J_OF(IJK)
638           K = K_OF(IJK)
639           IMJK = IM_OF(IJK)
640           IJMK = JM_OF(IJK)
641           IJKM = KM_OF(IJK)
642     
643           IF(FLOW_AT(IJK)) THEN
644              lVelGradS = ZERO
645              RETURN
646           ENDIF
647     
648           BCV = BC_ID(IJK)
649           IF(BCV > 0 ) THEN
650              BCT = BC_TYPE_ENUM(BCV)
651           ELSE
652              BCT = NONE
653           ENDIF
654           SELECT CASE (BCT)
655              CASE (CG_NSW)
656                 NOC_TRDS = .TRUE.
657                 UW_s = ZERO
658                 VW_s = ZERO
659                 WW_s = ZERO
660              CASE (CG_FSW)
661                 NOC_TRDS = .FALSE.
662                 UW_s = ZERO
663                 VW_s = ZERO
664                 WW_s = ZERO
665                 RETURN
666              CASE(CG_PSW)
667                 IF(BC_JJ_PS(BCV)==1) THEN   ! Johnson-Jackson partial slip bc
668                    NOC_TRDS = .FALSE.
669                    UW_s = BC_UW_S(BCV,M)
670                    VW_s = BC_VW_S(BCV,M)
671                    WW_s = BC_WW_S(BCV,M)
672                 ELSEIF(BC_HW_S(BCV,M)==UNDEFINED) THEN   ! same as NSW
673                    NOC_TRDS = .TRUE.
674                    UW_s = BC_UW_S(BCV,M)
675                    VW_s = BC_VW_S(BCV,M)
676                    WW_s = BC_WW_S(BCV,M)
677                 ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN   ! same as FSW
678                    NOC_TRDS = .FALSE.
679                    UW_s = ZERO
680                    VW_s = ZERO
681                    WW_s = ZERO
682                 ELSE                              ! partial slip
683                    NOC_TRDS = .FALSE.
684                 ENDIF
685              CASE (CG_MI)
686                 lVelGradS = ZERO
687                 RETURN
688              CASE (CG_PO)
689                 lVelGradS = ZERO
690                 RETURN
691              CASE (NONE)
692                 lVelGradS = ZERO
693                 RETURN
694           END SELECT
695     
696     
697     
698     ! du/dx, du/dy, du/dz
699     !=======================================================================
700           U_NODE_AT_E = ((.NOT.BLOCKED_U_CELL_AT(IJK)) .AND.&
701                          (.NOT.WALL_U_AT(IJK)))
702           U_NODE_AT_W = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.&
703                          (.NOT.WALL_U_AT(IMJK)))
704     
705           IF(U_NODE_AT_E.AND.U_NODE_AT_W) THEN
706              Ui = HALF * (U_S(IJK,M) + U_S(IMJK,M))
707              Xi = HALF * (X_U(IJK) + X_U(IMJK))
708              Yi = HALF * (Y_U(IJK) + Y_U(IMJK))
709              Zi = HALF * (Z_U(IJK) + Z_U(IMJK))
710              Sx = X_U(IJK) - X_U(IMJK)
711              Sy = Y_U(IJK) - Y_U(IMJK)
712              Sz = Z_U(IJK) - Z_U(IMJK)
713              CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
714              IF(abs(Sx) > ZERO) THEN
715                 dudx =  (U_S(IJK,M) - U_S(IMJK,M))/Sx
716                 dudy =  ZERO
717                 dudz =  ZERO
718                 IF(NOC_TRDS) THEN
719                    dudx = dudx - ((Ui-UW_s)/(Sx*DEL_H)*(Sy*Ny+Sz*Nz))
720                    dudy = (Ui-UW_s) / DEL_H * Ny
721                    dudz = (Ui-UW_s) / DEL_H * Nz
722                 ENDIF
723              ELSE
724                 dudx = ZERO
725                 dudy = ZERO
726                 dudz = ZERO
727              ENDIF
728           ELSEIF (U_NODE_AT_E.AND.(.NOT.U_NODE_AT_W).AND.NOC_TRDS) THEN
729              CALL GET_DEL_H(IJK,'SCALAR',X_U(IJK),Y_U(IJK),Z_U(IJK),&
730                             DEL_H,Nx,Ny,Nz)
731              dudx = (U_s(IJK,M) - UW_s) / DEL_H * Nx
732              dudy = (U_s(IJK,M) - UW_s) / DEL_H * Ny
733              dudz = (U_s(IJK,M) - UW_s) / DEL_H * Nz
734           ELSEIF ((.NOT.U_NODE_AT_E).AND.U_NODE_AT_W.AND.NOC_TRDS) THEN
735              CALL GET_DEL_H(IJK,'SCALAR',X_U(IMJK),Y_U(IMJK),Z_U(IMJK),&
736                             DEL_H,Nx,Ny,Nz)
737              dudx = (U_s(IMJK,M) - UW_s) / DEL_H * Nx
738              dudy = (U_s(IMJK,M) - UW_s) / DEL_H * Ny
739              dudz = (U_s(IMJK,M) - UW_s) / DEL_H * Nz
740           ELSE
741              dudx = ZERO
742              dudy = ZERO
743              dudz = ZERO
744           ENDIF
745           lVelGradS(1,1) = dudx
746           lVelGradS(1,2) = dudy
747           lVelGradS(1,3) = dudz
748     
749     ! dv/dx, dv/dy, dv/dz
750     !=======================================================================
751           V_NODE_AT_N = ((.NOT.BLOCKED_V_CELL_AT(IJK)) .AND.&
752                          (.NOT.WALL_V_AT(IJK)))
753           V_NODE_AT_S = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.&
754                          (.NOT.WALL_V_AT(IJMK)))
755     
756           IF(V_NODE_AT_N.AND.V_NODE_AT_S) THEN
757              Vi = HALF * (V_S(IJK,M) + V_S(IJMK,M))
758              Xi = HALF * (X_V(IJK) + X_V(IJMK))
759              Yi = HALF * (Y_V(IJK) + Y_V(IJMK))
760              Zi = HALF * (Z_V(IJK) + Z_V(IJMK))
761              Sx = X_V(IJK) - X_V(IJMK)
762              Sy = Y_V(IJK) - Y_V(IJMK)
763              Sz = Z_V(IJK) - Z_V(IJMK)
764              CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
765              IF(abs(Sy) > ZERO) THEN
766                 dvdx =  ZERO
767                 dvdy =  (V_S(IJK,M) - V_S(IJMK,M))/Sy
768                 dvdz =  ZERO
769                 IF(NOC_TRDS) THEN
770                    dvdx = (Vi-VW_s) / DEL_H * Nx
771                    dvdy = dvdy - ((Vi-VW_s)/(Sy*DEL_H)*(Sx*Nx+Sz*Nz))
772                    dvdz = (Vi-VW_s) / DEL_H * Nz
773                 ENDIF
774              ELSE
775                 dvdx =  ZERO
776                 dvdy =  ZERO
777                 dvdz =  ZERO
778              ENDIF
779           ELSEIF (V_NODE_AT_N.AND.(.NOT.V_NODE_AT_S).AND.NOC_TRDS) THEN
780              CALL GET_DEL_H(IJK,'SCALAR',X_V(IJK),Y_V(IJK),Z_V(IJK),&
781                             DEL_H,Nx,Ny,Nz)
782              dvdx = (V_s(IJK,M) - VW_s) / DEL_H * Nx
783              dvdy = (V_s(IJK,M) - VW_s) / DEL_H * Ny
784              dvdz = (V_s(IJK,M) - VW_s) / DEL_H * Nz
785           ELSEIF ((.NOT.V_NODE_AT_N).AND.V_NODE_AT_S.AND.NOC_TRDS) THEN
786              CALL GET_DEL_H(IJK,'SCALAR',X_V(IJMK),Y_V(IJMK),Z_V(IJMK),&
787                             DEL_H,Nx,Ny,Nz)
788              dvdx = (V_s(IJMK,M) - VW_s) / DEL_H * Nx
789              dvdy = (V_s(IJMK,M) - VW_s) / DEL_H * Ny
790              dvdz = (V_s(IJMK,M) - VW_s) / DEL_H * Nz
791           ELSE
792              dvdx =  ZERO
793              dvdy =  ZERO
794              dvdz =  ZERO
795           ENDIF
796           lVelGradS(2,1) = dvdx
797           lVelGradS(2,2) = dvdy
798           lVelGradS(2,3) = dvdz
799     
800     ! dw/dx, dw/dy, dw/dz
801     !=======================================================================
802           IF(NO_K) THEN
803              dwdx = ZERO
804              dwdy = ZERO
805              dwdz = ZERO
806           ELSE
807              W_NODE_AT_T = ((.NOT.BLOCKED_W_CELL_AT(IJK)) .AND.&
808                             (.NOT.WALL_W_AT(IJK)))
809              W_NODE_AT_B = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.&
810                             (.NOT.WALL_W_AT(IJKM)))
811     
812              IF(W_NODE_AT_T.AND.W_NODE_AT_B) THEN
813                 Wi = HALF * (W_S(IJK,M) + W_S(IJKM,M))
814                 Xi = HALF * (X_W(IJK) + X_W(IJKM))
815                 Yi = HALF * (Y_W(IJK) + Y_W(IJKM))
816                 Zi = HALF * (Z_W(IJK) + Z_W(IJKM))
817                 Sx = X_W(IJK) - X_W(IJKM)
818                 Sy = Y_W(IJK) - Y_W(IJKM)
819                 Sz = Z_W(IJK) - Z_W(IJKM)
820                 CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
821                 IF(abs(Sz) > ZERO) THEN
822                    dwdx =  ZERO
823                    dwdy =  ZERO
824                    dwdz =  (W_S(IJK,M) - W_S(IJKM,M))/Sz
825                    IF(NOC_TRDS) THEN
826                       dwdx = (Wi-WW_s) / DEL_H * Nx
827                       dwdy = (Wi-WW_s) / DEL_H * Ny
828                       dwdz = dwdz - ((Wi-WW_s)/(Sz*DEL_H)*(Sx*Nx+Sy*Ny))
829                    ENDIF
830                 ELSE
831                    dwdx = ZERO
832                    dwdy = ZERO
833                    dwdz = ZERO
834                 ENDIF
835              ELSEIF (W_NODE_AT_T.AND.(.NOT.W_NODE_AT_B).AND.NOC_TRDS) THEN
836                 CALL GET_DEL_H(IJK,'SCALAR',X_W(IJK),Y_W(IJK),Z_W(IJK),&
837                                DEL_H,Nx,Ny,Nz)
838                 dwdx = (W_s(IJK,M) - WW_s) / DEL_H * Nx
839                 dwdy = (W_s(IJK,M) - WW_s) / DEL_H * Ny
840                 dwdz = (W_s(IJK,M) - WW_s) / DEL_H * Nz
841              ELSEIF ((.NOT.W_NODE_AT_T).AND.W_NODE_AT_B.AND.NOC_TRDS) THEN
842                 CALL GET_DEL_H(IJK,'SCALAR',X_W(IJKM),Y_W(IJKM),Z_W(IJKM),&
843                                DEL_H,Nx,Ny,Nz)
844                 dwdx = (W_s(IJKM,M) - WW_s) / DEL_H * Nx
845                 dwdy = (W_s(IJKM,M) - WW_s) / DEL_H * Ny
846                 dwdz = (W_s(IJKM,M) - WW_s) / DEL_H * Nz
847              ELSE
848                 dwdx = ZERO
849                 dwdy = ZERO
850                 dwdz = ZERO
851              ENDIF
852           ENDIF
853           lVelGradS(3,1) = dwdx
854           lVelGradS(3,2) = dwdy
855           lVelGradS(3,3) = dwdz
856     
857           DO P = 1,3
858              DO Q = 1,3
859                 lRateStrainS(P,Q) = HALF * (lVelGradS(P,Q)+lVelGradS(Q,P))
860              ENDDO
861           ENDDO
862     
863     
864           RETURN
865           END SUBROUTINE CALC_CG_DERIV_VEL_SOLIDS
866