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