File: /nfs/home/0/users/jenkins/mfix.git/model/calc_trd_g.f

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