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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CALC_trD_s(trD_s, IER)                                 C
4     !  Purpose: Calculate the trace of solids phase rate of strain tensor  C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 19-DEC-96  C
7     !  Reviewer:                                          Date: dd-mmm-yy  C
8     !                                                                      C
9     !  Revision Number: 1                                                  C
10     !  Purpose: To incorporate Cartesian grid modifications                C
11     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
12     !                                                                      C
13     !  Literature/Document References:                                     C
14     !                                                                      C
15     !  Variables referenced:                                               C
16     !  Variables modified:                                                 C
17     !                                                                      C
18     !  Local variables:                                                    C
19     !                                                                      C
20     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
21     !
22           SUBROUTINE CALC_TRD_S(TRD_S)
23     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
24     !...Switches: -xf
25     !
26     !     Include param.inc file to specify parameter values
27     !
28     !-----------------------------------------------
29     !   M o d u l e s
30     !-----------------------------------------------
31           USE param
32           USE param1
33           USE parallel
34           USE geometry
35           USE fldvar
36           USE indices
37           USE physprop
38           USE compar
39           USE sendrecv
40           USE functions
41     !=======================================================================
42     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
43     !=======================================================================
44           USE bc
45           USE cutcell
46           USE quadric
47     !=======================================================================
48     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
49     !=======================================================================
50           IMPLICIT NONE
51     !-----------------------------------------------
52     !   G l o b a l   P a r a m e t e r s
53     !-----------------------------------------------
54     !-----------------------------------------------
55     !   D u m m y   A r g u m e n t s
56     !-----------------------------------------------
57     !
58     !                      Indices
59           INTEGER          I, J, K, IJK, IMJK, IJMK, IJKM, IM, M
60     !
61     !                      Strain rate tensor components for mth solids phase
62           DOUBLE PRECISION trD_s(DIMENSION_3, DIMENSION_M)
63     !=======================================================================
64     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
65     !=======================================================================
66           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
67           DOUBLE PRECISION :: dudx,dvdy,dwdz
68           DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Vi,Wi,Sx,Sy,Sz
69           DOUBLE PRECISION :: UW_s,VW_s,WW_s
70     
71           LOGICAL :: U_NODE_AT_E, U_NODE_AT_W
72           LOGICAL :: V_NODE_AT_N, V_NODE_AT_S
73           LOGICAL :: W_NODE_AT_T, W_NODE_AT_B
74           INTEGER :: BCV
75           CHARACTER(LEN=9) :: BCT
76     !=======================================================================
77     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
78     !=======================================================================
79     
80     !-----------------------------------------------
81           DO M = 1, MMAX
82     !!!!$omp    parallel do private(ijk,i,j,k,im,imjk,ijmk,ijkm)
83              DO IJK = ijkstart3, ijkend3
84                 IF (.NOT.WALL_AT(IJK)) THEN
85                    I = I_OF(IJK)
86                    J = J_OF(IJK)
87                    K = K_OF(IJK)
88                    IM = IM1(I)
89                    IMJK = IM_OF(IJK)
90                    IJMK = JM_OF(IJK)
91                    IJKM = KM_OF(IJK)
92     
93     !=======================================================================
94     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
95     !=======================================================================
96                    IF(.NOT.CUT_CELL_AT(IJK)) THEN
97     
98                       TRD_S(IJK,M) = (X_E(I)*U_S(IJK,M)-X_E(IM)*U_S(IMJK,M))*OX(I)*ODX&
99                          (I) + (V_S(IJK,M)-V_S(IJMK,M))*ODY(J) + (W_S(IJK,M)-W_S(IJKM,&
100                          M))*(OX(I)*ODZ(K))
101     
102                    ELSE  ! CUT CELL
103     
104                       BCV = BC_ID(IJK)
105     
106                       IF(BCV > 0 ) THEN
107                          BCT = BC_TYPE(BCV)
108                       ELSE
109                          BCT = 'NONE'
110                       ENDIF
111     
112                       SELECT CASE (BCT)
113                          CASE ('CG_NSW')
114                             NOC_TRDS = .TRUE.
115                             UW_s = ZERO
116                             VW_s = ZERO
117                             WW_s = ZERO
118                          CASE ('CG_FSW')
119                             NOC_TRDS = .FALSE.
120                             UW_s = ZERO
121                             VW_s = ZERO
122                             WW_s = ZERO
123                          CASE('CG_PSW')
124     
125                             IF(BC_JJ_PS(BCV)==1) THEN   ! Johnson-Jackson partial slip bc
126                                NOC_TRDS = .FALSE.
127                                UW_s = BC_UW_S(BCV,M)
128                                VW_s = BC_VW_S(BCV,M)
129                                WW_s = BC_WW_S(BCV,M)
130     
131                             ELSEIF(BC_HW_S(BCV,M)==UNDEFINED) THEN   ! same as NSW
132                                NOC_TRDS = .TRUE.
133                                UW_s = BC_UW_S(BCV,M)
134                                VW_s = BC_VW_S(BCV,M)
135                                WW_s = BC_WW_S(BCV,M)
136                             ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN    ! same as FSW
137                                NOC_TRDS = .FALSE.
138                             ELSE                                 ! partial slip
139                                NOC_TRDS = .FALSE.
140                                UW_s = ZERO
141                                VW_s = ZERO
142                                WW_s = ZERO
143                             ENDIF
144                          CASE ('CG_MI')
145                              TRD_S(IJK,M) = ZERO
146                             RETURN
147                          CASE ('CG_PO')
148                              TRD_S(IJK,M) = ZERO
149                             RETURN
150                       END SELECT
151     
152                       IF(FLOW_AT(IJK)) THEN
153                           TRD_S(IJK,M) = ZERO
154                          RETURN
155                       ENDIF
156     
157     !              du/dx
158     
159                       U_NODE_AT_E = ((.NOT.BLOCKED_U_CELL_AT(IJK)) .AND.(.NOT.WALL_U_AT(IJK)))
160                       U_NODE_AT_W = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.(.NOT.WALL_U_AT(IMJK)))
161     
162                       IF(U_NODE_AT_E.AND.U_NODE_AT_W) THEN
163     
164                          Ui = HALF * (U_S(IJK,M) + U_S(IMJK,M))
165                          Xi = HALF * (X_U(IJK) + X_U(IMJK))
166                          Yi = HALF * (Y_U(IJK) + Y_U(IMJK))
167                          Zi = HALF * (Z_U(IJK) + Z_U(IMJK))
168                          Sx = X_U(IJK) - X_U(IMJK)
169                          Sy = Y_U(IJK) - Y_U(IMJK)
170                          Sz = Z_U(IJK) - Z_U(IMJK)
171     
172                          CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
173     
174                          IF(Sx /= ZERO) THEN
175                             dudx =  (U_S(IJK,M) - U_S(IMJK,M))/Sx
176                             IF(NOC_TRDS) dudx = dudx - ((Ui-UW_s)/(Sx*DEL_H)*(Sy*Ny+Sz*Nz))
177                          ELSE
178                             dudx = ZERO
179                          ENDIF
180     
181                       ELSE
182                          dudx = ZERO
183                       ENDIF
184     
185     !              dv/dy
186     
187                       V_NODE_AT_N = ((.NOT.BLOCKED_V_CELL_AT(IJK)) .AND.(.NOT.WALL_V_AT(IJK)))
188                       V_NODE_AT_S = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.(.NOT.WALL_V_AT(IJMK)))
189     
190                       IF(V_NODE_AT_N.AND.V_NODE_AT_S) THEN
191     
192                          Vi = HALF * (V_S(IJK,M) + V_S(IJMK,M))
193                          Xi = HALF * (X_V(IJK) + X_V(IJMK))
194                          Yi = HALF * (Y_V(IJK) + Y_V(IJMK))
195                          Zi = HALF * (Z_V(IJK) + Z_V(IJMK))
196                          Sx = X_V(IJK) - X_V(IJMK)
197                          Sy = Y_V(IJK) - Y_V(IJMK)
198                          Sz = Z_V(IJK) - Z_V(IJMK)
199     
200                          CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
201     
202                          IF(Sy /= ZERO) THEN
203                             dvdy =  (V_S(IJK,M) - V_S(IJMK,M))/Sy
204                             IF(NOC_TRDS) dvdy = dvdy - ((Vi-VW_s)/(Sy*DEL_H)*(Sx*Nx+Sz*Nz))
205                          ELSE
206                             dvdy =  ZERO
207                          ENDIF
208     
209                       ELSE
210                          dvdy = ZERO
211                       ENDIF
212     
213     !              dw/dz
214     
215                       IF(NO_K) THEN
216     
217                          dwdz = ZERO
218     
219                       ELSE
220     
221                          W_NODE_AT_T = ((.NOT.BLOCKED_W_CELL_AT(IJK)) .AND.(.NOT.WALL_W_AT(IJK)))
222                          W_NODE_AT_B = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.(.NOT.WALL_W_AT(IJKM)))
223     
224                          IF(W_NODE_AT_T.AND.W_NODE_AT_B) THEN
225     
226                             Wi = HALF * (W_S(IJK,M) + W_S(IJKM,M))
227                             Xi = HALF * (X_W(IJK) + X_W(IJKM))
228                             Yi = HALF * (Y_W(IJK) + Y_W(IJKM))
229                             Zi = HALF * (Z_W(IJK) + Z_W(IJKM))
230                             Sx = X_W(IJK) - X_W(IJKM)
231                             Sy = Y_W(IJK) - Y_W(IJKM)
232                             Sz = Z_W(IJK) - Z_W(IJKM)
233     
234                             CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
235     
236                             IF(Sz /= ZERO) THEN
237                                dwdz =  (W_S(IJK,M) - W_S(IJKM,M))/Sz
238                                IF(NOC_TRDS) dwdz = dwdz - ((Wi-WW_s)/(Sz*DEL_H)*(Sx*Nx+Sy*Ny))
239                             ELSE
240                                dwdz = ZERO
241                             ENDIF
242     
243                          ELSE
244                             dwdz = ZERO
245                          ENDIF
246     
247                       ENDIF  ! NO_K
248     
249                       TRD_S(IJK,M) = dudx + dvdy + dwdz
250     
251                    ENDIF   ! CUT CELL
252     
253     ! Original term:
254     !               TRD_S(IJK,M) = (X_E(I)*U_S(IJK,M)-X_E(IM)*U_S(IMJK,M))*OX(I)*ODX&
255     !                  (I) + (V_S(IJK,M)-V_S(IJMK,M))*ODY(J) + (W_S(IJK,M)-W_S(IJKM,&
256     !                  M))*(OX(I)*ODZ(K))
257     !=======================================================================
258     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
259     !=======================================================================
260                 ENDIF
261              END DO
262           END DO
263     
264           END SUBROUTINE CALC_TRD_S
265     
266     !// Comments on the modifications for DMP version implementation
267     !// 001 Include header file and common declarations for parallelization
268     !// 350 Changed do loop limits: 1,ijkmax2-> ijkstart3, ijkend3
269     
270     
271     
272     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
273     !                                                                      C
274     !  Module name: CG_CALC_VEL_S_GRAD(IJK,M,D_s, IER)                     C
275     !  Purpose: Calculate velocity derivatives in scalar cut-cell          C
276     !           Solids phase                                               C
277     !                                                                      C
278     !  Author: Jeff Dietiker                              Date: 25-JAN-96  C
279     !  Reviewer:                                          Date: dd-mmm-yy  C
280     !                                                                      C
281     !                                                                      C
282     !  Literature/Document References:                                     C
283     !                                                                      C
284     !  Variables referenced:                                               C
285     !  Variables modified:                                                 C
286     !                                                                      C
287     !  Local variables:                                                    C
288     !                                                                      C
289     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
290     !
291           SUBROUTINE CG_CALC_VEL_S_GRAD(IJK,M,D_s)
292     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
293     !...Switches: -xf
294     !
295     !     Include param.inc file to specify parameter values
296     !
297     !-----------------------------------------------
298     !   M o d u l e s
299     !-----------------------------------------------
300           USE param
301           USE param1
302           USE parallel
303           USE geometry
304           USE fldvar
305           USE indices
306           USE compar
307           USE sendrecv
308           USE bc
309           USE cutcell
310           USE quadric
311           USE functions
312           IMPLICIT NONE
313     !-----------------------------------------------
314     !   G l o b a l   P a r a m e t e r s
315     !-----------------------------------------------
316     !-----------------------------------------------
317     !   D u m m y   A r g u m e n t s
318     !-----------------------------------------------
319     !
320     !                      Indices
321           INTEGER          I, J, K, IJK, IMJK, IJMK, IJKM, M
322     !
323           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
324           DOUBLE PRECISION :: dudx,dudy,dudz
325           DOUBLE PRECISION :: dvdx,dvdy,dvdz
326           DOUBLE PRECISION :: dwdx,dwdy,dwdz
327           DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Vi,Wi,Sx,Sy,Sz
328           DOUBLE PRECISION :: UW_s,VW_s,WW_s
329           DOUBLE PRECISION, DIMENSION (3,3) :: DELV
330           INTEGER :: P,Q
331           DOUBLE PRECISION, DIMENSION (3,3) :: D_s
332     
333     
334     !              |  du/dx    du/dy   du/dz  |
335     !      DELV =  |  dv/dx    dv/dy   dv/dz  |  =  dUi/dxj
336     !              |  dw/dx    dw/dy   dw/dz  |
337     
338           LOGICAL :: U_NODE_AT_E, U_NODE_AT_W
339           LOGICAL :: V_NODE_AT_N, V_NODE_AT_S
340           LOGICAL :: W_NODE_AT_T, W_NODE_AT_B
341           INTEGER :: BCV
342           CHARACTER(LEN=9) :: BCT
343     !-----------------------------------------------
344     !
345     !
346     !!!!$omp  parallel do private( IJK, I,J,K, IM,IMJK,IJMK,IJKM ) &
347     !!!!$omp& schedule(dynamic,chunk_size)
348     
349           DELV = ZERO
350           D_s  = ZERO
351     
352           IF (.NOT.CUT_CELL_AT(IJK)) RETURN
353     
354           I = I_OF(IJK)
355           J = J_OF(IJK)
356           K = K_OF(IJK)
357     
358           IMJK = IM_OF(IJK)
359           IJMK = JM_OF(IJK)
360           IJKM = KM_OF(IJK)
361     
362           BCV = BC_ID(IJK)
363     
364           IF(BCV > 0 ) THEN
365              BCT = BC_TYPE(BCV)
366           ELSE
367              BCT = 'NONE'
368           ENDIF
369     
370           SELECT CASE (BCT)
371              CASE ('CG_NSW')
372                 NOC_TRDS = .TRUE.
373                 UW_s = ZERO
374                 VW_s = ZERO
375                 WW_s = ZERO
376              CASE ('CG_FSW')
377                 NOC_TRDS = .FALSE.
378                 UW_s = ZERO
379                 VW_s = ZERO
380                 WW_s = ZERO
381                 RETURN
382     
383              CASE('CG_PSW')
384     
385                 IF(BC_JJ_PS(BCV)==1) THEN   ! Johnson-Jackson partial slip bc
386                    NOC_TRDS = .FALSE.
387                    UW_s = BC_UW_S(BCV,M)
388                    VW_s = BC_VW_S(BCV,M)
389                    WW_s = BC_WW_S(BCV,M)
390     
391                 ELSEIF(BC_HW_S(BCV,M)==UNDEFINED) THEN   ! same as NSW
392                    NOC_TRDS = .TRUE.
393                    UW_s = BC_UW_S(BCV,M)
394                    VW_s = BC_VW_S(BCV,M)
395                    WW_s = BC_WW_S(BCV,M)
396                 ELSEIF(BC_HW_S(BCV,M)==ZERO) THEN   ! same as FSW
397                    NOC_TRDS = .FALSE.
398                    UW_s = ZERO
399                    VW_s = ZERO
400                    WW_s = ZERO
401                 ELSE                              ! partial slip
402                    NOC_TRDS = .FALSE.
403                 ENDIF
404              CASE ('CG_MI')
405                 DELV = ZERO
406                 RETURN
407              CASE ('CG_PO')
408                 DELV = ZERO
409                 RETURN
410              CASE ('NONE')
411                 DELV = ZERO
412                 RETURN
413           END SELECT
414     
415           IF(FLOW_AT(IJK)) THEN
416              DELV = ZERO
417              RETURN
418           ENDIF
419     
420     !=======================================================================
421     !              du/dx, du/dy, du/dz
422     !=======================================================================
423     
424           U_NODE_AT_E = ((.NOT.BLOCKED_U_CELL_AT(IJK)) .AND.(.NOT.WALL_U_AT(IJK)))
425           U_NODE_AT_W = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.(.NOT.WALL_U_AT(IMJK)))
426     
427           IF(U_NODE_AT_E.AND.U_NODE_AT_W) THEN
428     
429              Ui = HALF * (U_S(IJK,M) + U_S(IMJK,M))
430              Xi = HALF * (X_U(IJK) + X_U(IMJK))
431              Yi = HALF * (Y_U(IJK) + Y_U(IMJK))
432              Zi = HALF * (Z_U(IJK) + Z_U(IMJK))
433              Sx = X_U(IJK) - X_U(IMJK)
434              Sy = Y_U(IJK) - Y_U(IMJK)
435              Sz = Z_U(IJK) - Z_U(IMJK)
436     
437              CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
438     
439              IF(abs(Sx) > ZERO) THEN
440                 dudx =  (U_S(IJK,M) - U_S(IMJK,M))/Sx
441                 dudy =  ZERO
442                 dudz =  ZERO
443                 IF(NOC_TRDS) THEN
444                    dudx = dudx - ((Ui-UW_s)/(Sx*DEL_H)*(Sy*Ny+Sz*Nz))
445                    dudy = (Ui-UW_s) / DEL_H * Ny
446                    dudz = (Ui-UW_s) / DEL_H * Nz
447                 ENDIF
448              ELSE
449                 dudx = ZERO
450                 dudy = ZERO
451                 dudz = ZERO
452              ENDIF
453     
454     
455           ELSE IF (U_NODE_AT_E.AND.(.NOT.U_NODE_AT_W).AND.NOC_TRDS) THEN
456              CALL GET_DEL_H(IJK,'SCALAR',X_U(IJK),Y_U(IJK),Z_U(IJK),DEL_H,Nx,Ny,Nz)
457              dudx = (U_s(IJK,M) - UW_s) / DEL_H * Nx
458              dudy = (U_s(IJK,M) - UW_s) / DEL_H * Ny
459              dudz = (U_s(IJK,M) - UW_s) / DEL_H * Nz
460           ELSE IF ((.NOT.U_NODE_AT_E).AND.U_NODE_AT_W.AND.NOC_TRDS) THEN
461              CALL GET_DEL_H(IJK,'SCALAR',X_U(IMJK),Y_U(IMJK),Z_U(IMJK),DEL_H,Nx,Ny,Nz)
462              dudx = (U_s(IMJK,M) - UW_s) / DEL_H * Nx
463              dudy = (U_s(IMJK,M) - UW_s) / DEL_H * Ny
464              dudz = (U_s(IMJK,M) - UW_s) / DEL_H * Nz
465           ELSE
466              dudx = ZERO
467              dudy = ZERO
468              dudz = ZERO
469           ENDIF
470     
471     
472           DELV(1,1) = dudx
473           DELV(1,2) = dudy
474           DELV(1,3) = dudz
475     !=======================================================================
476     !              dv/dx, dv/dy, dv/dz
477     !=======================================================================
478     
479           V_NODE_AT_N = ((.NOT.BLOCKED_V_CELL_AT(IJK)) .AND.(.NOT.WALL_V_AT(IJK)))
480           V_NODE_AT_S = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.(.NOT.WALL_V_AT(IJMK)))
481     
482           IF(V_NODE_AT_N.AND.V_NODE_AT_S) THEN
483     
484              Vi = HALF * (V_S(IJK,M) + V_S(IJMK,M))
485              Xi = HALF * (X_V(IJK) + X_V(IJMK))
486              Yi = HALF * (Y_V(IJK) + Y_V(IJMK))
487              Zi = HALF * (Z_V(IJK) + Z_V(IJMK))
488              Sx = X_V(IJK) - X_V(IJMK)
489              Sy = Y_V(IJK) - Y_V(IJMK)
490              Sz = Z_V(IJK) - Z_V(IJMK)
491     
492              CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
493     
494              IF(abs(Sy) > ZERO) THEN
495                 dvdx =  ZERO
496                 dvdy =  (V_S(IJK,M) - V_S(IJMK,M))/Sy
497                 dvdz =  ZERO
498                 IF(NOC_TRDS) THEN
499                    dvdx = (Vi-VW_s) / DEL_H * Nx
500                    dvdy = dvdy - ((Vi-VW_s)/(Sy*DEL_H)*(Sx*Nx+Sz*Nz))
501                    dvdz = (Vi-VW_s) / DEL_H * Nz
502                 ENDIF
503              ELSE
504                 dvdx =  ZERO
505                 dvdy =  ZERO
506                 dvdz =  ZERO
507              ENDIF
508     
509           ELSE IF (V_NODE_AT_N.AND.(.NOT.V_NODE_AT_S).AND.NOC_TRDS) THEN
510              CALL GET_DEL_H(IJK,'SCALAR',X_V(IJK),Y_V(IJK),Z_V(IJK),DEL_H,Nx,Ny,Nz)
511              dvdx = (V_s(IJK,M) - VW_s) / DEL_H * Nx
512              dvdy = (V_s(IJK,M) - VW_s) / DEL_H * Ny
513              dvdz = (V_s(IJK,M) - VW_s) / DEL_H * Nz
514           ELSE IF ((.NOT.V_NODE_AT_N).AND.V_NODE_AT_S.AND.NOC_TRDS) THEN
515              CALL GET_DEL_H(IJK,'SCALAR',X_V(IJMK),Y_V(IJMK),Z_V(IJMK),DEL_H,Nx,Ny,Nz)
516              dvdx = (V_s(IJMK,M) - VW_s) / DEL_H * Nx
517              dvdy = (V_s(IJMK,M) - VW_s) / DEL_H * Ny
518              dvdz = (V_s(IJMK,M) - VW_s) / DEL_H * Nz
519           ELSE
520              dvdx =  ZERO
521              dvdy =  ZERO
522              dvdz =  ZERO
523           ENDIF
524     
525     
526           DELV(2,1) = dvdx
527           DELV(2,2) = dvdy
528           DELV(2,3) = dvdz
529     
530     !=======================================================================
531     !              dw/dx, dw/dy, dw/dz
532     !=======================================================================
533     
534           IF(NO_K) THEN
535     
536              dwdx = ZERO
537              dwdy = ZERO
538              dwdz = ZERO
539     
540           ELSE
541     
542              W_NODE_AT_T = ((.NOT.BLOCKED_W_CELL_AT(IJK)) .AND.(.NOT.WALL_W_AT(IJK)))
543              W_NODE_AT_B = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.(.NOT.WALL_W_AT(IJKM)))
544     
545              IF(W_NODE_AT_T.AND.W_NODE_AT_B) THEN
546     
547                 Wi = HALF * (W_S(IJK,M) + W_S(IJKM,M))
548                 Xi = HALF * (X_W(IJK) + X_W(IJKM))
549                 Yi = HALF * (Y_W(IJK) + Y_W(IJKM))
550                 Zi = HALF * (Z_W(IJK) + Z_W(IJKM))
551                 Sx = X_W(IJK) - X_W(IJKM)
552                 Sy = Y_W(IJK) - Y_W(IJKM)
553                 Sz = Z_W(IJK) - Z_W(IJKM)
554     
555                 CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
556     
557                 IF(abs(Sz) > ZERO) THEN
558                    dwdx =  ZERO
559                    dwdy =  ZERO
560                    dwdz =  (W_S(IJK,M) - W_S(IJKM,M))/Sz
561                    IF(NOC_TRDS) THEN
562                       dwdx = (Wi-WW_s) / DEL_H * Nx
563                       dwdy = (Wi-WW_s) / DEL_H * Ny
564                       dwdz = dwdz - ((Wi-WW_s)/(Sz*DEL_H)*(Sx*Nx+Sy*Ny))
565                    ENDIF
566                 ELSE
567                    dwdx = ZERO
568                    dwdy = ZERO
569                    dwdz = ZERO
570                 ENDIF
571     
572              ELSE IF (W_NODE_AT_T.AND.(.NOT.W_NODE_AT_B).AND.NOC_TRDS) THEN
573                 CALL GET_DEL_H(IJK,'SCALAR',X_W(IJK),Y_W(IJK),Z_W(IJK),DEL_H,Nx,Ny,Nz)
574                 dwdx = (W_s(IJK,M) - WW_s) / DEL_H * Nx
575                 dwdy = (W_s(IJK,M) - WW_s) / DEL_H * Ny
576                 dwdz = (W_s(IJK,M) - WW_s) / DEL_H * Nz
577     
578              ELSE IF ((.NOT.W_NODE_AT_T).AND.W_NODE_AT_B.AND.NOC_TRDS) THEN
579                 CALL GET_DEL_H(IJK,'SCALAR',X_W(IJKM),Y_W(IJKM),Z_W(IJKM),DEL_H,Nx,Ny,Nz)
580                 dwdx = (W_s(IJKM,M) - WW_s) / DEL_H * Nx
581                 dwdy = (W_s(IJKM,M) - WW_s) / DEL_H * Ny
582                 dwdz = (W_s(IJKM,M) - WW_s) / DEL_H * Nz
583     
584              ELSE
585                 dwdx = ZERO
586                 dwdy = ZERO
587                 dwdz = ZERO
588              ENDIF
589     
590           ENDIF
591     
592           DELV(3,1) = dwdx
593           DELV(3,2) = dwdy
594           DELV(3,3) = dwdz
595     
596           DO P = 1,3
597              DO Q = 1,3
598                 D_s(P,Q) = HALF * (DELV(P,Q)+DELV(Q,P))
599              ENDDO
600           ENDDO
601     
602     
603           RETURN
604           END SUBROUTINE CG_CALC_VEL_S_GRAD
605     
606     !// Comments on the modifications for DMP version implementation
607     !// 001 Include header file and common declarations for parallelization
608     !// 350 Changed do loop limits: 1,ijkmax2-> ijkstart3, ijkend3
609     
610     
611     
612     
613     
614