MFIX  2016-1
xsi_mod.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Module name: XSI C
4 ! Purpose: Determine convection weighting factors for higher order C
5 ! discretization. C
6 ! C
7 ! Author: M. Syamlal Date: 6-MAR-97 C
8 ! C
9 ! Comments: contains following functions subroutines & functions: C
10 ! calc_xsi, cxs, dw, C
11 ! xsi_func C
12 ! C
13 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
14 
15  MODULE xsi
16  IMPLICIT NONE
17 
18  CONTAINS
19 
20 
21  SUBROUTINE calc_xsi(DISCR, PHI, U, V, W, XSI_E, XSI_N, XSI_T, &
22  incr)
23 
24 ! Modules
25 !---------------------------------------------------------------------//
26  USE chischeme, only: chi_flag
27  USE chischeme, only: chi_e, chi_n, chi_t
28 
29  USE compar, only: ijkstart3, ijkend3
30 
31  USE discretization, only: phi_c_of
32  USE discretization, only: superbee
33  USE discretization, only: chi_smart, smart
34  USE discretization, only: ultra_quick
35  USE discretization, only: quickest
36  USE discretization, only: chi_muscl, muscl
37  USE discretization, only: vanleer
38  USE discretization, only: minmod
39  USE discretization, only: central_scheme
40 
41  USE functions, only: east_of, west_of, north_of, south_of
42  USE functions, only: top_of, bottom_of
43 
44  USE geometry, only: do_k
45  USE geometry, only: odx, odx_e, ody, ody_n, odz, odz_t
46  USE geometry, only: ox
47 
48  USE indices, only: i_of, j_of, k_of
49  USE indices, only: im1, ip1, jm1, jp1, km1, kp1
50 
51  USE param, only: dimension_3
52 
53  USE param1, only: zero
54 
55  USE run, only: shear, dt
56 
57  USE sendrecv, only: send_recv
58 
59  USE error_manager, only: err_msg, init_err_msg, finl_err_msg
60  USE error_manager, only: ival, flush_err_msg
61  IMPLICIT NONE
62 
63 ! Dummy arguments
64 !---------------------------------------------------------------------//
65 ! discretization method
66  INTEGER, INTENT(IN) :: DISCR
67 ! convected quantity
68  DOUBLE PRECISION, INTENT(IN) :: PHI(dimension_3)
69 ! Velocity components
70  DOUBLE PRECISION, INTENT(IN) :: U(dimension_3)
71  DOUBLE PRECISION, INTENT(IN) :: V(dimension_3)
72  DOUBLE PRECISION, INTENT(IN) :: W(dimension_3)
73 ! Convection weighting factors
74  DOUBLE PRECISION, INTENT(OUT) :: XSI_e(dimension_3)
75  DOUBLE PRECISION, INTENT(OUT) :: XSI_n(dimension_3)
76  DOUBLE PRECISION, INTENT(OUT) :: XSI_t(dimension_3)
77 ! shear indicator
78  INTEGER, INTENT(IN) :: incr
79 
80 ! Local variables
81 !---------------------------------------------------------------------//
82 ! Indices
83  INTEGER :: IJK, IJKC, IJKD, IJKU, I, J, K
84 !
85  DOUBLE PRECISION :: PHI_C
86 ! down wind factor
87  DOUBLE PRECISION :: dwf
88 ! Courant number
89  DOUBLE PRECISION :: cf
90 ! cell widths for QUICKEST
91  DOUBLE PRECISION :: oDXc, oDXuc, oDYc, oDYuc, oDZc, oDZuc
92 !---------------------------------------------------------------------//
93 
94  IF (shear) THEN
95 ! calculate XSI_E, XSI_N, XSI_T when periodic shear BCs are used
96  call cxs(incr, discr, u, v, w, phi, xsi_e, xsi_n, xsi_t)
97  ELSE
98 
99  SELECT CASE (discr) !first order upwinding
100  CASE (:1)
101 
102 !$omp parallel do default(none) private(IJK) shared(ijkstart3, ijkend3, u, v, w, xsi_e, xsi_n, xsi_t, do_k)
103  DO ijk = ijkstart3, ijkend3
104  xsi_e(ijk) = xsi_func(u(ijk),zero)
105  xsi_n(ijk) = xsi_func(v(ijk),zero)
106  IF (do_k) xsi_t(ijk) = xsi_func(w(ijk),zero)
107  ENDDO
108 
109 
110  CASE (2) !Superbee
111 
112 !!!$omp parallel do private(IJK, IJKC,IJKD,IJKU, PHI_C,DWF)
113  DO ijk = ijkstart3, ijkend3
114 
115  IF (u(ijk) >= zero) THEN
116  ijkc = ijk
117  ijkd = east_of(ijk)
118  ijku = west_of(ijkc)
119  ELSE
120  ijkc = east_of(ijk)
121  ijkd = ijk
122  ijku = east_of(ijkc)
123  ENDIF
124  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
125  dwf = superbee(phi_c)
126  xsi_e(ijk) = xsi_func(u(ijk),dwf)
127 
128  IF (v(ijk) >= zero) THEN
129  ijkc = ijk
130  ijkd = north_of(ijk)
131  ijku = south_of(ijkc)
132  ELSE
133  ijkc = north_of(ijk)
134  ijkd = ijk
135  ijku = north_of(ijkc)
136  ENDIF
137  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
138  dwf = superbee(phi_c)
139  xsi_n(ijk) = xsi_func(v(ijk),dwf)
140 
141  IF (do_k) THEN
142  IF (w(ijk) >= zero) THEN
143  ijkc = ijk
144  ijkd = top_of(ijk)
145  ijku = bottom_of(ijkc)
146  ELSE
147  ijkc = top_of(ijk)
148  ijkd = ijk
149  ijku = top_of(ijkc)
150  ENDIF
151  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
152  dwf = superbee(phi_c)
153  xsi_t(ijk) = xsi_func(w(ijk),dwf)
154  ENDIF
155  ENDDO
156 
157 
158  CASE (3) !SMART
159 
160 !!!$omp parallel do private(IJK, IJKC,IJKD,IJKU, PHI_C,DWF)
161  DO ijk = ijkstart3, ijkend3
162 
163  IF (u(ijk) >= zero) THEN
164  ijkc = ijk
165  ijkd = east_of(ijk)
166  ijku = west_of(ijkc)
167  ELSE
168  ijkc = east_of(ijk)
169  ijkd = ijk
170  ijku = east_of(ijkc)
171  ENDIF
172  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
173  if(chi_flag) then
174  dwf = chi_smart(phi_c, chi_e(ijk))
175  else
176  dwf = smart(phi_c)
177  endif
178  xsi_e(ijk) = xsi_func(u(ijk),dwf)
179 
180  IF (v(ijk) >= zero) THEN
181  ijkc = ijk
182  ijkd = north_of(ijk)
183  ijku = south_of(ijkc)
184  ELSE
185  ijkc = north_of(ijk)
186  ijkd = ijk
187  ijku = north_of(ijkc)
188  ENDIF
189  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
190  if(chi_flag) then
191  dwf = chi_smart(phi_c, chi_n(ijk))
192  else
193  dwf = smart(phi_c)
194  endif
195  xsi_n(ijk) = xsi_func(v(ijk),dwf)
196 
197  IF (do_k) THEN
198  IF (w(ijk) >= zero) THEN
199  ijkc = ijk
200  ijkd = top_of(ijk)
201  ijku = bottom_of(ijkc)
202  ELSE
203  ijkc = top_of(ijk)
204  ijkd = ijk
205  ijku = top_of(ijkc)
206  ENDIF
207  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
208  if(chi_flag) then
209  dwf = chi_smart(phi_c, chi_t(ijk))
210  else
211  dwf = smart(phi_c)
212  endif
213  xsi_t(ijk) = xsi_func(w(ijk),dwf)
214  ENDIF
215  ENDDO
216 
217 
218  CASE (4) !ULTRA-QUICK
219 
220 !!!$omp parallel do private(IJK, I,J,K, IJKC,IJKD,IJKU, PHI_C,DWF,CF)
221  DO ijk = ijkstart3, ijkend3
222 
223  i = i_of(ijk)
224  IF (u(ijk) >= zero) THEN
225  ijkc = ijk
226  ijkd = east_of(ijk)
227  ijku = west_of(ijkc)
228  ELSE
229  ijkc = east_of(ijk)
230  ijkd = ijk
231  ijku = east_of(ijkc)
232  ENDIF
233  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
234  cf = abs(u(ijk))*dt*odx_e(i)
235  dwf = ultra_quick(phi_c,cf)
236  xsi_e(ijk) = xsi_func(u(ijk),dwf)
237 
238  j = j_of(ijk)
239  IF (v(ijk) >= zero) THEN
240  ijkc = ijk
241  ijkd = north_of(ijk)
242  ijku = south_of(ijkc)
243  ELSE
244  ijkc = north_of(ijk)
245  ijkd = ijk
246  ijku = north_of(ijkc)
247  ENDIF
248  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
249  cf = abs(v(ijk))*dt*ody_n(j)
250  dwf = ultra_quick(phi_c,cf)
251  xsi_n(ijk) = xsi_func(v(ijk),dwf)
252 
253  IF (do_k) THEN
254  k = k_of(ijk)
255  IF (w(ijk) >= zero) THEN
256  ijkc = ijk
257  ijkd = top_of(ijk)
258  ijku = bottom_of(ijkc)
259  ELSE
260  ijkc = top_of(ijk)
261  ijkd = ijk
262  ijku = top_of(ijkc)
263  ENDIF
264  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
265  cf = abs(w(ijk))*dt*ox(i)*odz_t(k)
266  dwf = ultra_quick(phi_c,cf)
267  xsi_t(ijk) = xsi_func(w(ijk),dwf)
268  ENDIF
269  ENDDO
270 
271 
272  CASE (5) !QUICKEST
273 
274 !!!$omp parallel do &
275 !!!$omp& private(IJK,I,J,K, IJKC,IJKD,IJKU, &
276 !!!$omp& ODXC,ODXUC, PHI_C,CF,DWF, &
277 !!!$omp& ODYC,ODYUC, ODZC,ODZUC )
278  DO ijk = ijkstart3, ijkend3
279 
280  i = i_of(ijk)
281  IF (u(ijk) >= zero) THEN
282  ijkc = ijk
283  ijkd = east_of(ijk)
284  ijku = west_of(ijkc)
285  odxc = odx(i)
286  odxuc = odx_e(im1(i))
287  ELSE
288  ijkc = east_of(ijk)
289  ijkd = ijk
290  ijku = east_of(ijkc)
291  odxc = odx(ip1(i))
292  odxuc = odx_e(ip1(i))
293  ENDIF
294  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
295  cf = abs(u(ijk))*dt*odx_e(i)
296  dwf = quickest(phi_c,cf,odxc,odxuc,odx_e(i))
297  xsi_e(ijk) = xsi_func(u(ijk),dwf)
298 
299  j = j_of(ijk)
300  IF (v(ijk) >= zero) THEN
301  ijkc = ijk
302  ijkd = north_of(ijk)
303  ijku = south_of(ijkc)
304  odyc = ody(j)
305  odyuc = ody_n(jm1(j))
306  ELSE
307  ijkc = north_of(ijk)
308  ijkd = ijk
309  ijku = north_of(ijkc)
310  odyc = ody(jp1(j))
311  odyuc = ody_n(jp1(j))
312  ENDIF
313  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
314  cf = abs(v(ijk))*dt*ody_n(j)
315  dwf = quickest(phi_c,cf,odyc,odyuc,ody_n(j))
316  xsi_n(ijk) = xsi_func(v(ijk),dwf)
317 
318  IF (do_k) THEN
319  k = k_of(ijk)
320  IF (w(ijk) >= zero) THEN
321  ijkc = ijk
322  ijkd = top_of(ijk)
323  ijku = bottom_of(ijkc)
324  odzc = odz(k)
325  odzuc = odz_t(km1(k))
326  ELSE
327  ijkc = top_of(ijk)
328  ijkd = ijk
329  ijku = top_of(ijkc)
330  odzc = odz(kp1(k))
331  odzuc = odz_t(kp1(k))
332  ENDIF
333  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
334  cf = abs(w(ijk))*dt*ox(i)*odz_t(k)
335  dwf = quickest(phi_c,cf,odzc,odzuc,odz_t(k))
336  xsi_t(ijk) = xsi_func(w(ijk),dwf)
337  ENDIF
338  ENDDO
339 
340 
341  CASE (6) !MUSCL
342 
343 !!!$omp parallel do private(IJK, IJKC,IJKD,IJKU, PHI_C,DWF )
344  DO ijk = ijkstart3, ijkend3
345 
346  IF (u(ijk) >= zero) THEN
347  ijkc = ijk
348  ijkd = east_of(ijk)
349  ijku = west_of(ijkc)
350  ELSE
351  ijkc = east_of(ijk)
352  ijkd = ijk
353  ijku = east_of(ijkc)
354  ENDIF
355  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
356  if(chi_flag) then
357  dwf = chi_muscl(phi_c, chi_e(ijk))
358  else
359  dwf = muscl(phi_c)
360  endif
361  xsi_e(ijk) = xsi_func(u(ijk),dwf)
362 
363  IF (v(ijk) >= zero) THEN
364  ijkc = ijk
365  ijkd = north_of(ijk)
366  ijku = south_of(ijkc)
367  ELSE
368  ijkc = north_of(ijk)
369  ijkd = ijk
370  ijku = north_of(ijkc)
371  ENDIF
372  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
373  if(chi_flag) then
374  dwf = chi_muscl(phi_c, chi_n(ijk))
375  else
376  dwf = muscl(phi_c)
377  endif
378  xsi_n(ijk) = xsi_func(v(ijk),dwf)
379 
380  IF (do_k) THEN
381  IF (w(ijk) >= zero) THEN
382  ijkc = ijk
383  ijkd = top_of(ijk)
384  ijku = bottom_of(ijkc)
385  ELSE
386  ijkc = top_of(ijk)
387  ijkd = ijk
388  ijku = top_of(ijkc)
389  ENDIF
390  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
391  if(chi_flag) then
392  dwf = chi_muscl(phi_c, chi_t(ijk))
393  else
394  dwf = muscl(phi_c)
395  endif
396  xsi_t(ijk) = xsi_func(w(ijk),dwf)
397  ENDIF
398  ENDDO
399 
400 
401  CASE (7) !Van Leer
402 
403 !!!$omp parallel do private( IJK, IJKC,IJKD,IJKU, PHI_C,DWF )
404  DO ijk = ijkstart3, ijkend3
405 
406  IF (u(ijk) >= zero) THEN
407  ijkc = ijk
408  ijkd = east_of(ijk)
409  ijku = west_of(ijkc)
410  ELSE
411  ijkc = east_of(ijk)
412  ijkd = ijk
413  ijku = east_of(ijkc)
414  ENDIF
415  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
416  dwf = vanleer(phi_c)
417  xsi_e(ijk) = xsi_func(u(ijk),dwf)
418 
419  IF (v(ijk) >= zero) THEN
420  ijkc = ijk
421  ijkd = north_of(ijk)
422  ijku = south_of(ijkc)
423  ELSE
424  ijkc = north_of(ijk)
425  ijkd = ijk
426  ijku = north_of(ijkc)
427  ENDIF
428  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
429  dwf = vanleer(phi_c)
430  xsi_n(ijk) = xsi_func(v(ijk),dwf)
431 
432  IF (do_k) THEN
433  IF (w(ijk) >= zero) THEN
434  ijkc = ijk
435  ijkd = top_of(ijk)
436  ijku = bottom_of(ijkc)
437  ELSE
438  ijkc = top_of(ijk)
439  ijkd = ijk
440  ijku = top_of(ijkc)
441  ENDIF
442  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
443  dwf = vanleer(phi_c)
444  xsi_t(ijk) = xsi_func(w(ijk),dwf)
445  ENDIF
446  ENDDO
447 
448 
449  CASE (8) !Minmod
450 
451 !!!$omp parallel do private(IJK, IJKC,IJKD,IJKU, PHI_C,DWF )
452  DO ijk = ijkstart3, ijkend3
453 
454  IF (u(ijk) >= zero) THEN
455  ijkc = ijk
456  ijkd = east_of(ijk)
457  ijku = west_of(ijkc)
458  ELSE
459  ijkc = east_of(ijk)
460  ijkd = ijk
461  ijku = east_of(ijkc)
462  ENDIF
463  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
464  dwf = minmod(phi_c)
465  xsi_e(ijk) = xsi_func(u(ijk),dwf)
466 
467  IF (v(ijk) >= zero) THEN
468  ijkc = ijk
469  ijkd = north_of(ijk)
470  ijku = south_of(ijkc)
471  ELSE
472  ijkc = north_of(ijk)
473  ijkd = ijk
474  ijku = north_of(ijkc)
475  ENDIF
476  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
477  dwf = minmod(phi_c)
478  xsi_n(ijk) = xsi_func(v(ijk),dwf)
479 
480  IF (do_k) THEN
481  IF (w(ijk) >= zero) THEN
482  ijkc = ijk
483  ijkd = top_of(ijk)
484  ijku = bottom_of(ijkc)
485  ELSE
486  ijkc = top_of(ijk)
487  ijkd = ijk
488  ijku = top_of(ijkc)
489  ENDIF
490 
491  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
492  dwf = minmod(phi_c)
493  xsi_t(ijk) = xsi_func(w(ijk),dwf)
494  ENDIF
495  ENDDO
496 
497 
498  CASE (9) ! Central
499 
500 !!!$omp parallel do private(IJK, IJKC,IJKD,IJKU, PHI_C,DWF)
501  DO ijk = ijkstart3, ijkend3
502 
503  IF (u(ijk) >= zero) THEN
504  ijkc = ijk
505  ijkd = east_of(ijk)
506  ijku = west_of(ijkc)
507  ELSE
508  ijkc = east_of(ijk)
509  ijkd = ijk
510  ijku = east_of(ijkc)
511  ENDIF
512  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
513  dwf = central_scheme(phi_c)
514  xsi_e(ijk) = xsi_func(u(ijk),dwf)
515 
516  IF (v(ijk) >= zero) THEN
517  ijkc = ijk
518  ijkd = north_of(ijk)
519  ijku = south_of(ijkc)
520  ELSE
521  ijkc = north_of(ijk)
522  ijkd = ijk
523  ijku = north_of(ijkc)
524  ENDIF
525  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
526  dwf = central_scheme(phi_c)
527  xsi_n(ijk) = xsi_func(v(ijk),dwf)
528 
529  IF (do_k) THEN
530  IF (w(ijk) >= zero) THEN
531  ijkc = ijk
532  ijkd = top_of(ijk)
533  ijku = bottom_of(ijkc)
534  ELSE
535  ijkc = top_of(ijk)
536  ijkd = ijk
537  ijku = top_of(ijkc)
538  ENDIF
539  phi_c = phi_c_of(phi(ijku),phi(ijkc),phi(ijkd))
540  dwf = central_scheme(phi_c)
541  xsi_t(ijk) = xsi_func(w(ijk),dwf)
542  ENDIF
543  ENDDO
544 
545 
546  CASE DEFAULT !Error
547 ! should never hit this
548  CALL init_err_msg("CALC_XSI")
549  WRITE(err_msg, 1100) ival(discr)
550  1100 FORMAT('ERROR 1100: Invalid DISCRETIZE= ', a,' The check_data ',&
551  'routines should',/, 'have already caught this error and ',&
552  'prevented the simulation from ',/,'running. Please notify ',&
553  'the MFIX developers.')
554  CALL flush_err_msg(abort=.true.)
555  CALL finl_err_msg
556 
557  END SELECT
558 
559  ENDIF ! end if/else shear
560 
561  call send_recv(xsi_e,2)
562  call send_recv(xsi_n,2)
563  call send_recv(xsi_t,2)
564 
565  RETURN
566  END SUBROUTINE calc_xsi
567 
568 
569 
570 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
571 ! C
572 ! Subroutine: CXS C
573 ! Purpose: Calculates XSI_E, XSI_N,XSI_T when using periodic shear C
574 ! BC's. Uses true velocities. C
575 ! C
576 ! C
577 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
578  SUBROUTINE cxs(INCR, DISCR, U, V, WW, PHI, XSI_E, XSI_N, XSI_T)
579 
580 ! Modules
581 !---------------------------------------------------------------------//
582  USE param
583  USE param1
584  USE run
585  USE geometry
586  USE indices
587  USE vshear
588  USE compar
589  USE sendrecv
590  USE functions
591  IMPLICIT NONE
592 
593 ! Dummy argments
594 !---------------------------------------------------------------------//
595  INTEGER, INTENT(IN) :: incr
596  INTEGER, INTENT(IN) :: DISCR
597  DOUBLE PRECISION, INTENT(IN) :: V(dimension_3)
598  DOUBLE PRECISION, INTENT(IN) :: U(dimension_3)
599  DOUBLE PRECISION, INTENT(IN) :: WW(dimension_3)
600  DOUBLE PRECISION, INTENT(IN) :: PHI(dimension_3)
601  DOUBLE PRECISION, INTENT(OUT) :: XSI_E(dimension_3)
602  DOUBLE PRECISION, INTENT(OUT) :: XSI_N(dimension_3)
603  DOUBLE PRECISION, INTENT(OUT) :: XSI_T(dimension_3)
604 
605 ! Local variables
606 !---------------------------------------------------------------------//
607  DOUBLE PRECISION :: VV(dimension_3)
608  DOUBLE PRECISION :: SRT
609  DOUBLE PRECISION :: DWFE,DWFN,DWFT
610  DOUBLE PRECISION :: PHICU, PHIDU, PHIUU
611  DOUBLE PRECISION :: PHICV, PHIDV, PHIUV
612  DOUBLE PRECISION :: PHICW, PHIDW, PHIUW
613  INTEGER :: IJK, IJKC, IJKD, IJKU, I
614 !---------------------------------------------------------------------//
615 
616  IF (incr .eq. 2) THEN !V momentum balance
617  srt=(2d0*v_sh/xlength)
618 
619 !!!$omp parallel do private(IJK, IJKC,IJKD,IJKU,DWFE,DWFN,DWFT,&
620 !!!$omp& PHICU,PHIDU,PHIUU,PHICV,PHIDV,PHIUV,PHICW,PHIDW,PHIUW,I)&
621 !!!$omp& shared(DISCR)
622  DO ijk = ijkstart3, ijkend3
623 
624  i=i_of(ijk)
625  vv(ijk)=v(ijk)+vsh(ijk)
626 
627  IF (u(ijk) >= zero) THEN
628  ijkc = ijk
629  ijkd = east_of(ijk)
630  ijku = west_of(ijkc)
631  phicu=phi(ijkc)
632  phidu=phi(ijkd)+srt*1d0/odx_e(i)
633  phiuu=phi(ijku)-srt*1d0/odx_e(im1(i))
634  ELSE
635  ijkc = east_of(ijk)
636  ijkd = ijk
637  ijku = east_of(ijkc)
638  phicu=phi(ijkc)+srt*1d0/odx_e(i)
639  phidu=phi(ijkd)
640  phiuu=phi(ijku)+srt*1d0/odx_e(i)&
641  +srt*1d0/odx_e(ip1(i))
642  ENDIF
643 
644  IF (vv(ijk) >= zero) THEN
645  ijkc = ijk
646  ijkd = north_of(ijk)
647  ijku= south_of(ijkc)
648  phicv=phi(ijkc)
649  phidv=phi(ijkd)
650  phiuv=phi(ijku)
651  ELSE
652  ijkc= north_of(ijk)
653  ijkd= ijk
654  ijku= north_of(ijkc)
655  phicv=phi(ijkc)
656  phidv=phi(ijkd)
657  phiuv=phi(ijku)
658  ENDIF
659 
660  IF (do_k) THEN
661  IF (ww(ijk) >= zero) THEN
662  ijkc = ijk
663  ijkd = top_of(ijk)
664  ijku = bottom_of(ijkc)
665  ELSE
666  ijkc = top_of(ijk)
667  ijkd = ijk
668  ijku = top_of(ijkc)
669  ENDIF
670  phicw=phi(ijkc)
671  phidw=phi(ijkd)
672  phiuw=phi(ijku)
673  ELSE
674  phicw=0d0
675  phidw=0d0
676  phiuw=0d0
677  dwft=0d0
678  ENDIF
679 
680  CALL dw(u(ijk), vv(ijk), ww(ijk), ijk, phicu, phidu, phiuu, &
681  phicv, phidv, phiuv, phicw, phidw, phiuw, &
682  dwfe, dwfn, dwft, discr)
683 
684  xsi_e(ijk) = xsi_func(u(ijk),dwfe)
685  xsi_n(ijk) = xsi_func(vv(ijk),dwfn)
686  xsi_t(ijk) = xsi_func(ww(ijk),dwft)
687 
688  vv(ijk)=v(ijk)-vsh(ijk)
689  ENDDO
690 
691 
692  ELSEIF (incr .eq. 1) THEN !u momentum balance
693 
694 !!!$omp parallel do private(IJK, IJKC,IJKD,IJKU,DWFE,DWFN,DWFT,&
695 !!!$omp& PHICU,PHIDU,PHIUU,PHICV,PHIDV,PHIUV,PHICW,PHIDW,PHIUW,I)&
696 !!!$omp& shared(DISCR)
697  DO ijk = ijkstart3, ijkend3
698 
699  vv(ijk)=v(ijk)+vshe(ijk)
700  IF (u(ijk) >= zero) THEN
701  ijkc = ijk
702  ijkd = east_of(ijk)
703  ijku = west_of(ijkc)
704  phicu=phi(ijkc)
705  phidu=phi(ijkd)
706  phiuu=phi(ijku)
707  ELSE
708  ijkc = east_of(ijk)
709  ijkd = ijk
710  ijku = east_of(ijkc)
711  phicu=phi(ijkc)
712  phidu=phi(ijkd)
713  phiuu=phi(ijku)
714  ENDIF
715 
716  IF (vv(ijk) >= zero) THEN
717  ijkc = ijk
718  ijkd = north_of(ijk)
719  ijku= south_of(ijkc)
720  phicv=phi(ijkc)
721  phidv=phi(ijkd)
722  phiuv=phi(ijku)
723  ELSE
724  ijkc= north_of(ijk)
725  ijkd= ijk
726  ijku= north_of(ijkc)
727  phicv=phi(ijkc)
728  phidv=phi(ijkd)
729  phiuv=phi(ijku)
730  ENDIF
731 
732  IF (do_k) THEN
733  IF (ww(ijk) >= zero) THEN
734  ijkc = ijk
735  ijkd = top_of(ijk)
736  ijku = bottom_of(ijkc)
737  ELSE
738  ijkc = top_of(ijk)
739  ijkd = ijk
740  ijku = top_of(ijkc)
741  ENDIF
742  phicw=phi(ijkc)
743  phidw=phi(ijkd)
744  phiuw=phi(ijku)
745  ELSE
746  phicw=0d0
747  phidw=0d0
748  phiuw=0d0
749  dwft=0d0
750  ENDIF
751 
752  CALL dw(u(ijk), vv(ijk), ww(ijk), ijk, phicu, phidu, phiuu, &
753  phicv, phidv, phiuv, phicw, phidw, phiuw, &
754  dwfe, dwfn, dwft, discr)
755 
756  xsi_e(ijk) = xsi_func(u(ijk),dwfe)
757  xsi_n(ijk) = xsi_func(vv(ijk),dwfn)
758  xsi_t(ijk) = xsi_func(ww(ijk),dwft)
759 
760  vv(ijk)=v(ijk)-vshe(ijk)
761  ENDDO
762 
763 
764  ELSEIF (incr .eq. 0) THEN !scalars and w momentum
765 
766 !!!$omp parallel do private(IJK, IJKC,IJKD,IJKU,DWFE,DWFN,DWFT,&
767 !!!$omp& PHICU,PHIDU,PHIUU,PHICV,PHIDV,PHIUV,PHICW,PHIDW,PHIUW,I)&
768 !!!$omp& shared(DISCR)
769  DO ijk = ijkstart3, ijkend3
770 
771  vv(ijk)=v(ijk)+vsh(ijk)
772  IF (u(ijk) >= zero) THEN
773  ijkc = ijk
774  ijkd = east_of(ijk)
775  ijku = west_of(ijkc)
776  phicu=phi(ijkc)
777  phidu=phi(ijkd)
778  phiuu=phi(ijku)
779  ELSE
780  ijkc = east_of(ijk)
781  ijkd = ijk
782  ijku = east_of(ijkc)
783  phicu=phi(ijkc)
784  phidu=phi(ijkd)
785  phiuu=phi(ijku)
786  ENDIF
787 
788  IF (vv(ijk) >= zero) THEN
789  ijkc = ijk
790  ijkd = north_of(ijk)
791  ijku= south_of(ijkc)
792  phicv=phi(ijkc)
793  phidv=phi(ijkd)
794  phiuv=phi(ijku)
795  ELSE
796  ijkc= north_of(ijk)
797  ijkd= ijk
798  ijku= north_of(ijkc)
799  phicv=phi(ijkc)
800  phidv=phi(ijkd)
801  phiuv=phi(ijku)
802  ENDIF
803 
804  IF (do_k) THEN
805  IF (ww(ijk) >= zero) THEN
806  ijkc = ijk
807  ijkd = top_of(ijk)
808  ijku = bottom_of(ijkc)
809  ELSE
810  ijkc = top_of(ijk)
811  ijkd = ijk
812  ijku = top_of(ijkc)
813  ENDIF
814  phicw=phi(ijkc)
815  phidw=phi(ijkd)
816  phiuw=phi(ijku)
817  ELSE
818  phicw=0d0
819  phidw=0d0
820  phiuw=0d0
821  dwft=0d0
822  ENDIF
823 
824  CALL dw(u(ijk), vv(ijk), ww(ijk), ijk, phicu, phidu, phiuu, &
825  phicv, phidv, phiuv, phicw, phidw, phiuw, &
826  dwfe, dwfn, dwft, discr)
827 
828  xsi_e(ijk) = xsi_func(u(ijk),dwfe)
829  xsi_n(ijk) = xsi_func(vv(ijk),dwfn)
830  xsi_t(ijk) = xsi_func(ww(ijk),dwft)
831 
832  vv(ijk)=v(ijk)-vsh(ijk)
833  ENDDO
834 
835  ELSE
836  write(*,*) 'INCR ERROR'
837  ENDIF
838 
839  call send_recv(xsi_e,2)
840  call send_recv(xsi_n,2)
841  call send_recv(xsi_t,2)
842  RETURN
843  END SUBROUTINE cxs
844 
845 
846 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
847 ! C
848 ! Subroutine: DW C
849 ! Purpose: Calculates DWFs for various discretization methods when C
850 ! period shear BCs are used C
851 ! C
852 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
853  SUBROUTINE dw(UU, VV, WW, IJK, PHICU, PHIDU, PHIUU, &
854  phicv, phidv, phiuv, &
855  phicw, phidw, phiuw, &
856  dwfe, dwfn, dwft, discr)
857 
858 ! Modules
859 !---------------------------------------------------------------------//
860  USE compar
861  USE discretization
862  USE exit, only: mfix_exit
863  USE functions
864  USE geometry
865  USE indices
866  USE param
867  USE param1
868  USE run
869  USE sendrecv
870  USE vshear
871  IMPLICIT NONE
872 
873 ! Dummy argments
874 !---------------------------------------------------------------------//
875 !
876  DOUBLE PRECISION, INTENT(IN) :: UU, VV, WW
877 !
878  INTEGER, INTENT(IN) :: IJK
879 !
880  DOUBLE PRECISION, INTENT(IN) :: PHICU, PHIDU, PHIUU
881  DOUBLE PRECISION, INTENT(IN) :: PHICV, PHIDV, PHIUV
882  DOUBLE PRECISION, INTENT(IN) :: PHICW, PHIDW, PHIUW
883 !
884  DOUBLE PRECISION, INTENT(OUT) :: DWFE, DWFN, DWFT
885 ! discretization method
886  INTEGER, INTENT(IN) :: DISCR
887 
888 ! Local variables
889 !---------------------------------------------------------------------//
890 !
891  DOUBLE PRECISION :: PHI_C
892 ! cell widths for QUICKEST
893  DOUBLE PRECISION oDXc, oDXuc, oDYc, oDYuc, oDZc, oDZuc
894  DOUBLE PRECISION CF, DZUC
895 ! Indices
896  INTEGER :: I, J, K
897 !---------------------------------------------------------------------//
898 
899  SELECT CASE (discr) !first order upwinding
900  CASE (:1)
901  dwfe=0d0
902  dwfn=0d0
903  dwft=0d0
904 
905 
906  CASE (2) !Superbee
907  phi_c = phi_c_of(phiuu,phicu,phidu)
908  dwfe = superbee(phi_c)
909  phi_c = phi_c_of(phiuv,phicv,phidv)
910  dwfn = superbee(phi_c)
911  IF (do_k) THEN
912  phi_c = phi_c_of(phiuw,phicw,phidw)
913  dwft = superbee(phi_c)
914  ENDIF
915 
916 
917  CASE (3) !SMART
918  phi_c = phi_c_of(phiuu,phicu,phidu)
919  dwfe = smart(phi_c)
920  phi_c = phi_c_of(phiuv,phicv,phidv)
921  dwfn = smart(phi_c)
922  IF (do_k) THEN
923  phi_c = phi_c_of(phiuw,phicw,phidw)
924  dwft = smart(phi_c)
925  ENDIF
926 
927 
928  CASE (4) !ULTRA-QUICK
929  i=i_of(ijk)
930  phi_c = phi_c_of(phiuu,phicu,phidu)
931  cf = abs(uu)*dt*odx_e(i)
932  dwfe = ultra_quick(phi_c,cf)
933 
934  j=j_of(ijk)
935  phi_c = phi_c_of(phiuv,phicv,phidv)
936  cf = abs(vv)*dt*ody_n(j)
937  dwfn = ultra_quick(phi_c,cf)
938 
939  IF (do_k) THEN
940  k=k_of(ijk)
941  phi_c = phi_c_of(phiuw,phicw,phidw)
942  cf = abs(ww)*dt*ox(i)*odz_t(k)
943  dwft = ultra_quick(phi_c,cf)
944  ENDIF
945 
946 
947  CASE(5) !QUICKEST
948  i=i_of(ijk)
949  IF (uu >= zero) THEN
950  odxc = odx(i)
951  odxuc = odx_e(im1(i))
952  ELSE
953  odxc = odx(ip1(i))
954  odxuc = odx_e(ip1(i))
955  ENDIF
956  phi_c = phi_c_of(phiuu,phicu,phidu)
957  cf = abs(uu)*dt*odx_e(i)
958  dwfe = quickest(phi_c,cf,odxc,odxuc,odx_e(i))
959 
960  j=j_of(ijk)
961  IF (vv >= zero) THEN
962  odyc = ody(j)
963  odyuc = ody_n(jm1(j))
964  ELSE
965  odyc = ody(jp1(j))
966  odyuc = ody_n(jp1(j))
967  ENDIF
968  phi_c = phi_c_of(phiuv,phicv,phidv)
969  cf = abs(vv)*dt*ody_n(j)
970  dwfn = quickest(phi_c,cf,odyc,odyuc,ody_n(j))
971 
972  IF (do_k) THEN
973  k=k_of(ijk)
974  IF (ww >= zero) THEN
975  odzc = odz(k)
976  odzuc = odz_t(km1(k))
977  ELSE
978  odzc = odz(kp1(k))
979  dzuc = odz_t(kp1(k))
980  ENDIF
981  phi_c = phi_c_of(phiuw,phicw,phidw)
982  cf = abs(ww)*dt*ox(i)*odz_t(k)
983  dwft = quickest(phi_c,cf,odzc,odzuc,odz_t(k))
984  ENDIF
985 
986 
987  CASE (6) !MUSCL
988  phi_c = phi_c_of(phiuu,phicu,phidu)
989  dwfe = muscl(phi_c)
990  phi_c = phi_c_of(phiuv,phicv,phidv)
991  dwfn = muscl(phi_c)
992  IF (do_k) THEN
993  phi_c = phi_c_of(phiuw,phicw,phidw)
994  dwft = muscl(phi_c)
995  ENDIF
996 
997  CASE (7) !Van Leer
998  phi_c = phi_c_of(phiuu,phicu,phidu)
999  dwfe = vanleer(phi_c)
1000  phi_c = phi_c_of(phiuv,phicv,phidv)
1001  dwfn = vanleer(phi_c)
1002  IF (do_k) THEN
1003  phi_c = phi_c_of(phiuw,phicw,phidw)
1004  dwft = vanleer(phi_c)
1005  ENDIF
1006 
1007 
1008  CASE (8) !Minmod
1009  phi_c = phi_c_of(phiuu,phicu,phidu)
1010  dwfe = minmod(phi_c)
1011  phi_c = phi_c_of(phiuv,phicv,phidv)
1012  dwfn = minmod(phi_c)
1013  IF (do_k) THEN
1014  phi_c = phi_c_of(phiuw,phicw,phidw)
1015  dwft = minmod(phi_c)
1016  ENDIF
1017 
1018 
1019  CASE DEFAULT !Error
1020  WRITE (*,*) 'DISCRETIZE = ', discr, ' not supported.'
1021  CALL mfix_exit(mype)
1022 
1023  END SELECT
1024  RETURN
1025  END SUBROUTINE dw
1026 
1027 
1028 
1029 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1030 ! C
1031 ! Function: Xsi_func C
1032 ! Purpose: Special function for xsi that should be similar to: C
1033 ! xsi(v,dw) = merge( v >= 0, dwf, one-dwf) C
1034 ! C
1035 ! Slight difference when v is exactly zero but may not be C
1036 ! significant. C
1037 ! C
1038 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1039  DOUBLE PRECISION FUNCTION xsi_func(XXXv,XXXdwf)
1040 
1041  IMPLICIT NONE
1042  DOUBLE PRECISION, INTENT(IN) :: XXXv, XXXdwf
1043  xsi_func = (sign(1d0,(-xxxv))+1d0)/(2d0) + &
1044  sign(1d0,xxxv)*xxxdwf
1045  END FUNCTION xsi_func
1046 
1047  END MODULE xsi
indices::ip1
integer, dimension(:), allocatable ip1
Definition: indices_mod.f:50
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
xsi::cxs
subroutine cxs(INCR, DISCR, U, V, WW, PHI, XSI_E, XSI_N, XSI_T)
Definition: xsi_mod.f:579
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
error_manager::finl_err_msg
subroutine finl_err_msg
Definition: error_manager_mod.f:216
geometry::ody
double precision, dimension(:), allocatable ody
Definition: geometry_mod.f:116
run::shear
logical shear
Definition: run_mod.f:175
geometry::odx
double precision, dimension(:), allocatable odx
Definition: geometry_mod.f:114
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
indices::im1
integer, dimension(:), allocatable im1
Definition: indices_mod.f:50
error_manager::ival
Definition: error_manager_mod.f:17
discretization::chi_smart
double precision function chi_smart(PHI_C, Chi)
Definition: discretization_mod.f:66
run::dt
double precision dt
Definition: run_mod.f:51
discretization::superbee
double precision function superbee(PHI_C)
Definition: discretization_mod.f:20
xsi::calc_xsi
subroutine calc_xsi(DISCR, PHI, U, V, W, XSI_E, XSI_N, XSI_T, incr)
Definition: xsi_mod.f:23
discretization::quickest
double precision function quickest(PHI_C, CF, ODXC, ODXUC, ODXCD)
Definition: discretization_mod.f:147
xsi
Definition: xsi_mod.f:15
discretization::muscl
double precision function muscl(PHI_C)
Definition: discretization_mod.f:174
vshear::vsh
double precision, dimension(:), allocatable vsh
Definition: vshear_mod.f:3
discretization::phi_c_of
double precision function phi_c_of(PHI_U, PHI_C, PHI_D)
Definition: discretization_mod.f:342
xsi::dw
subroutine dw(UU, VV, WW, IJK, PHICU, PHIDU, PHIUU, PHICV, PHIDV, PHIUV, PHICW, PHIDW, PHIUW, DWFE, DWFN, DWFT, DISCR)
Definition: xsi_mod.f:857
error_manager::init_err_msg
subroutine init_err_msg(CALLER)
Definition: error_manager_mod.f:171
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
geometry::ody_n
double precision, dimension(:), allocatable ody_n
Definition: geometry_mod.f:123
sendrecv::send_recv
Definition: sendrecv_mod.f:75
exit::mfix_exit
subroutine mfix_exit(myID, normal_termination)
Definition: exit.f:5
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
indices::jm1
integer, dimension(:), allocatable jm1
Definition: indices_mod.f:51
exit
Definition: exit.f:2
discretization::ultra_quick
double precision function ultra_quick(PHI_C, CF)
Definition: discretization_mod.f:116
geometry::ox
double precision, dimension(:), allocatable ox
Definition: geometry_mod.f:140
indices::jp1
integer, dimension(:), allocatable jp1
Definition: indices_mod.f:51
discretization::minmod
double precision function minmod(PHI_C)
Definition: discretization_mod.f:265
discretization::vanleer
double precision function vanleer(PHI_C)
Definition: discretization_mod.f:248
discretization::central_scheme
double precision function central_scheme(PHI_C)
Definition: discretization_mod.f:289
indices::kp1
integer, dimension(:), allocatable kp1
Definition: indices_mod.f:52
chischeme::chi_t
double precision, dimension(:), allocatable chi_t
Definition: chischeme_mod.f:30
geometry::xlength
double precision xlength
Definition: geometry_mod.f:33
run
Definition: run_mod.f:13
vshear::vshe
double precision, dimension(:), allocatable vshe
Definition: vshear_mod.f:4
param
Definition: param_mod.f:2
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
indices::km1
integer, dimension(:), allocatable km1
Definition: indices_mod.f:52
compar::mype
integer mype
Definition: compar_mod.f:24
geometry::odz
double precision, dimension(:), allocatable odz
Definition: geometry_mod.f:118
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
chischeme::chi_n
double precision, dimension(:), allocatable chi_n
Definition: chischeme_mod.f:29
error_manager::err_msg
character(len=line_length), dimension(line_count) err_msg
Definition: error_manager_mod.f:29
discretization::smart
double precision function smart(PHI_C)
Definition: discretization_mod.f:42
error_manager
Definition: error_manager_mod.f:8
discretization::chi_muscl
double precision function chi_muscl(PHI_C, Chi)
Definition: discretization_mod.f:198
xsi::xsi_func
double precision function xsi_func(XXXv, XXXdwf)
Definition: xsi_mod.f:1040
geometry::odz_t
double precision, dimension(:), allocatable odz_t
Definition: geometry_mod.f:125
chischeme::chi_e
double precision, dimension(:), allocatable chi_e
Definition: chischeme_mod.f:28
chischeme::chi_flag
logical chi_flag
Definition: chischeme_mod.f:34
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
error_manager::flush_err_msg
subroutine flush_err_msg(DEBUG, HEADER, FOOTER, ABORT, LOG, CALL_TREE)
Definition: error_manager_mod.f:305
run::v_sh
double precision v_sh
Definition: run_mod.f:177