MFIX  2016-1
conv_dif_v_g.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Module name: CONV_DIF_V_g(A_m, B_m, IER) C
4 ! Purpose: Determine convection diffusion terms for V_g momentum eqs C
5 ! The off-diagonal coefficients calculated here must be positive. The C
6 ! center coefficient and the source vector are negative; C
7 ! See source_v_g C
8 ! C
9 ! Author: M. Syamlal Date: 24-DEC-96 C
10 ! C
11 ! C
12 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
13  SUBROUTINE conv_dif_v_g(A_M, B_M, IER)
14 
15 ! Modules
16 !---------------------------------------------------------------------//
17  USE param, only: dimension_3, dimension_m
18  USE run, only: momentum_y_eq
19  USE run, only: discretize
20  USE run, only: def_cor
21  USE visc_g, only: epmu_gt
22  IMPLICIT NONE
23 
24 ! Dummy arguments
25 !---------------------------------------------------------------------//
26 ! Septadiagonal matrix A_m
27  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
28 ! Vector b_m
29  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
30 ! Error index
31  INTEGER, INTENT(INOUT) :: IER
32 !---------------------------------------------------------------------//
33 
34  IF (.NOT.momentum_y_eq(0)) RETURN
35 
36  IF (def_cor) THEN
37 ! USE DEFERRED CORRECTION TO SOLVE U_G
38  CALL store_a_v_g0 (a_m)
39  IF (discretize(4) > 1) CALL store_a_v_gdc(b_m(1,0))
40 
41  ELSE
42 ! DO NOT USE DEFERRED CORRECTION TO SOLVE V_G
43  IF (discretize(4) == 0) THEN ! 0 & 1 => FOUP
44  CALL store_a_v_g0(a_m)
45  ELSE
46  CALL store_a_v_g1(a_m)
47  ENDIF
48  ENDIF
49 
50  CALL dif_v_is(epmu_gt, a_m, 0, ier)
51 
52  RETURN
53  END SUBROUTINE conv_dif_v_g
54 
55 
56 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
57 ! C
58 ! Purpose: Calculate the components of velocity on the east, north, C
59 ! and top face of a v-momentum cell C
60 ! C
61 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
62  SUBROUTINE get_vcell_gvterms(U, V, WW)
63 
64 ! Modules
65 !---------------------------------------------------------------------//
66  USE compar, only: ijkstart3, ijkend3
67 
68  USE cutcell, only: cut_v_treatment_at
69  USE cutcell, only: theta_vn, theta_vn_bar
70  USE cutcell, only: theta_v_ne, theta_v_se
71  USE cutcell, only: theta_v_nt, theta_v_st
72  USE cutcell, only: alpha_ve_c, alpha_vn_c, alpha_vt_c
73 
74  USE fldvar, only: u_g, v_g, w_g
75 
76  USE fun_avg, only: avg_y_n, avg_y
77  USE functions, only: jp_of
78  USE geometry, only: do_k
79  USE indices, only: j_of
80 
81  USE param, only: dimension_3
82  IMPLICIT NONE
83 
84 ! Dummy arguments
85 !---------------------------------------------------------------------//
86  DOUBLE PRECISION, INTENT(OUT) :: U(dimension_3)
87  DOUBLE PRECISION, INTENT(OUT) :: V(dimension_3)
88  DOUBLE PRECISION, INTENT(OUT) :: WW(dimension_3)
89 
90 ! Local variables
91 !---------------------------------------------------------------------//
92 ! indices
93  INTEGER :: IJK, J, IJPK
94 ! for cartesian grid
95  DOUBLE PRECISION :: AW, HW, VELW
96 !---------------------------------------------------------------------//
97 
98 !!!$omp parallel do private(IJK,J,IJPK)
99  DO ijk = ijkstart3, ijkend3
100  j = j_of(ijk)
101  ijpk = jp_of(ijk)
102 
103  IF(cut_v_treatment_at(ijk)) THEN
104 
105 ! East face (i+1/2, j+1/2, k)
106  u(ijk) = (theta_v_se(ijk) * u_g(ijk) + &
107  theta_v_ne(ijk) * u_g(ijpk))
108  CALL get_interpolation_terms_g(ijk,'U_MOMENTUM', &
109  alpha_ve_c(ijk), aw, hw, velw)
110  u(ijk) = u(ijk) * aw
111 
112 ! North face (i, j+1, k)
113  v(ijk) = (theta_vn_bar(ijk) * v_g(ijk) + &
114  theta_vn(ijk) * v_g(ijpk))
115  CALL get_interpolation_terms_g(ijk,'U_MOMENTUM',&
116  alpha_vn_c(ijk), aw, hw, velw)
117  v(ijk) = v(ijk) * aw
118 
119 ! Top face (i, j+1/2, k+1/2)
120  IF (do_k) THEN
121  ww(ijk) = (theta_v_nt(ijk) * w_g(ijk) + &
122  theta_v_st(ijk) * w_g(ijpk))
123  CALL get_interpolation_terms_g(ijk,'U_MOMENTUM',&
124  alpha_vt_c(ijk), aw, hw, velw)
125  ww(ijk) = ww(ijk) * aw
126  ENDIF
127 
128  ELSE
129  u(ijk) = avg_y(u_g(ijk),u_g(ijpk),j)
130  v(ijk) = avg_y_n(v_g(ijk),v_g(ijpk))
131  IF (do_k) ww(ijk) = avg_y(w_g(ijk),w_g(ijpk),j)
132  ENDIF
133  ENDDO ! end do ijk
134 
135  RETURN
136  END SUBROUTINE get_vcell_gvterms
137 
138 
139 
140 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
141 ! C
142 ! Purpose: Calculate the convective fluxes through the faces of a C
143 ! v-momentum cell. Note the fluxes are calculated at all faces C
144 ! regardless of flow_at_n of condition of the west, south, or C
145 ! bottom face. C
146 ! C
147 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
148  SUBROUTINE get_vcell_gcflux_terms(FLUX_E, FLUX_W, FLUX_N, &
149  flux_s, flux_t, flux_b, ijk)
150 
151 ! Modules
152 !---------------------------------------------------------------------//
153  USE cutcell, only: cut_v_treatment_at
154  USE cutcell, only: theta_vn, theta_vn_bar
155  USE cutcell, only: theta_v_ne, theta_v_se
156  USE cutcell, only: theta_v_nt, theta_v_st
157  USE cutcell, only: alpha_ve_c, alpha_vn_c, alpha_vt_c
158 
159  USE functions, only: jp_of, im_of, jm_of, km_of
160 
161  USE geometry, only: do_k
162 
163  USE mflux, only: flux_ge, flux_gn, flux_gt
164 
165  USE param1, only: half
166  IMPLICIT NONE
167 
168 ! Dummy arguments
169 !---------------------------------------------------------------------//
170 ! fluxes through faces of given ijk u-momentum cell
171  DOUBLE PRECISION, INTENT(OUT) :: flux_e, flux_w
172  DOUBLE PRECISION, INTENT(OUT) :: flux_n, flux_s
173  DOUBLE PRECISION, INTENT(OUT) :: flux_t, flux_b
174 ! indices
175  INTEGER, INTENT(IN) :: ijk
176 
177 ! Local variables
178 !---------------------------------------------------------------------//
179  INTEGER :: imjk, ijmk, ijkm
180  INTEGER :: ijpk, imjpk, ijpkm
181 
182 ! for cartesian grid
183  DOUBLE PRECISION :: AW, HW, VELW
184 !---------------------------------------------------------------------//
185 ! indices
186  ijpk = jp_of(ijk)
187  imjk = im_of(ijk)
188  ijmk = jm_of(ijk)
189  ijkm = km_of(ijk)
190  ijpkm = jp_of(ijkm)
191  imjpk = jp_of(imjk)
192 
193  IF(cut_v_treatment_at(ijk)) THEN
194 ! East face (i+1/2, j+1/2, k)
195  flux_e = (theta_v_se(ijk) * flux_ge(ijk) + &
196  theta_v_ne(ijk) * flux_ge(ijpk))
197  CALL get_interpolation_terms_g(ijk,'V_MOMENTUM',&
198  alpha_ve_c(ijk), aw, hw, velw)
199 ! West face (i-1/2, j+1/2, k)
200  flux_e = flux_e * aw
201  flux_w = (theta_v_se(imjk) * flux_ge(imjk) + &
202  theta_v_ne(imjk) * flux_ge(imjpk))
203  CALL get_interpolation_terms_g(ijk,'V_MOMENTUM',&
204  alpha_ve_c(imjk), aw, hw, velw)
205  flux_w = flux_w * aw
206 
207 ! North face (i, j+1, k)
208  flux_n = (theta_vn_bar(ijk) * flux_gn(ijk) + &
209  theta_vn(ijk) * flux_gn(ijpk))
210  CALL get_interpolation_terms_g(ijk,'V_MOMENTUM',&
211  alpha_vn_c(ijk), aw, hw, velw)
212  flux_n = flux_n * aw
213 ! South face (i, j, k)
214  flux_s = (theta_vn_bar(ijmk) * flux_gn(ijmk) + &
215  theta_vn(ijmk) * flux_gn(ijk))
216  CALL get_interpolation_terms_g(ijk,'V_MOMENTUM',&
217  alpha_vn_c(ijmk), aw, hw, velw)
218  flux_s = flux_s * aw
219 
220  IF (do_k) THEN
221 ! Top face (i, j+1/2, k+1/2)
222  flux_t = (theta_v_nt(ijk) * flux_gt(ijk) + &
223  theta_v_st(ijk) * flux_gt(ijpk))
224  CALL get_interpolation_terms_g(ijk,'V_MOMENTUM',&
225  alpha_vt_c(ijk), aw, hw, velw)
226  flux_t = flux_t * aw
227 ! Bottom face (i, j+1/2, k-1/2)
228  flux_b = (theta_v_nt(ijkm) * flux_gt(ijkm) + &
229  theta_v_st(ijkm) * flux_gt(ijpkm))
230  CALL get_interpolation_terms_g(ijk,'V_MOMENTUM',&
231  alpha_vt_c(ijkm), aw, hw, velw)
232  flux_b = flux_b * aw
233  ENDIF
234  ELSE
235  flux_e = half * (flux_ge(ijk) + flux_ge(ijpk))
236  flux_w = half * (flux_ge(imjk) + flux_ge(imjpk))
237 
238  flux_n = half * (flux_gn(ijk) + flux_gn(ijpk))
239  flux_s = half * (flux_gn(ijmk) + flux_gn(ijk))
240 
241  IF (do_k) THEN
242  flux_t = half * (flux_gt(ijk) + flux_gt(ijpk))
243  flux_b = half * (flux_gt(ijkm) + flux_gt(ijpkm))
244  ENDIF
245  ENDIF ! end if/else cut_v_treatmeant_at
246 
247  RETURN
248  END SUBROUTINE get_vcell_gcflux_terms
249 
250 
251 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
252 ! C
253 ! Purpose: Calculate the components of diffusive flux through the C
254 ! the faces of a v-momentum cell. Note the fluxes are calculated at C
255 ! all faces regardless of flow_at_n of condition of the west, south C
256 ! or bottom face. C
257 ! C
258 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
259  SUBROUTINE get_vcell_gdiff_terms(D_FE, D_FW, D_FN, D_FS, &
260  d_ft, d_fb, ijk)
261 
262 ! Modules
263 !---------------------------------------------------------------------//
264  USE cutcell, only: cut_v_treatment_at
265  USE cutcell, only: oneodx_e_v, oneody_n_v, oneodz_t_v
266 
267  USE fldvar, only: epg_jfac
268 
269  USE functions, only: wall_at
270  USE functions, only: east_of, north_of, top_of
271  USE functions, only: west_of, bottom_of
272  USE functions, only: im_of, jm_of, km_of
273 
274  USE geometry, only: odx_e, ody, odz_t
275  USE geometry, only: do_k
276  USE geometry, only: ox
277  USE geometry, only: ayz_v, axz_v, axy_v
278 
279  USE indices, only: i_of, j_of, k_of
280  USE indices, only: jp1, im1, km1
281 
282  USE param
283  USE param1, only: zero
284  USE visc_g, only: epmu_gt, df_gv
285  IMPLICIT NONE
286 
287 ! Dummy arguments
288 !---------------------------------------------------------------------//
289 ! diffusion through faces of given ijk v-momentum cell
290  DOUBLE PRECISION, INTENT(OUT) :: d_fe, d_fw
291  DOUBLE PRECISION, INTENT(OUT) :: d_fn, d_fs
292  DOUBLE PRECISION, INTENT(OUT) :: d_ft, d_fb
293 ! ijk index
294  INTEGER, INTENT(IN) :: ijk
295 
296 ! Local variables
297 !---------------------------------------------------------------------//
298 ! indices
299  INTEGER :: imjk, ijmk, ijkm
300  INTEGER :: i, j, k, jp, im, km
301  INTEGER :: ijkc, ijkn, ijke, ijkne, ijkw, ijknw
302  INTEGER :: ijkt, ijktn, ijkb, ijkbn
303 ! length terms
304  DOUBLE PRECISION :: C_AE, C_AW, C_AN, C_AS, C_AT, C_AB
305 ! avg voidage
306  DOUBLE PRECISION :: EPGA
307 !---------------------------------------------------------------------//
308 
309  imjk = im_of(ijk)
310  ijmk = jm_of(ijk)
311  ijkm = km_of(ijk)
312 
313  i = i_of(ijk)
314  j = j_of(ijk)
315  k = k_of(ijk)
316  im = im1(i)
317  jp = jp1(j)
318  km = km1(k)
319 
320  ijkn = north_of(ijk)
321  IF (wall_at(ijk)) THEN
322  ijkc = ijkn
323  ELSE
324  ijkc = ijk
325  ENDIF
326  ijke = east_of(ijk)
327  ijkne = east_of(ijkn)
328  ijkw = west_of(ijk)
329  ijknw = north_of(ijkw)
330 
331  IF(cut_v_treatment_at(ijk)) THEN
332  c_ae = oneodx_e_v(ijk)
333  c_aw = oneodx_e_v(imjk)
334  c_an = oneody_n_v(ijk)
335  c_as = oneody_n_v(ijmk)
336  c_at = oneodz_t_v(ijk)
337  c_ab = oneodz_t_v(ijkm)
338  ELSE
339  c_ae = odx_e(i)
340  c_aw = odx_e(im)
341  c_an = ody(jp)
342  c_as = ody(j)
343  c_at = odz_t(k)
344  c_ab = odz_t(km)
345  ENDIF
346 
347 ! East face (i+1/2, j+1/2, k)
348  d_fe = avg_y_h(avg_x_h(epmu_gt(ijkc),epmu_gt(ijke),i),&
349  avg_x_h(epmu_gt(ijkn),epmu_gt(ijkne),i),j)*&
350  c_ae*ayz_v(ijk)
351 ! West face (i-1/2, j+1/2, k)
352  d_fw = avg_y_h(avg_x_h(epmu_gt(ijkw),epmu_gt(ijkc),im),&
353  avg_x_h(epmu_gt(ijknw),epmu_gt(ijkn),im),j)*&
354  c_aw*ayz_v(imjk)
355 
356 ! North face (i, j+1, k)
357  d_fn = epmu_gt(ijkn)*c_an*axz_v(ijk)
358 ! South face (i, j, k)
359  d_fs = epmu_gt(ijkc)*c_as*axz_v(ijmk)
360 
361  d_ft = zero
362  d_fb = zero
363  IF (do_k) THEN
364  ijkt = top_of(ijk)
365  ijktn = north_of(ijkt)
366  ijkb = bottom_of(ijk)
367  ijkbn = north_of(ijkb)
368 
369 ! Top face (i, j+1/2, k+1/2)
370  d_ft = avg_y_h(avg_z_h(epmu_gt(ijkc),epmu_gt(ijkt),k),&
371  avg_z_h(epmu_gt(ijkn),epmu_gt(ijktn),k),j)*&
372  ox(i)*c_at*axy_v(ijk)
373 ! Bottom face (i, j+1/2, k-1/2)
374  d_fb = avg_y_h(avg_z_h(epmu_gt(ijkb),epmu_gt(ijkc),km),&
375  avg_z_h(epmu_gt(ijkbn),epmu_gt(ijkn),km),j)*&
376  ox(i)*c_ab*axy_v(ijkm)
377  ENDIF ! end if (do_k)
378 
379  df_gv(ijk,east) = d_fe
380  df_gv(ijk,west) = d_fw
381  df_gv(ijk,north) = d_fn
382  df_gv(ijk,south) = d_fs
383  df_gv(ijk,top) = d_ft
384  df_gv(ijk,bottom) = d_fb
385 
386 ! if jackson, implement jackson style governing equations: multiply by
387 ! the void fraction otherwise multiply by 1
388  epga = avg_y(epg_jfac(ijkc), epg_jfac(ijkn), j)
389  d_fe = epga*d_fe
390  d_fw = epga*d_fw
391  d_fn = epga*d_fn
392  d_fs = epga*d_fs
393  d_ft = epga*d_ft
394  d_fb = epga*d_fb
395 
396  RETURN
397 
398  CONTAINS
399 
400  include 'fun_avg.inc'
401 
402  END SUBROUTINE get_vcell_gdiff_terms
403 
404 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
405 ! C
406 ! Subroutine: STORE_A_V_g0 C
407 ! Purpose: Determine convection diffusion terms for V_g momentum eqs. C
408 ! The off-diagonal coefficients calculated here must be positive. C
409 ! The center coefficient and the source vector are negative. See C
410 ! source_v_g. C
411 ! Implement FOUP discretization C
412 ! C
413 ! Author: M. Syamlal Date: 6-JUN-96 C
414 ! C
415 ! Revision Number: 1 C
416 ! Purpose: To incorporate Cartesian grid modifications C
417 ! Author: Jeff Dietiker Date: 01-Jul-09 C
418 ! C
419 ! C
420 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
421  SUBROUTINE store_a_v_g0(A_V_G)
422 
423 ! Modules
424 !---------------------------------------------------------------------//
425  USE compar, only: ijkstart3, ijkend3
426 
427  USE functions, only: flow_at_n
428  USE functions, only: ip_of, jp_of, kp_of
429  USE functions, only: im_of, jm_of, km_of
430 
431  USE geometry, only: do_k
432 
433  USE param
434  USE param1, only: zero
435 
436  IMPLICIT NONE
437 
438 ! Dummy arguments
439 !---------------------------------------------------------------------//
440 ! Septadiagonal matrix A_V_g
441  DOUBLE PRECISION, INTENT(INOUT) :: A_V_g(dimension_3, -3:3, 0:dimension_m)
442 
443 ! Local variables
444 !---------------------------------------------------------------------//
445 ! Indices
446  INTEGER :: IJK
447  INTEGER :: IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
448 ! Face mass flux
449  DOUBLE PRECISION :: flux_e, flux_w, flux_n, flux_s
450  DOUBLE PRECISION :: flux_t, flux_b
451 ! Diffusion parameter
452  DOUBLE PRECISION :: D_fe, d_fw, d_fn, d_fs, d_ft, d_fb
453 
454 !---------------------------------------------------------------------//
455 
456 !$omp parallel do default(none) &
457 !$omp private(IJK, IMJK, IJMK, IPJK, IJPK, IJKM, IJKP, &
458 !$omp D_fe, d_fw, d_fn, d_fs, d_ft, d_fb, &
459 !$omp flux_e, flux_w, flux_n, flux_s, flux_t, &
460 !$omp flux_b) &
461 !$omp shared(ijkstart3, ijkend3, do_k, a_v_g)
462  DO ijk = ijkstart3, ijkend3
463 
464  IF (flow_at_n(ijk)) THEN
465 
466 ! Calculate convection-diffusion fluxes through each of the faces
467  CALL get_vcell_gcflux_terms(flux_e, flux_w, flux_n, &
468  flux_s, flux_t, flux_b, ijk)
469  CALL get_vcell_gdiff_terms(d_fe, d_fw, d_fn, d_fs, &
470  d_ft, d_fb, ijk)
471 
472  ipjk = ip_of(ijk)
473  ijpk = jp_of(ijk)
474  imjk = im_of(ijk)
475  ijmk = jm_of(ijk)
476 
477 ! East face (i+1/2, j+1/2, k)
478  IF (flux_e >= zero) THEN
479  a_v_g(ijk,east,0) = d_fe
480  a_v_g(ipjk,west,0) = d_fe + flux_e
481  ELSE
482  a_v_g(ijk,east,0) = d_fe - flux_e
483  a_v_g(ipjk,west,0) = d_fe
484  ENDIF
485 ! West face (i-1/2, j+1/2, k)
486  IF (.NOT.flow_at_n(imjk)) THEN
487  IF (flux_w >= zero) THEN
488  a_v_g(ijk,west,0) = d_fw + flux_w
489  ELSE
490  a_v_g(ijk,west,0) = d_fw
491  ENDIF
492  ENDIF
493 
494 
495 ! North face (i, j+1, k)
496  IF (flux_n >= zero) THEN
497  a_v_g(ijk,north,0) = d_fn
498  a_v_g(ijpk,south,0) = d_fn + flux_n
499  ELSE
500  a_v_g(ijk,north,0) = d_fn - flux_n
501  a_v_g(ijpk,south,0) = d_fn
502  ENDIF
503 ! South face (i, j, k)
504  IF (.NOT.flow_at_n(ijmk)) THEN
505  IF (flux_s >= zero) THEN
506  a_v_g(ijk,south,0) = d_fs + flux_s
507  ELSE
508  a_v_g(ijk,south,0) = d_fs
509  ENDIF
510  ENDIF
511 
512 
513  IF (do_k) THEN
514  ijkp = kp_of(ijk)
515  ijkm = km_of(ijk)
516 
517 ! Top face (i, j+1/2, k+1/2)
518  IF (flux_t >= zero) THEN
519  a_v_g(ijk,top,0) = d_ft
520  a_v_g(ijkp,bottom,0) = d_ft + flux_t
521  ELSE
522  a_v_g(ijk,top,0) = d_ft - flux_t
523  a_v_g(ijkp,bottom,0) = d_ft
524  ENDIF
525 ! Bottom face (i, j+1/2, k-1/2)
526  IF (.NOT.flow_at_n(ijkm)) THEN
527  IF (flux_b >= zero) THEN
528  a_v_g(ijk,bottom,0) = d_fb + flux_b
529  ELSE
530  a_v_g(ijk,bottom,0) = d_fb
531  ENDIF
532  ENDIF
533  ENDIF ! end if (do_k)
534 
535  ENDIF ! end if flow_at_n
536  END DO ! end do ijk
537 !$omp end parallel do
538 
539  RETURN
540  END SUBROUTINE store_a_v_g0
541 
542 
543 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
544 ! C
545 ! Subroutine: STORE_A_V_GDC C
546 ! Purpose: Use deferred correction method to solve the v-momentum C
547 ! equation. This method combines first order upwind and a user C
548 ! specified higher order method C
549 ! C
550 ! Author: C. GUENTHER Date: 8-APR-99 C
551 ! C
552 ! Revision Number: 1 C
553 ! Purpose: To incorporate Cartesian grid modifications C
554 ! Author: Jeff Dietiker Date: 01-Jul-09 C
555 ! C
556 ! C
557 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
558  SUBROUTINE store_a_v_gdc(B_M)
559 
560 ! Modules
561 !---------------------------------------------------------------------//
562  USE compar, only: ijkstart3, ijkend3
563 
564  USE discretization, only: fpfoi_of
565 
566  USE fldvar, only: v_g
567 
568  USE function3, only: funijk3
569  USE functions, only: flow_at_n
570  USE functions, only: ip_of, jp_of, kp_of
571  USE functions, only: im_of, jm_of, km_of
572 
573  USE geometry, only: do_k
574 
575  USE indices, only: i_of, j_of, k_of
576 
577  USE param, only: dimension_3, dimension_4
578  USE param1, only: zero
579 
580  USE run, only: discretize, fpfoi
581  USE sendrecv3, only: send_recv3
582  USE xsi, only: calc_xsi
583  IMPLICIT NONE
584 
585 ! Dummy arguments
586 !---------------------------------------------------------------------//
587 ! Vector b_m
588  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3)
589 
590 ! Local variables
591 !---------------------------------------------------------------------//
592 ! Indices
593  INTEGER :: I, J, K, IJK
594  INTEGER :: IPJK, IMJK, IJMK, IJPK, IJKM, IJKP
595  INTEGER :: IJK4, IPPP, IPPP4, JPPP, JPPP4, KPPP, KPPP4
596  INTEGER :: IMMM, IMMM4, JMMM, JMMM4, KMMM, KMMM4
597 ! indication for shear
598  INTEGER :: incr
599 ! deferred corrction contribution form high order method
600  DOUBLE PRECISION :: MOM_HO
601 ! low order approximation
602  DOUBLE PRECISION :: MOM_LO
603 ! convection factor at each face
604  DOUBLE PRECISION :: flux_e, flux_w, flux_n, flux_s
605  DOUBLE PRECISION :: flux_t, flux_b
606 ! deferred correction contributions from each face
607  DOUBLE PRECISION :: EAST_DC
608  DOUBLE PRECISION :: WEST_DC
609  DOUBLE PRECISION :: NORTH_DC
610  DOUBLE PRECISION :: SOUTH_DC
611  DOUBLE PRECISION :: TOP_DC
612  DOUBLE PRECISION :: BOTTOM_DC
613 
614 ! temporary use of global arrays:
615 ! array1 (locally u) - the x directional velocity
616  DOUBLE PRECISION :: U(dimension_3)
617 ! array2 (locally v) - the y directional velocity
618  DOUBLE PRECISION :: V(dimension_3)
619 ! array3 (locally ww) - the z directional velocity
620  DOUBLE PRECISION :: WW(dimension_3)
621 !---------------------------------------------------------------------//
622  DOUBLE PRECISION :: TMP4(dimension_4)
623  DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: XSI_e, XSI_n, XSI_t
624 
625  CALL get_vcell_gvterms(u, v, ww)
626 
627 ! Send recv the third ghost layer
628  IF (fpfoi) THEN
629  Do ijk = ijkstart3, ijkend3
630  i = i_of(ijk)
631  j = j_of(ijk)
632  k = k_of(ijk)
633  ijk4 = funijk3(i,j,k)
634  tmp4(ijk4) = v_g(ijk)
635  ENDDO
636  CALL send_recv3(tmp4)
637  ENDIF
638 
639 ! shear indicator: v-momentum
640  incr=2
641  CALL calc_xsi (discretize(4), v_g, u, v, ww, xsi_e, xsi_n, &
642  xsi_t, incr)
643 
644 !!!$omp parallel do &
645 !!!$omp& private( I, J, K, IJK, IPJK, IMJK, IJPK, IJMK, &
646 !!!$omp& flux_e, flux_w, flux_n, flux_s, flux_b, flux_t, &
647 !!!$omp& MOM_HO, MOM_LO, EAST_DC, WEST_DC, NORTH_DC, &
648 !!!$omp& SOUTH_DC, TOP_DC, BOTTOM_DC)
649  DO ijk = ijkstart3, ijkend3
650 
651  IF (flow_at_n(ijk)) THEN
652 
653 ! Calculate convection fluxes through each of the faces
654  CALL get_vcell_gcflux_terms(flux_e, flux_w, flux_n, &
655  flux_s, flux_t, flux_b, ijk)
656 
657  ipjk = ip_of(ijk)
658  imjk = im_of(ijk)
659  ijpk = jp_of(ijk)
660  ijmk = jm_of(ijk)
661  ijkp = kp_of(ijk)
662  ijkm = km_of(ijk)
663 
664 ! Third Ghost layer information
665  ippp = ip_of(ip_of(ipjk))
666  ippp4 = funijk3(i_of(ippp), j_of(ippp), k_of(ippp))
667  immm = im_of(im_of(imjk))
668  immm4 = funijk3(i_of(immm), j_of(immm), k_of(immm))
669  jppp = jp_of(jp_of(ijpk))
670  jppp4 = funijk3(i_of(jppp), j_of(jppp), k_of(jppp))
671  jmmm = jm_of(jm_of(ijmk))
672  jmmm4 = funijk3(i_of(jmmm), j_of(jmmm), k_of(jmmm))
673  kppp = kp_of(kp_of(ijkp))
674  kppp4 = funijk3(i_of(kppp), j_of(kppp), k_of(kppp))
675  kmmm = km_of(km_of(ijkm))
676  kmmm4 = funijk3(i_of(kmmm), j_of(kmmm), k_of(kmmm))
677 
678 
679 ! East face (i+1/2, j+1/2, k)
680  IF(u(ijk) >= zero)THEN
681  mom_lo = v_g(ijk)
682  IF (fpfoi) mom_ho = fpfoi_of(v_g(ipjk), v_g(ijk), &
683  v_g(imjk), v_g(im_of(imjk)))
684  ELSE
685  mom_lo = v_g(ipjk)
686  IF (fpfoi) mom_ho = fpfoi_of(v_g(ijk), v_g(ipjk), &
687  v_g(ip_of(ipjk)), tmp4(ippp4))
688  ENDIF
689  IF (.NOT. fpfoi) mom_ho = xsi_e(ijk)*v_g(ipjk)+&
690  (1.0-xsi_e(ijk))*v_g(ijk)
691  east_dc = flux_e*(mom_lo-mom_ho)
692 
693 
694 ! North face (i, j+1, k)
695  IF(v(ijk) >= zero)THEN
696  mom_lo = v_g(ijk)
697  IF (fpfoi) mom_ho = fpfoi_of(v_g(ijpk), v_g(ijk), &
698  v_g(ijmk), v_g(jm_of(ijmk)))
699  ELSE
700  mom_lo = v_g(ijpk)
701  IF (fpfoi) mom_ho = fpfoi_of(v_g(ijk), v_g(ijpk), &
702  v_g(jp_of(ijpk)), tmp4(jppp4))
703  ENDIF
704  IF (.NOT. fpfoi) mom_ho = xsi_n(ijk)*v_g(ijpk)+&
705  (1.0-xsi_n(ijk))*v_g(ijk)
706  north_dc = flux_n*(mom_lo-mom_ho)
707 
708 
709 ! Top face (i, j+1/2, k+1/2)
710  IF (do_k) THEN
711  IF(ww(ijk) >= zero)THEN
712  mom_lo = v_g(ijk)
713  IF (fpfoi) mom_ho = fpfoi_of(v_g(ijkp), v_g(ijk), &
714  v_g(ijkm), v_g(km_of(ijkm)))
715  ELSE
716  mom_lo = v_g(ijkp)
717  IF (fpfoi) mom_ho = fpfoi_of(v_g(ijk), v_g(ijkp), &
718  v_g(kp_of(ijkp)), tmp4(kppp4))
719  ENDIF
720  IF (.NOT. fpfoi) mom_ho = xsi_t(ijk)*v_g(ijkp)+&
721  (1.0-xsi_t(ijk))*v_g(ijk)
722  top_dc = flux_t*(mom_lo-mom_ho)
723  ELSE
724  top_dc = zero
725  ENDIF
726 
727 
728 ! West face (i-1/2, j+1/2, k)
729  IF(u(imjk) >= zero)THEN
730  mom_lo = v_g(imjk)
731  IF (fpfoi) mom_ho = fpfoi_of(v_g(ijk), v_g(imjk), &
732  v_g(im_of(imjk)), tmp4(immm4))
733  ELSE
734  mom_lo = v_g(ijk)
735  IF (fpfoi) mom_ho = fpfoi_of(v_g(imjk), v_g(ijk), &
736  v_g(ipjk), v_g(ip_of(ipjk)))
737  ENDIF
738  IF (.NOT. fpfoi)mom_ho = xsi_e(imjk)*v_g(ijk)+&
739  (1.0-xsi_e(imjk))*v_g(imjk)
740  west_dc = flux_w*(mom_lo-mom_ho)
741 
742 
743 ! South face (i, j-1, k)
744  IF(v(ijmk) >= zero)THEN
745  mom_lo = v_g(ijmk)
746  IF (fpfoi) mom_ho = fpfoi_of(v_g(ijk), v_g(ijmk), &
747  v_g(jm_of(ijmk)), tmp4(jmmm4))
748  ELSE
749  mom_lo = v_g(ijk)
750  IF (fpfoi) mom_ho = fpfoi_of(v_g(ijmk), v_g(ijk), &
751  v_g(ijpk), v_g(jp_of(ijpk)))
752  ENDIF
753  IF (.NOT. fpfoi) mom_ho = xsi_n(ijmk)*v_g(ijk)+&
754  (1.0-xsi_n(ijmk))*v_g(ijmk)
755  south_dc = flux_s*(mom_lo-mom_ho)
756 
757 
758 ! Bottom face (i, j+1/2, k-1/2)
759  IF (do_k) THEN
760  IF(ww(ijk) >= zero)THEN
761  mom_lo = v_g(ijkm)
762  IF (fpfoi) mom_ho = fpfoi_of(v_g(ijk), v_g(ijkm), &
763  v_g(km_of(ijkm)), tmp4(kmmm4))
764  ELSE
765  mom_lo = v_g(ijk)
766  IF (fpfoi) mom_ho = fpfoi_of(v_g(ijkm), v_g(ijk), &
767  v_g(ijkp), v_g(kp_of(ijkp)))
768  ENDIF
769  IF (.NOT. fpfoi) mom_ho = xsi_t(ijkm)*v_g(ijk)+&
770  (1.0-xsi_t(ijkm))*v_g(ijkm)
771  bottom_dc = flux_b*(mom_lo-mom_ho)
772  ELSE
773  bottom_dc = zero
774  ENDIF ! end if (do_k)
775 
776 
777 ! CONTRIBUTION DUE TO DEFERRED CORRECTION
778  b_m(ijk) = b_m(ijk)+west_dc-east_dc+south_dc-north_dc&
779  +bottom_dc-top_dc
780 
781  ENDIF !end if flow_at_n
782  ENDDO ! end do ijk
783 
784  RETURN
785  END SUBROUTINE store_a_v_gdc
786 
787 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
788 ! C
789 ! Subroutine: STORE_A_V_g1 C
790 ! Purpose: Determine convection diffusion terms for V_g momentum eqs C
791 ! The off-diagonal coefficients calculated here must be positive. C
792 ! The center coefficient and the source vector are negative. C
793 ! Implements higher order discretization. C
794 ! See source_v_g C
795 ! C
796 ! Author: M. Syamlal Date: 20-MAR-97 C
797 ! C
798 ! Revision Number: 1 C
799 ! Purpose: To incorporate Cartesian grid modifications C
800 ! Author: Jeff Dietiker Date: 01-Jul-09 C
801 ! C
802 ! C
803 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
804  SUBROUTINE store_a_v_g1(A_V_G)
805 
806 ! Modules
807 !---------------------------------------------------------------------//
808  USE compar, only: ijkstart3, ijkend3
809  USE fldvar, only: v_g
810 
811  USE functions, only: flow_at_n
812  USE functions, only: ip_of, jp_of, kp_of
813  USE functions, only: im_of, jm_of, km_of
814  USE geometry, only: do_k
815  USE param
816  USE param1, only: one
817  USE run, only: discretize
818  USE xsi, only: calc_xsi
819  IMPLICIT NONE
820 
821 ! Dummy arguments
822 !---------------------------------------------------------------------//
823 ! Septadiagonal matrix A_V_g
824  DOUBLE PRECISION, INTENT(INOUT) :: A_V_g(dimension_3, -3:3)
825 
826 ! Local variables
827 !---------------------------------------------------------------------//
828 ! Indices
829  INTEGER :: IJK, IPJK, IMJK, IJPK, IJMK, IJKP, IJKM
830 ! indicator for shear
831  INTEGER :: incr
832 ! Diffusion parameter
833  DOUBLE PRECISION :: d_fe, d_fw, d_fn, d_fs, d_ft, d_fb
834 ! Face mass flux
835  DOUBLE PRECISION :: Flux_e, flux_w, flux_n, flux_s
836  DOUBLE PRECISION :: flux_t, flux_b
837 
838 ! temporary use of global arrays:
839 ! array1 (locally u) - the x directional velocity
840  DOUBLE PRECISION :: U(dimension_3)
841 ! array2 (locally v) - the y directional velocity
842  DOUBLE PRECISION :: V(dimension_3)
843 ! array3 (locally ww) - the z directional velocity
844  DOUBLE PRECISION :: WW(dimension_3)
845 !---------------------------------------------------------------------//
846  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: TMP4
847  DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: XSI_e, XSI_n, XSI_t
848 
849  CALL get_vcell_gvterms(u, v, ww)
850 
851 ! shear indicator: y-momentum
852  incr=2
853  CALL calc_xsi (discretize(4), v_g, u, v, ww, xsi_e, xsi_n, &
854  xsi_t, incr)
855 
856 !!!$omp parallel do &
857 !!!$omp& private(IJK, IMJK, IJMK, IPJK, IJPK, IJKM, IJKP, &
858 !!!$omp& d_fe, d_fw, d_fn, d_fs, d_ft, d_fb, &
859 !!!$omp& flux_e, flux_w, flux_n, flux_s, flux_t, flux_b)
860  DO ijk = ijkstart3, ijkend3
861 
862  IF (flow_at_n(ijk)) THEN
863 
864 ! Calculate convection-diffusion fluxes through each of the faces
865  CALL get_vcell_gcflux_terms(flux_e, flux_w, flux_n, &
866  flux_s, flux_t, flux_b, ijk)
867  CALL get_vcell_gdiff_terms(d_fe, d_fw, d_fn, d_fs, &
868  d_ft, d_fb, ijk)
869 
870  imjk = im_of(ijk)
871  ijmk = jm_of(ijk)
872  ipjk = ip_of(ijk)
873  ijpk = jp_of(ijk)
874 
875 ! East face (i+1/2, j+1/2, k)
876  a_v_g(ijk,east) = d_fe - xsi_e(ijk)*flux_e
877  a_v_g(ipjk,west) = d_fe + (one - xsi_e(ijk))*flux_e
878 ! West face (i-1/2, j+1/2, k)
879  IF (.NOT.flow_at_n(imjk)) THEN
880  a_v_g(ijk,west) = d_fw + (one - xsi_e(imjk))*flux_w
881  ENDIF
882 
883 
884 ! North face (i, j+1, k)
885  a_v_g(ijk,north) = d_fn - xsi_n(ijk)*flux_n
886  a_v_g(ijpk,south) = d_fn + (one - xsi_n(ijk))*flux_n
887 ! South face (i, j, k)
888  IF (.NOT.flow_at_n(ijmk)) THEN
889  a_v_g(ijk,south) = d_fs + (one - xsi_n(ijmk))*flux_s
890  ENDIF
891 
892 
893  IF (do_k) THEN
894  ijkp = kp_of(ijk)
895  ijkm = km_of(ijk)
896 ! Top face (i, j+1/2, k+1/2)
897  a_v_g(ijk,top) = d_ft - xsi_t(ijk)*flux_t
898  a_v_g(ijkp,bottom) = d_ft + (one - xsi_t(ijk))*flux_t
899 ! Bottom face (i, j+1/2, k-1/2)
900  IF (.NOT.flow_at_n(ijkm)) THEN
901  a_v_g(ijk,bottom) = d_fb + (one - xsi_t(ijkm))*flux_b
902  ENDIF
903  ENDIF ! end if (do_k)
904 
905  ENDIF ! end if (flow_at_n)
906  ENDDO ! end do (ijk)
907 
908  RETURN
909  END SUBROUTINE store_a_v_g1
mflux::flux_ge
double precision, dimension(:), allocatable flux_ge
Definition: mflux_mod.f:13
run::momentum_y_eq
logical, dimension(0:dim_m) momentum_y_eq
Definition: run_mod.f:77
get_vcell_gdiff_terms
subroutine get_vcell_gdiff_terms(D_FE, D_FW, D_FN, D_FS, D_FT, D_FB, IJK)
Definition: conv_dif_v_g.f:261
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
store_a_v_g0
subroutine store_a_v_g0(A_V_G)
Definition: conv_dif_v_g.f:422
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
cutcell::alpha_vt_c
double precision, dimension(:), allocatable alpha_vt_c
Definition: cutcell_mod.f:269
cutcell::oneodx_e_v
double precision, dimension(:), allocatable oneodx_e_v
Definition: cutcell_mod.f:320
geometry::ody
double precision, dimension(:), allocatable ody
Definition: geometry_mod.f:116
store_a_v_gdc
subroutine store_a_v_gdc(B_M)
Definition: conv_dif_v_g.f:559
param1::one
double precision, parameter one
Definition: param1_mod.f:29
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
store_a_v_g1
subroutine store_a_v_g1(A_V_G)
Definition: conv_dif_v_g.f:805
get_interpolation_terms_g
subroutine get_interpolation_terms_g(IJK, TYPE_OF_CELL, ALPHA_CUT, AW, HW, VELW)
Definition: get_cut_cell_volume_area.f:1448
conv_dif_v_g
subroutine conv_dif_v_g(A_M, B_M, IER)
Definition: conv_dif_v_g.f:14
indices::im1
integer, dimension(:), allocatable im1
Definition: indices_mod.f:50
cutcell::alpha_ve_c
double precision, dimension(:), allocatable alpha_ve_c
Definition: cutcell_mod.f:260
visc_g::df_gv
double precision, dimension(:,:), allocatable df_gv
Definition: visc_g_mod.f:32
fldvar::epg_jfac
double precision, dimension(:), allocatable epg_jfac
Definition: fldvar_mod.f:18
xsi::calc_xsi
subroutine calc_xsi(DISCR, PHI, U, V, W, XSI_E, XSI_N, XSI_T, incr)
Definition: xsi_mod.f:23
param::east
integer east
Definition: param_mod.f:29
xsi
Definition: xsi_mod.f:15
geometry::ayz_v
double precision, dimension(:), allocatable ayz_v
Definition: geometry_mod.f:227
get_vcell_gvterms
subroutine get_vcell_gvterms(U, V, WW)
Definition: conv_dif_v_g.f:63
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
discretization::fpfoi_of
double precision function fpfoi_of(PHI_D, PHI_C, PHI_U, PHI_UU)
Definition: discretization_mod.f:366
fun_avg
Definition: fun_avg_mod.f:1
cutcell::theta_vn_bar
double precision, dimension(:), allocatable theta_vn_bar
Definition: cutcell_mod.f:252
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
geometry::odx_e
double precision, dimension(:), allocatable odx_e
Definition: geometry_mod.f:121
get_vcell_gcflux_terms
subroutine get_vcell_gcflux_terms(FLUX_E, FLUX_W, FLUX_N, FLUX_S, FLUX_T, FLUX_B, IJK)
Definition: conv_dif_v_g.f:150
mflux::flux_gn
double precision, dimension(:), allocatable flux_gn
Definition: mflux_mod.f:15
cutcell::theta_v_ne
double precision, dimension(:), allocatable theta_v_ne
Definition: cutcell_mod.f:248
geometry::axy_v
double precision, dimension(:), allocatable axy_v
Definition: geometry_mod.f:231
geometry::ox
double precision, dimension(:), allocatable ox
Definition: geometry_mod.f:140
indices::jp1
integer, dimension(:), allocatable jp1
Definition: indices_mod.f:51
cutcell::oneody_n_v
double precision, dimension(:), allocatable oneody_n_v
Definition: cutcell_mod.f:321
param::dimension_4
integer dimension_4
Definition: param_mod.f:16
visc_g::epmu_gt
double precision, dimension(:), allocatable epmu_gt
Definition: visc_g_mod.f:13
cutcell::theta_vn
double precision, dimension(:), allocatable theta_vn
Definition: cutcell_mod.f:251
param::north
integer north
Definition: param_mod.f:37
fldvar::v_g
double precision, dimension(:), allocatable v_g
Definition: fldvar_mod.f:99
sendrecv3::send_recv3
Definition: sendrecv3_mod.f:75
param::south
integer south
Definition: param_mod.f:41
fldvar::w_g
double precision, dimension(:), allocatable w_g
Definition: fldvar_mod.f:111
param1::half
double precision, parameter half
Definition: param1_mod.f:28
run
Definition: run_mod.f:13
cutcell::theta_v_nt
double precision, dimension(:), allocatable theta_v_nt
Definition: cutcell_mod.f:254
param
Definition: param_mod.f:2
cutcell::theta_v_se
double precision, dimension(:), allocatable theta_v_se
Definition: cutcell_mod.f:249
run::def_cor
logical def_cor
Definition: run_mod.f:103
mflux::flux_gt
double precision, dimension(:), allocatable flux_gt
Definition: mflux_mod.f:17
dif_v_is
subroutine dif_v_is(DIF, A_M, M)
Definition: dif_v_is.f:11
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
indices::km1
integer, dimension(:), allocatable km1
Definition: indices_mod.f:52
cutcell::cut_v_treatment_at
logical, dimension(:), allocatable cut_v_treatment_at
Definition: cutcell_mod.f:351
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
visc_g
Definition: visc_g_mod.f:1
param::west
integer west
Definition: param_mod.f:33
fldvar::u_g
double precision, dimension(:), allocatable u_g
Definition: fldvar_mod.f:87
param::top
integer top
Definition: param_mod.f:45
run::discretize
integer, dimension(dim_eqs) discretize
Definition: run_mod.f:67
cutcell::oneodz_t_v
double precision, dimension(:), allocatable oneodz_t_v
Definition: cutcell_mod.f:322
cutcell::theta_v_st
double precision, dimension(:), allocatable theta_v_st
Definition: cutcell_mod.f:255
param::dimension_m
integer dimension_m
Definition: param_mod.f:18
geometry::odz_t
double precision, dimension(:), allocatable odz_t
Definition: geometry_mod.f:125
function3::funijk3
integer function funijk3(LI3, LJ3, LK3)
Definition: function3_mod.f:49
run::fpfoi
logical fpfoi
Definition: run_mod.f:106
param::bottom
integer bottom
Definition: param_mod.f:49
geometry::axz_v
double precision, dimension(:), allocatable axz_v
Definition: geometry_mod.f:229
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
cutcell::alpha_vn_c
double precision, dimension(:), allocatable alpha_vn_c
Definition: cutcell_mod.f:263