MFIX  2016-1
ghdmassflux.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: GHDMASSFLUX C
4 ! Purpose: Calculate the species mass flux 3-components of Joi at cellC
5 ! faces to compute species velocities and source terms in T C
6 ! equation. C
7 ! C
8 ! Author: S. Benyahia Date: 13-MAR-09 C
9 ! Reviewer: Date: C
10 ! C
11 ! Literature/Document References: C
12 ! C. Hrenya handnotes and Garzo, Hrenya, Dufty papers (PRE, 2007) C
13 ! C
14 ! Variables modified: JoiX JoiY JoiZ C
15 ! C
16 ! Local variables: all terms in mass flux C
17 ! C
18 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19 
20  SUBROUTINE ghdmassflux ()
21 
22 !-----------------------------------------------
23 ! Modules
24 !-----------------------------------------------
25  USE param
26  USE param1
27  USE geometry
28  USE compar
29  USE fldvar
30  USE indices
31  USE visc_s
32  USE ghdtheory
33  USE physprop
34  USE run
35  USE constant
36  USE toleranc
37  USE drag
38  USE fun_avg
39  USE functions
40  IMPLICIT NONE
41 !-----------------------------------------------
42 ! Local variables
43 !-----------------------------------------------
44 ! Index
45  INTEGER :: IJK, I, J, K
46 ! Index
47  INTEGER :: IJKE, IJKN, IJKT
48 ! Solids phase
49  INTEGER :: M, L
50 ! number densities to compute del(Nj)
51  DOUBLE PRECISION NjC, NjE, NjN, NjT
52 
53 ! mass of species
54  DOUBLE PRECISION :: Mi, Mj
55 ! number density of species
56  DOUBLE PRECISION :: Ni
57 ! mixture density and temperature at cell faces
58  DOUBLE PRECISION :: ropsE, ropsN, ropsT, ThetaE, ThetaN, ThetaT
59  DOUBLE PRECISION :: EPSA
60 ! transport coefficient at cell faces
61  DOUBLE PRECISION :: DiTE, DiTN, DiTT
62  DOUBLE PRECISION :: DijE, DijN, DijT
63  DOUBLE PRECISION :: DijFE, DijFN, DijFT
64 ! Terms in the calculation of Joi-X,Y,Z
65  DOUBLE PRECISION :: ordinDiffTermX, ordinDiffTermY, ordinDiffTermZ
66  DOUBLE PRECISION :: massMobilityTermX, massMobilityTermY, &
67  massMobilityTermZ
68  DOUBLE PRECISION :: massMobilityTermXvelupdate, &
69  massMobilityTermYvelupdate, &
70  massMobilityTermZvelupdate
71  DOUBLE PRECISION :: thermalDiffTermX, thermalDiffTermY, &
72  thermalDiffTermZ
73  DOUBLE PRECISION :: ropsme, ropsmn, ropsmt
74 
75  DOUBLE PRECISION :: addtermx, &
76  addtermy, addtermz
77  DOUBLE PRECISION :: massMobilityTermNoDragX, &
78  massMobilityTermNoDragY, &
79  massMobilityTermNoDragZ
80  DOUBLE PRECISION :: DiTE_H, DiTE_A, DiTN_H, DiTN_A, DiTT_H, DiTT_A
81  DOUBLE PRECISION :: DijE_H, DijE_A, DijN_H, DijN_A, DijT_H, DijT_A
82  DOUBLE PRECISION :: DijFE_H, DijFE_A, DijFN_H, DijFN_A, DijFT_H, DijFT_A
83 !-----------------------------------------------
84 ! Function subroutines
85 !-----------------------------------------------
86 
87  DO m = 1, smax
88  DO 200 ijk = ijkstart3, ijkend3
89  i = i_of(ijk)
90  j = j_of(ijk)
91  k = k_of(ijk)
92 
93  IF ( fluid_at(ijk) ) THEN
94  mi = (pi/6.d0)*d_p(ijk,m)**3 * ro_s(ijk,m)
95  ni = rop_s(ijk,m) / mi
96 
97  ijke = east_of(ijk)
98  ijkn = north_of(ijk)
99  ijkt = top_of(ijk)
100 
101 ! mixture density and temperature evaluated at cell faces
102  ropse = avg_x(rop_s(ijk,mmax),rop_s(ijke,mmax),i)
103  ropsn = avg_y(rop_s(ijk,mmax),rop_s(ijkn,mmax),j)
104  ropst = avg_z(rop_s(ijk,mmax),rop_s(ijkt,mmax),k)
105 
106  thetae = avg_x(theta_m(ijk,mmax),theta_m(ijke,mmax),i)
107  thetan = avg_y(theta_m(ijk,mmax),theta_m(ijkn,mmax),j)
108  thetat = avg_z(theta_m(ijk,mmax),theta_m(ijkt,mmax),k)
109 
110 ! Thermal diffusion evaluated at cell faces (all used transport coef.
111 ! will be evaluated this way)
112  dite_h = avg_x_s(dit(ijk,m)*rop_s(ijk,mmax)/theta_m(ijk,mmax),&
113  dit(ijke,m)*rop_s(ijke,mmax)/theta_m(ijke,mmax),i)
114  dite_a = avg_x(dit(ijk,m)*rop_s(ijk,mmax)/theta_m(ijk,mmax),&
115  dit(ijke,m)*rop_s(ijke,mmax)/theta_m(ijke,mmax),i)
116 
117  ditn_h = avg_y_s(dit(ijk,m)*rop_s(ijk,mmax)/theta_m(ijk,mmax),&
118  dit(ijkn,m)*rop_s(ijkn,mmax)/theta_m(ijkn,mmax),j)
119  ditn_a = avg_y(dit(ijk,m)*rop_s(ijk,mmax)/theta_m(ijk,mmax),&
120  dit(ijkn,m)*rop_s(ijkn,mmax)/theta_m(ijkn,mmax),j)
121 
122  ditt_h = avg_z_s(dit(ijk,m)*rop_s(ijk,mmax)/theta_m(ijk,mmax),&
123  dit(ijkt,m)*rop_s(ijkt,mmax)/theta_m(ijkt,mmax),k)
124  ditt_a = avg_z(dit(ijk,m)*rop_s(ijk,mmax)/theta_m(ijk,mmax),&
125  dit(ijkt,m)*rop_s(ijkt,mmax)/theta_m(ijkt,mmax),k)
126 
127  IF(m .eq. 1)THEN
128  IF(min(abs(dite_h),abs(dite_a)) .eq. abs(dite_h))THEN
129  dite = dite_h
130  dit_harme(ijk) = .true.
131  ELSE
132  dite = dite_a
133  dit_harme(ijk) = .false.
134  ENDIF
135 
136  IF(min(abs(ditn_h),abs(ditn_a)) .eq. abs(ditn_h))THEN
137  ditn = ditn_h
138  dit_harmn(ijk) = .true.
139  ELSE
140  ditn = ditn_a
141  dit_harmn(ijk) = .false.
142  ENDIF
143 
144  IF(min(abs(ditt_h),abs(ditt_a)) .eq. abs(ditt_h))THEN
145  ditt = ditt_h
146  dit_harmt(ijk) = .true.
147  ELSE
148  ditt = ditt_a
149  dit_harmt(ijk) = .false.
150  ENDIF
151  ELSE
152  IF(dit_harme(ijk))THEN
153  dite = dite_h
154  ELSE
155  dite = dite_a
156  ENDIF
157 
158  IF(dit_harmn(ijk))THEN
159  ditn = ditn_h
160  ELSE
161  ditn = ditn_a
162  ENDIF
163 
164  IF(dit_harmt(ijk))THEN
165  ditt = ditt_h
166  ELSE
167  ditt = ditt_a
168  ENDIF
169  ENDIF ! end if/else M=1
170 
171 ! initializing variables for summation over L
172  ordindifftermx = zero
173  ordindifftermy = zero
174  ordindifftermz = zero
175  massmobilitytermx = zero
176  massmobilitytermy = zero
177  massmobilitytermz = zero
178  massmobilitytermxvelupdate = zero
179  massmobilitytermyvelupdate = zero
180  massmobilitytermzvelupdate = zero
181  addtermx = zero
182  addtermy = zero
183  addtermz = zero
184  massmobilitytermnodragx = zero
185  massmobilitytermnodragy = zero
186  massmobilitytermnodragz = zero
187 
188  DO l = 1, smax
189  mj = (pi/6.d0)*d_p(ijk,l)**3 * ro_s(ijk,l)
190 
191  njc = rop_s(ijk,l) / mj
192  nje = rop_s(ijke,l) / mj
193  njn = rop_s(ijkn,l) / mj
194  njt = rop_s(ijkt,l) / mj
195 
196  IF((rop_s(ijk,mmax)/ro_s(ijk,m) > dil_ep_s) .and. &
197  (rop_s(ijke,mmax)/ro_s(ijk,m) > dil_ep_s))THEN
198  dije_h = avg_x_s(dij(ijk,m,l)*mi*mj/rop_s(ijk,mmax),&
199  dij(ijke,m,l)*mi*mj/rop_s(ijke,mmax),i)
200  dije_a = avg_x(dij(ijk,m,l)*mi*mj/rop_s(ijk,mmax),&
201  dij(ijke,m,l)*mi*mj/rop_s(ijke,mmax),i)
202 
203  IF(m .eq. 1)THEN
204  IF(min(abs(dije_h),abs(dije_a)) .eq. abs(dije_h))THEN
205  dije = dije_h
206  dij_harme(ijk,l) = .true.
207  ELSE
208  dije = dije_a
209  dij_harme(ijk,l) = .false.
210  ENDIF
211  ELSE
212  IF(dij_harme(ijk,l))THEN
213  dije = dije_h
214  ELSE
215  dije = dije_a
216  ENDIF
217  ENDIF
218  ELSE
219  dije = zero
220  ENDIF
221 
222  IF((rop_s(ijk,mmax)/ro_s(ijk,m) > dil_ep_s) .and. &
223  (rop_s(ijkn,mmax)/ro_s(ijk,m) > dil_ep_s))THEN
224  dijn_h = avg_y_s(dij(ijk,m,l)*mi*mj/rop_s(ijk,mmax),&
225  dij(ijkn,m,l)*mi*mj/rop_s(ijkn,mmax),j)
226  dijn_a = avg_y(dij(ijk,m,l)*mi*mj/rop_s(ijk,mmax),&
227  dij(ijkn,m,l)*mi*mj/rop_s(ijkn,mmax),j)
228  IF(m .eq. 1)THEN
229  IF(min(abs(dijn_h),abs(dijn_a)) .eq. abs(dijn_h))THEN
230  dijn = dijn_h
231  dij_harmn(ijk,l) = .true.
232  ELSE
233  dijn = dijn_a
234  dij_harmn(ijk,l) = .false.
235  ENDIF
236  ELSE
237  IF(dij_harmn(ijk,l))THEN
238  dijn = dijn_h
239  ELSE
240  dijn = dijn_a
241  ENDIF
242  ENDIF
243  ELSE
244  dijn = zero
245  ENDIF
246 
247  IF((rop_s(ijk,mmax)/ro_s(ijk,m) > dil_ep_s) .and. &
248  (rop_s(ijkt,mmax)/ro_s(ijk,m) > dil_ep_s))THEN
249  dijt_h = avg_z_s(dij(ijk,m,l)*mi*mj/rop_s(ijk,mmax),&
250  dij(ijkt,m,l)*mi*mj/rop_s(ijkt,mmax),k)
251  dijt_a = avg_z(dij(ijk,m,l)*mi*mj/rop_s(ijk,mmax),&
252  dij(ijkt,m,l)*mi*mj/rop_s(ijkt,mmax),k)
253  IF(m .eq. 1)THEN
254  IF(min(abs(dijt_h),abs(dijt_a)) .eq. abs(dijt_h))THEN
255  dijt = dijt_h
256  dij_harmt(ijk,l) = .true.
257  ELSE
258  dijt = dijt_a
259  dij_harmt(ijk,l) = .false.
260  ENDIF
261  ELSE
262  IF(dij_harmt(ijk,l))THEN
263  dijt = dijt_h
264  ELSE
265  dijt = dijt_a
266  ENDIF
267  ENDIF
268  ELSE
269  dijt = zero
270  ENDIF
271 
272 
273  dijfe_h = avg_x_s(dijf(ijk,m,l),dijf(ijke,m,l),i)
274  dijfe_a = avg_x(dijf(ijk,m,l),dijf(ijke,m,l),i)
275  dijfn_h = avg_y_s(dijf(ijk,m,l),dijf(ijkn,m,l),j)
276  dijfn_a = avg_y(dijf(ijk,m,l),dijf(ijkn,m,l),j)
277  dijft_h = avg_z_s(dijf(ijk,m,l),dijf(ijkt,m,l),k)
278  dijft_a = avg_z(dijf(ijk,m,l),dijf(ijkt,m,l),k)
279 
280  IF(m .eq. 1)THEN
281  IF(min(abs(dijfe_h),abs(dijfe_a)) .eq. abs(dijfe_h))THEN
282  dijfe = dijfe_h
283  dijf_harme(ijk,l) = .true.
284  ELSE
285  dijfe = dijfe_a
286  dijf_harme(ijk,l) = .false.
287  ENDIF
288 
289  IF(min(abs(dijfn_h),abs(dijfn_a)) .eq. abs(dijfn_h))THEN
290  dijfn = dijfn_h
291  dijf_harmn(ijk,l) = .true.
292  ELSE
293  dijfn = dijfn_a
294  dijf_harmn(ijk,l) = .false.
295  ENDIF
296 
297  IF(min(abs(dijft_h),abs(dijft_a)) .eq. abs(dijft_h))THEN
298  dijft = dijft_h
299  dijf_harmt(ijk,l) = .true.
300  ELSE
301  dijft = dijft_a
302  dijf_harmt(ijk,l) = .false.
303  ENDIF
304  ELSE
305  IF(dijf_harme(ijk,l))THEN
306  dijfe = dijfe_h
307  ELSE
308  dijfe = dijfe_a
309  ENDIF
310 
311  IF(dijf_harmn(ijk,l))THEN
312  dijfn = dijfn_h
313  ELSE
314  dijfn = dijfn_a
315  ENDIF
316 
317  IF(dijf_harmt(ijk,l))THEN
318  dijft = dijft_h
319  ELSE
320  dijft = dijft_a
321  ENDIF
322  ENDIF ! end if/else (M=1)
323 
324  ordindifftermx = ordindifftermx + dije * (nje - njc) * odx_e(i)
325  ordindifftermy = ordindifftermy + dijn * (njn - njc) * ody_n(j)
326  ordindifftermz = ordindifftermz + dijt * (njt - njc) * (ox_e(i)*odz_t(k))
327 
328  massmobilitytermx = massmobilitytermx + dijfe * fix(ijk,l)
329  massmobilitytermy = massmobilitytermy + dijfn * fiy(ijk,l)
330  massmobilitytermz = massmobilitytermz + dijft * fiz(ijk,l)
331 
332  massmobilitytermxvelupdate = massmobilitytermxvelupdate + dijfe * fixvel(ijk,l)
333  massmobilitytermyvelupdate = massmobilitytermyvelupdate + dijfn * fiyvel(ijk,l)
334  massmobilitytermzvelupdate = massmobilitytermzvelupdate + dijft * fizvel(ijk,l)
335 
336  massmobilitytermnodragx = massmobilitytermnodragx + dijfe * fiminusdragx(ijk,l)
337  massmobilitytermnodragy = massmobilitytermnodragy + dijfn * fiminusdragy(ijk,l)
338  massmobilitytermnodragz = massmobilitytermnodragz + dijft * fiminusdragz(ijk,l)
339 
340  ENDDO ! end do (l=1,smax)
341 
342 
343  thermaldifftermx = dite * ( theta_m(ijke,mmax) - theta_m(ijk,mmax) ) * odx_e(i)
344  thermaldifftermy = ditn * ( theta_m(ijkn,mmax) - theta_m(ijk,mmax) ) * ody_n(j)
345  thermaldifftermz = ditt * ( theta_m(ijkt,mmax) - theta_m(ijk,mmax) ) * (ox_e(i)*odz_t(k))
346 
347  joix(ijk,m) = -(ordindifftermx + thermaldifftermx + massmobilitytermx)
348  joiy(ijk,m) = -(ordindifftermy + thermaldifftermy + massmobilitytermy)
349  joiz(ijk,m) = -(ordindifftermz + thermaldifftermz + massmobilitytermz)
350 
351  joiminusdragx(ijk,m) = (ordindifftermx + thermaldifftermx + massmobilitytermnodragx)
352  joiminusdragy(ijk,m) = (ordindifftermy + thermaldifftermy + massmobilitytermnodragy)
353  joiminusdragz(ijk,m) = (ordindifftermz + thermaldifftermz + massmobilitytermnodragz)
354 
355  ropsme=avg_x(rop_s(ijk,m),rop_s(ijke,m),i)
356  ropsmn=avg_y(rop_s(ijk,m),rop_s(ijkn,m),j)
357  ropsmt=avg_z(rop_s(ijk,m),rop_s(ijkt,m),k)
358 
359  deltau(ijk,m) = -(ordindifftermx+thermaldifftermx+massmobilitytermxvelupdate)
360  deltav(ijk,m) = -(ordindifftermy+thermaldifftermy+massmobilitytermyvelupdate)
361  deltaw(ijk,m) = -(ordindifftermz+thermaldifftermz+massmobilitytermzvelupdate)
362 
363 ! set fluxes to zero in case very dilute conditions
364  epsa = avg_x(rop_s(ijk,m),rop_s(ijke,m),i) / ro_s(ijk,m)
365  IF(epsa <= zero_ep_s) joix(ijk,m) = zero
366 
367  epsa = avg_y(rop_s(ijk,m),rop_s(ijkn,m),j) / ro_s(ijk,m)
368  IF(epsa <= zero_ep_s) joiy(ijk,m) = zero
369 
370  epsa = avg_z(rop_s(ijk,m),rop_s(ijkt,m),k) / ro_s(ijk,m)
371  IF(epsa <= zero_ep_s) joiz(ijk,m) = zero
372 
373 ! set flux components to zero in case of walls
374  IF (ip_at_e(ijk)) joix(ijk,m) = zero
375  IF (ip_at_n(ijk)) joiy(ijk,m) = zero
376  IF (ip_at_t(ijk)) joiz(ijk,m) = zero
377 
378  ELSE
379  joix(ijk,m) = zero
380  joiy(ijk,m) = zero
381  joiz(ijk,m) = zero
382  ENDIF ! end if/else (fluid_at(ijk))
383 
384  200 CONTINUE ! outer IJK loop
385  ENDDO ! for m = 1,smax
386 
387  RETURN
388  END SUBROUTINE ghdmassflux
389 
390 
391 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
392 ! C
393 ! Subroutine: UpdateSpeciesVelocities C
394 ! Purpose: Update solids velocities at celll faces based on the C
395 ! formula Ui = U + Joi/(mi ni); also calculate averaged C
396 ! velocities for dilute conditions. C
397 ! C
398 ! Author: S. Benyahia Date: 6-MAR-09 C
399 ! Reviewer: Date: C
400 ! C
401 ! Literature/Document References: C
402 ! C. Hrenya handnotes and Garzo, Hrenya, Dufty papers (PRE, 2007) C
403 ! C
404 ! Variables modified: solid species velocity components: Us, Vs, Ws C
405 ! C
406 ! C
407 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
408 
409  SUBROUTINE updatespeciesvelocities ()
410 
411 !-----------------------------------------------
412 ! Modules
413 !-----------------------------------------------
414  USE param
415  USE param1
416  USE geometry
417  USE compar
418  USE fldvar
419  USE indices
420  USE is
421  USE drag
422  USE visc_s
423  USE ghdtheory
424  USE physprop
425  USE run
426  USE constant
427  USE toleranc
428  USE fun_avg
429  USE functions
430  IMPLICIT NONE
431 !-----------------------------------------------
432 ! Local variables
433 !-----------------------------------------------
434 ! Index
435  INTEGER :: IJK, I, J, K
436 ! Index
437  INTEGER :: IJKE, IJKN, IJKT, IJKW, IJKS, IJKB, IMJK, IJMK, &
438  IJKM
439 ! Solids phase
440  INTEGER :: SS
441 ! species density at cell faces
442  integer :: kk, maxFluxS
443  double precision :: epgN, rogN, mugN, Vg
444  double precision :: Ur(smax), vrelSq(smax), vel, velup(smax)
445  double precision :: maxFlux
446  double precision :: rosN(smax), dp(smax)
447  double precision :: DijN(smax,smax), JoiM(smax), &
448  DijN_H(smax,smax), DijN_A(smax,smax)
449  double precision :: beta_cell(smax), beta_ij_cell(smax,smax)
450 
451  integer :: ntrial
452  double precision tolx, tolf
453 !-----------------------------------------------
454 ! Function subroutines
455 !-----------------------------------------------
456 
457 ! for Newton method
458  ntrial = 100
459  tolx = 1d-14
460  tolf = 1d-14
461 
462  DO 200 ijk = ijkstart3, ijkend3
463  i = i_of(ijk)
464  j = j_of(ijk)
465  k = k_of(ijk)
466 
467  IF ( fluid_at(ijk) ) THEN
468 
469  ijke = east_of(ijk)
470  ijkw = west_of(ijk)
471  ijkn = north_of(ijk)
472  ijks = south_of(ijk)
473  ijkt = top_of(ijk)
474  ijkb = bottom_of(ijk)
475  imjk = im_of(ijk)
476  ijmk = jm_of(ijk)
477  ijkm = km_of(ijk)
478 
479 ! First Compute Us
480 ! ---------------------------------------------------------------->>>
481  IF(.NOT.ip_at_e(ijk) .OR. .NOT.sip_at_e(ijk)) THEN
482  IF(ro_g0 == zero) THEN ! special case of a granular flow
483  do ss = 1, smax
484  rosn(ss) = avg_x(rop_s(ijk,ss),rop_s(ijke,ss),i)/ ro_s(ijk,ss) ! this is ep_s
485  if(rosn(ss) > zero_ep_s) then
486  u_s(ijk,ss) = u_s(ijk,mmax) + joix(ijk,ss)/(rosn(ss)*ro_s(ijk,ss))
487  else
488  u_s(ijk,ss) = u_s(ijk,mmax) ! case of zero flux
489  endif
490  enddo
491  ELSE
492  vg = u_g(ijk) - u_s(ijk,mmax)
493  epgn = avg_x(ep_g(ijk),ep_g(ijke),i)
494  rogn = avg_x(rop_g(ijk),rop_g(ijke),i)
495  mugn = avg_x(mu_g(ijk),mu_g(ijke),i)
496 
497  do ss = 1, smax
498  ur(ss) = u_g(ijk)-u_s(ijk,ss)
499 ! vrel must not include Ur, which is being solved for iteratively and must be updated.
500  vrelsq(ss) = (v_g(ijk)-v_s(ijk,ss))**2 + (w_g(ijk)-w_s(ijk,ss))**2
501  rosn(ss) = avg_x(rop_s(ijk,ss),rop_s(ijke,ss),i)
502  velup(ss) = 0.d0
503  beta_cell(ss) = beta_cell_x(ijk,ss)
504  dp(ss) = d_p(ijk,ss)
505 
506  IF(drag_type_enum .eq. hys)THEN
507  joim(ss) = deltau(ijk,ss)
508  ELSE
509  joim(ss) = joiminusdragx(ijk,ss)
510  ENDIF
511 
512  do kk = 1, smax
513  dijn_h(ss,kk) = avg_x_s(dijf(ijk,ss,kk),dijf(ijke,ss,kk),i)
514  dijn_a(ss,kk) = avg_x(dijf(ijk,ss,kk),dijf(ijke,ss,kk),i)
515  if(dijf_harme(ijk,kk))THEN
516  dijn(ss,kk) = dijn_h(ss,kk)
517  ELSE
518  dijn(ss,kk) = dijn_a(ss,kk)
519  ENDIF
520  if(ss .eq. kk)then
521  beta_ij_cell(ss,kk)=0.d0
522  else
523  beta_ij_cell(ss,kk)=beta_ij_cell_x(ijk,ss,kk)
524  endif
525  enddo
526  enddo
527 
528  IF(drag_type_enum .eq. hys)THEN
529  vel=u_s(ijk,mmax)
530  CALL velocity_update(velup, smax, rosn, dijn, &
531  joim, beta_cell, beta_ij_cell,vel)
532  ELSE
533  CALL urnewt(ntrial, ur, smax, ijk, tolx, tolf, &
534  epgn, rogn, mugn, vg, vrelsq, rosn, dp, dijn, joim)
535  ENDIF
536 
537 ! species velocity and flux update.
538  do ss = 1, smax
539  IF(drag_type_enum .eq. hys)THEN
540  u_s(ijk,ss)=velup(ss)
541  joix(ijk,ss) = rosn(ss) * (u_s(ijk,ss)-u_s(ijk,mmax))
542  ELSE
543  u_s(ijk,ss) = u_g(ijk) - ur(ss)
544  joix(ijk,ss) = rosn(ss) * (u_s(ijk,ss)-u_s(ijk,mmax))
545  ENDIF
546  enddo
547  ENDIF ! for a granular case (no gas and no drag)
548  ENDIF ! end if (.not.ip_at_e or .not.sip_at_e)
549 
550  if(smax==2) then ! only for binary, how to implement for smax > 2?
551  joix(ijk,2)=-joix(ijk,1)
552  elseif(smax > 2) then
553  maxflux = joix(ijk,1)
554  maxfluxs = 1
555  do ss = 2, smax ! finding species with maximum flux in a cell
556  if( abs(joix(ijk,ss)) > abs(maxflux) ) then
557  maxflux = joix(ijk,ss)
558  maxfluxs = ss
559  endif
560  enddo
561  joix(ijk,maxfluxs) = 0d0 ! reset max. flux to zero
562  do ss = 1, smax ! re-calc species with max. flux to satisfy SUM(fluxes) = 0
563  if(ss /= maxfluxs) joix(ijk,maxfluxs) = joix(ijk,maxfluxs) - joix(ijk,ss)
564  enddo
565  endif
566 ! End Compute Us
567 ! ----------------------------------------------------------------<<<
568 
569 
570 
571 ! Now Compute Vs
572 ! ---------------------------------------------------------------->>>
573  IF (.NOT.ip_at_n(ijk) .OR. .NOT.sip_at_n(ijk)) THEN
574  IF(ro_g0 == zero) THEN ! special case of a granular flow
575  do ss = 1, smax
576  rosn(ss) = avg_y(rop_s(ijk,ss),rop_s(ijkn,ss),j)/ ro_s(ijk,ss) ! this is ep_s
577  if(rosn(ss) > zero_ep_s) then
578  v_s(ijk,ss) = v_s(ijk,mmax) + joiy(ijk,ss)/(rosn(ss)*ro_s(ijk,ss))
579  else
580  v_s(ijk,ss) = v_s(ijk,mmax) ! case of zero flux
581  endif
582  enddo
583  ELSE
584  vg = v_g(ijk) - v_s(ijk,mmax)
585  epgn = avg_y(ep_g(ijk),ep_g(ijkn),j)
586  rogn = avg_y(rop_g(ijk),rop_g(ijkn),j)
587  mugn = avg_y(mu_g(ijk),mu_g(ijkn),j)
588 
589  do ss = 1, smax
590  ur(ss) = v_g(ijk)-v_s(ijk,ss)
591 ! vrel must not include Ur, which is being solved for iteratively and must be updated.
592  vrelsq(ss) = (u_g(ijk)-u_s(ijk,ss))**2 + (w_g(ijk)-w_s(ijk,ss))**2
593  rosn(ss) = avg_y(rop_s(ijk,ss),rop_s(ijkn,ss),j)
594  velup(ss) = 0.d0
595  beta_cell(ss) = beta_cell_y(ijk,ss)
596  dp(ss) = d_p(ijk,ss)
597 
598  IF(drag_type_enum .eq. hys)THEN
599  joim(ss) = deltav(ijk,ss)
600  ELSE
601  joim(ss) = joiminusdragy(ijk,ss)
602  ENDIF
603 
604  do kk = 1, smax
605 
606  dijn_h(ss,kk) = avg_y_s(dijf(ijk,ss,kk),dijf(ijkn,ss,kk),j)
607  dijn_a(ss,kk) = avg_y(dijf(ijk,ss,kk),dijf(ijkn,ss,kk),j)
608 
609  if(dijf_harmn(ijk,kk))THEN
610  dijn(ss,kk) = dijn_h(ss,kk)
611  ELSE
612  dijn(ss,kk) = dijn_a(ss,kk)
613  ENDIF
614 
615  if(ss .eq. kk)then
616  beta_ij_cell(ss,kk)=0.d0
617  else
618  beta_ij_cell(ss,kk)=beta_ij_cell_y(ijk,ss,kk)
619  endif
620  enddo
621  enddo
622 
623  IF(drag_type_enum .eq. hys)THEN
624  vel=v_s(ijk,mmax)
625  CALL velocity_update(velup, smax, rosn, dijn, joim, &
626  beta_cell, beta_ij_cell,vel)
627  ELSE
628  CALL urnewt(ntrial, ur, smax, ijk, tolx, tolf, &
629  epgn, rogn, mugn, vg, vrelsq, rosn, dp, dijn, joim)
630  ENDIF
631 
632 ! species velocity and flux update
633  do ss = 1, smax
634  IF(drag_type_enum .eq. hys)THEN
635  v_s(ijk,ss)=velup(ss)
636  joiy(ijk,ss) = rosn(ss) * (v_s(ijk,ss)-v_s(ijk,mmax))
637  ELSE
638  v_s(ijk,ss) = v_g(ijk) - ur(ss)
639  joiy(ijk,ss) = rosn(ss) * (v_s(ijk,ss)-v_s(ijk,mmax))
640  ENDIF
641  enddo
642  ENDIF ! for a granular case (no gas and no drag)
643  ENDIF ! end if (.not.ip_at_n or .not.sip_at_n)
644 
645  if(smax==2) then ! only for binary, how to implement for smax > 2?
646  joiy(ijk,2)=-joiy(ijk,1)
647  elseif(smax > 2) then
648  maxflux = joiy(ijk,1)
649  maxfluxs = 1
650  do ss = 2, smax ! finding species with maximum flux in a cell
651  if( abs(joiy(ijk,ss)) > abs(maxflux) ) then
652  maxflux = joiy(ijk,ss)
653  maxfluxs = ss
654  endif
655  enddo
656  joiy(ijk,maxfluxs) = 0d0 ! reset max. flux to zero
657  do ss = 1, smax ! re-calc species with max. flux to satisfy SUM(fluxes) = 0
658  if(ss /= maxfluxs) joiy(ijk,maxfluxs) = joiy(ijk,maxfluxs) - joiy(ijk,ss)
659  enddo
660  endif
661 ! End Compute Vs
662 ! ----------------------------------------------------------------<<<
663 
664 
665 ! Finaly Compute Ws
666 ! ---------------------------------------------------------------->>>
667  IF(.NOT.no_k .AND. (.NOT.ip_at_t(ijk) .OR. .NOT.sip_at_t(ijk))) THEN
668  IF(ro_g0 == zero) THEN ! special case of a granular flow
669  do ss = 1, smax
670  rosn(ss) = avg_z(rop_s(ijk,ss),rop_s(ijkt,ss),k)/ ro_s(ijk,ss) ! this is ep_s
671  if(rosn(ss) > zero_ep_s) then
672  w_s(ijk,ss) = w_s(ijk,mmax) + joiz(ijk,ss)/(rosn(ss)*ro_s(ijk,ss))
673  else
674  w_s(ijk,ss) = w_s(ijk,mmax) ! case of zero flux
675  endif
676  enddo
677  ELSE
678  vg = w_g(ijk) - w_s(ijk,mmax)
679  epgn = avg_z(ep_g(ijk),ep_g(ijkt),k)
680  rogn = avg_z(rop_g(ijk),rop_g(ijkt),k)
681  mugn = avg_z(mu_g(ijk),mu_g(ijkt),k)
682 
683  do ss = 1, smax
684  ur(ss) = w_g(ijk)-w_s(ijk,ss)
685 ! vrel must not include Ur, which is being solved for iteratively and must be updated.
686  vrelsq(ss) = (u_g(ijk)-u_s(ijk,ss))**2 + (v_g(ijk)-v_s(ijk,ss))**2
687  rosn(ss) = avg_z(rop_s(ijk,ss),rop_s(ijkt,ss),k)
688  velup(ss) = 0.d0
689  beta_cell(ss) = beta_cell_z(ijk,ss)
690  dp(ss) = d_p(ijk,ss)
691 
692  IF(drag_type_enum .eq. hys)THEN
693  joim(ss) = deltaw(ijk,ss)
694  ELSE
695  joim(ss) = joiminusdragz(ijk,ss)
696  ENDIF
697 
698  do kk = 1, smax
699  dijn_h(ss,kk) = avg_z_s(dijf(ijk,ss,kk),dijf(ijkt,ss,kk),k)
700  dijn_a(ss,kk) = avg_z(dijf(ijk,ss,kk),dijf(ijkt,ss,kk),k)
701  if(dijf_harmt(ijk,kk))THEN
702  dijn(ss,kk) = dijn_h(ss,kk)
703  ELSE
704  dijn(ss,kk) = dijn_a(ss,kk)
705  ENDIF
706  if(ss .eq. kk)then
707  beta_ij_cell(ss,kk)=0.d0
708  else
709  beta_ij_cell(ss,kk)=beta_ij_cell_z(ijk,ss,kk)
710  endif
711  enddo
712  enddo
713 
714  IF(drag_type_enum .eq. hys)THEN
715  vel=w_s(ijk,mmax)
716  CALL velocity_update(velup, smax, rosn, dijn, &
717  joim, beta_cell, beta_ij_cell,vel)
718  ELSE
719  CALL urnewt(ntrial, ur, smax, ijk, tolx, tolf, &
720  epgn, rogn, mugn, vg, vrelsq, rosn, dp, dijn, joim)
721  ENDIF
722 
723 ! species velocity and flux update
724  do ss = 1, smax
725  IF(drag_type_enum .eq. hys)THEN
726  w_s(ijk,ss)=velup(ss)
727  joiz(ijk,ss) = rosn(ss) * (w_s(ijk,ss)-w_s(ijk,mmax))
728  ELSE
729  w_s(ijk,ss) = w_g(ijk) - ur(ss)
730  joiz(ijk,ss) = rosn(ss) * (w_s(ijk,ss)-w_s(ijk,mmax))
731  ENDIF
732  enddo
733  ENDIF ! for a granular case (no gas and no drag)
734  ENDIF ! end if (.not.ip_at_t or .not.sip_at_t)
735 
736  if(smax==2) then ! only for binary, how to implement for smax > 2?
737  joiz(ijk,2)=-joiz(ijk,1)
738  elseif(smax > 2) then
739  maxflux = joiz(ijk,1)
740  maxfluxs = 1
741  do ss = 2, smax ! finding species with maximum flux in a cell
742  if( abs(joiz(ijk,ss)) > abs(maxflux) ) then
743  maxflux = joiz(ijk,ss)
744  maxfluxs = ss
745  endif
746  enddo
747  joiz(ijk,maxfluxs) = 0d0 ! reset max. flux to zero
748  do ss = 1, smax ! re-calc species with max. flux to satisfy SUM(fluxes) = 0
749  if(ss /= maxfluxs) joiz(ijk,maxfluxs) = joiz(ijk,maxfluxs) - joiz(ijk,ss)
750  enddo
751  endif
752 ! End Compute Ws
753 ! ----------------------------------------------------------------<<<
754 
755 ! if .not. fluid_at(ijk), do nothing (keep old values of velocities).
756 
757  ENDIF ! Fluid_at
758  200 CONTINUE ! outer IJK loop
759 
760  RETURN
761  END SUBROUTINE updatespeciesvelocities
762 
763 
764 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
765 ! C
766 ! Subroutine: UrNEWT C
767 ! C
768 ! C
769 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
770 
771  subroutine urnewt(ntrial, x, s, ijk, tolx, tolf, epgN, rogN, &
772  mugn, vg, vrelsq, rosn, dp, dijn, joim)
773 
774  Implicit NONE
775 
776 !-----------------------------------------------
777 ! Dummy arguments
778 !-----------------------------------------------
779  integer, intent(in) :: s
780  integer, intent(in) :: ijk
781  integer, intent(in) :: ntrial
782  double precision, intent(in) :: tolx, tolf
783 
784  DOUBLE PRECISION, intent(inout) :: X(s)
785  double precision, intent(in) :: epgN, rogN, mugN, Vg
786  double precision, intent(in) :: vrelSq(s), rosN(s), dp(s)
787  double precision, intent(in) :: DijN(s,s), JoiM(s)
788 !-----------------------------------------------
789 ! Local parameters
790 !-----------------------------------------------
791  INTEGER :: NP
792  parameter(np=15) ! solves up to NP variable/equations;
793 !-----------------------------------------------
794 ! Local variables
795 !-----------------------------------------------
796  INTEGER :: I, K, INDX(s)
797  DOUBLE PRECISION :: D, ERRF, ERRX, FJAC(s,s), FVEC(s), P(s)
798 !-----------------------------------------------
799 
800  DO k = 1, ntrial
801  CALL ur_jacobi_eval(x, s, ijk, fvec, fjac, epgn, rogn, mugn, vg, &
802  vrelsq, rosn, dp, dijn, joim)
803  errf = 0d0
804  DO i = 1, s
805  errf = errf + dabs(fvec(i))
806  ENDDO
807 
808  IF(errf <= tolf) RETURN
809  DO i = 1, s
810  p(i) = -fvec(i) ! RHS of linear equations.
811  IF(rosn(i) .eq. 0.d0)THEN
812  p(i) = 0.d0
813  ENDIF
814  ENDDO
815 
816  CALL ludcmp(fjac, s, np, indx, d, 'UrNewt') ! solve system of s linear equations using
817  CALL lubksb(fjac, s, np, indx, p) ! LU decomposition
818 
819  errx = 0d0
820  DO i = 1, s
821  errx = errx + dabs(p(i))
822  x(i) = x(i) + p(i)
823  ENDDO
824  IF(errx <= tolx) RETURN
825  ENDDO ! for K trials
826 
827  RETURN
828  END SUBROUTINE urnewt
829 
830 
831 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
832 ! C
833 ! Subroutine: Ur_JACOBI_EVAL C
834 ! C
835 ! C
836 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
837 
838  SUBROUTINE ur_jacobi_eval(X, s, ijk, FVEC, FJAC, epgN, rogN, &
839  mugn, vg, vrelsq, rosn, dp, dijn, joim)
840  Implicit NONE
841 
842 !-----------------------------------------------
843 ! Dummy arguments
844 !-----------------------------------------------
845  integer, intent(in) :: s
846  integer, intent(in) :: ijk
847  double precision, intent(in) :: X(s)
848 ! vector function and matrix jacobian
849  DOUBLE PRECISION, INTENT(OUT) :: FVEC(s), FJAC(s,s)
850  DOUBLE PRECISION, INTENT(IN) :: epgN, rogN, mugN, Vg
851  DOUBLE PRECISION, INTENT(IN) :: vrelSq(s), rosN(s), dp(s)
852  DOUBLE PRECISION, INTENT(IN) :: DijN(s,s), JoiM(s)
853 !-----------------------------------------------
854 ! Local parameters
855 !-----------------------------------------------
856  double precision :: pi
857  parameter(pi=3.14159265458979323846d0)
858  double precision :: one
859  parameter(one=1.d0)
860  double precision :: zero
861  parameter(zero=0.d0)
862 !-----------------------------------------------
863 ! Local variables
864 !-----------------------------------------------
865 ! Indices
866  INTEGER :: I, J
867 ! various quantities
868  DOUBLE PRECISION :: RE_G, C_d, DgA, Vi, vrel
869  DOUBLE PRECISION :: FgsOni(s), dFgsdVi(s), sum(s)
870 !-----------------------------------------------
871 
872  DO i = 1, s
873  vrel = dsqrt(vrelsq(i) + x(i)**2)
874  vi = pi * dp(i)**3 / 6d0
875  re_g = dp(i)*vrel*rogn/mugn
876  IF(re_g <= 1000d0 .and. re_g> zero)THEN
877  c_d = (24.d0/re_g) * (one + 0.15d0 * re_g**0.687d0)
878  ELSE
879  c_d = 0.44d0
880  ENDIF
881  dga = 0.75d0 * c_d * vrel * rogn / (epgn**2.65d0 * dp(i))
882  fgsoni(i) = dga * vi ! this is the wen-Yu drag
883  IF(vrel == zero)THEN
884  dfgsdvi(i) = zero
885  ELSEIF(re_g <= 1000d0)THEN
886  dfgsdvi(i) = 1.8549d0*mugn*re_g**0.687d0*vi*dabs(x(i))/ &
887  (dp(i)**2*epgn**2.65d0*vrel**2)
888  ELSE
889  dfgsdvi(i) = 0.33d0*rogn*vi*dabs(x(i))/ &
890  (dp(i)*epgn**2.65d0*vrel)
891  ENDIF
892  ENDDO
893 
894 ! Start computing values of the function FVEC(i)
895  DO i = 1, s
896  sum(i) = zero
897  DO j = 1, s
898  sum(i) = sum(i) + dijn(i,j) * fgsoni(j) * x(j)
899  ENDDO
900  ENDDO
901  DO i = 1, s
902  fvec(i) = sum(i) - rosn(i)*x(i) + rosn(i)*vg + joim(i)
903  ENDDO
904 
905 !
906 ! Evaluation of functions FVEC(i) is done above, and their Jacobian is computed below.
907 !
908  DO i = 1, s
909  DO j = 1, s
910  if(j == i) THEN
911  fjac(i,j) = dijn(i,i) * (fgsoni(i) + dfgsdvi(i)*x(i)) - rosn(i)
912  IF(rosn(i) .eq. 0.d0)THEN
913  fjac(i,j) = 1.d0
914  ENDIF
915  else
916  fjac(i,j) = dijn(i,j) * (fgsoni(j) + dfgsdvi(j)*x(j))
917  endif
918  ENDDO ! for j
919  ENDDO ! for i
920 !
921 ! End of computing jacobian (J_ij).
922 !
923  RETURN
924  END SUBROUTINE ur_jacobi_eval
925 
926 
927 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
928 ! C
929 ! Subroutine: Velocity_update C
930 ! C
931 ! C
932 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
933 
934  SUBROUTINE velocity_update(FVEC, s, rosi, Diji, Joii, &
935  beta_celli, beta_ij_celli, &
936  velocity)
937  USE toleranc
938  USE fldvar
939  Implicit NONE
940 
941 !-----------------------------------------------
942 ! Dummy arguments
943 !-----------------------------------------------
944  integer, intent (in) :: s
945 ! various quantities
946  DOUBLE PRECISION, INTENT(IN) :: rosi(s), Diji(s,s), &
947  Joii(s), beta_celli(s), &
948  beta_ij_celli(s,s), velocity
949 ! vector function
950  DOUBLE PRECISION, INTENT(OUT) :: FVEC(s)
951 !-----------------------------------------------
952 ! Local parameters
953 !-----------------------------------------------
954  integer :: NP
955  parameter(np=15) ! solves up to NP variable/equations;
956 !-----------------------------------------------
957 ! Local variables
958 !-----------------------------------------------
959 ! Indices
960  INTEGER :: i, k, l
961  INTEGER :: indx(s)
962  DOUBLE PRECISION :: D
963 ! matrix jacobian
964  DOUBLE PRECISION :: FJAC(s,s)
965 !-----------------------------------------------
966 
967  DO i=1,s
968  fvec(i)= 0.d0
969  DO k=1,s
970  fjac(i,k) = 0.d0
971  ENDDO
972  ENDDO
973 
974  DO i=1,s
975 
976  IF(rosi(i) .ne. 0.d0)THEN
977  fvec(i) = rosi(i)*velocity+joii(i)
978  ELSE
979  fvec(i) = 0.d0
980  ENDIF
981 
982  DO l=1,s
983 
984  if(i .eq. l)then
985  fjac(i,l)=fjac(i,l)+rosi(i)
986  IF(rosi(i) .eq. 0.d0)THEN ! This is done to avoid a singular matrix if ros(i) is ZERO
987  fjac(i,l) = 1.d0
988  ENDIF
989  endif
990 
991  fjac(i,l)=fjac(i,l)-diji(i,l)*beta_celli(l)
992 
993  DO k=1,s
994 
995  if(l .ne. k)then
996  fjac(i,k)=fjac(i,k)+diji(i,l)*beta_ij_celli(l,k)
997  endif
998 
999  enddo
1000 
1001  enddo
1002  enddo
1003 
1004  CALL ludcmp(fjac, s, np, indx, d, 'velocity_update') ! solve system of s linear equations using
1005  CALL lubksb(fjac, s, np, indx, fvec) ! LU decomposition
1006 
1007  RETURN
1008 
1009  END SUBROUTINE velocity_update
1010 
ghdtheory::beta_ij_cell_z
double precision, dimension(:,:,:), allocatable beta_ij_cell_z
Definition: ghdtheory_mod.f:118
ghdtheory::joix
double precision, dimension(:,:), allocatable joix
Definition: ghdtheory_mod.f:30
fldvar::v_s
double precision, dimension(:,:), allocatable v_s
Definition: fldvar_mod.f:105
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
visc_s
Definition: visc_s_mod.f:1
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
fldvar::ep_g
double precision, dimension(:), allocatable ep_g
Definition: fldvar_mod.f:15
ghdtheory::beta_ij_cell_x
double precision, dimension(:,:,:), allocatable beta_ij_cell_x
Definition: ghdtheory_mod.f:112
ghdtheory::dit_harme
logical, dimension(:), allocatable dit_harme
Definition: ghdtheory_mod.f:121
geometry::ox_e
double precision, dimension(:), allocatable ox_e
Definition: geometry_mod.f:143
ghdtheory::dijf_harmt
logical, dimension(:,:), allocatable dijf_harmt
Definition: ghdtheory_mod.f:132
ghdtheory::dit_harmt
logical, dimension(:), allocatable dit_harmt
Definition: ghdtheory_mod.f:125
ghdtheory::beta_ij_cell_y
double precision, dimension(:,:,:), allocatable beta_ij_cell_y
Definition: ghdtheory_mod.f:115
ghdtheory::dij_harmn
logical, dimension(:,:), allocatable dij_harmn
Definition: ghdtheory_mod.f:137
velocity_update
subroutine velocity_update(FVEC, s, rosi, Diji, Joii, beta_celli, beta_ij_celli, velocity)
Definition: ghdmassflux.f:937
fldvar::w_s
double precision, dimension(:,:), allocatable w_s
Definition: fldvar_mod.f:117
ghdtheory::joiminusdragx
double precision, dimension(:,:), allocatable joiminusdragx
Definition: ghdtheory_mod.f:60
ghdmassflux
subroutine ghdmassflux()
Definition: ghdmassflux.f:21
ghdtheory::fiminusdragy
double precision, dimension(:,:), allocatable fiminusdragy
Definition: ghdtheory_mod.f:63
drag
Definition: drag_mod.f:11
updatespeciesvelocities
subroutine updatespeciesvelocities()
Definition: ghdmassflux.f:410
ghdtheory::deltaw
double precision, dimension(:,:), allocatable deltaw
Definition: ghdtheory_mod.f:82
ghdtheory::fix
double precision, dimension(:,:), allocatable fix
Definition: ghdtheory_mod.f:39
is
Definition: is_mod.f:11
ur_jacobi_eval
subroutine ur_jacobi_eval(X, s, ijk, FVEC, FJAC, epgN, rogN, mugN, Vg, vrelSq, rosN, dp, DijN, JoiM)
Definition: ghdmassflux.f:840
fldvar::u_s
double precision, dimension(:,:), allocatable u_s
Definition: fldvar_mod.f:93
ghdtheory::beta_cell_x
double precision, dimension(:,:), allocatable beta_cell_x
Definition: ghdtheory_mod.f:104
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
ghdtheory::fiz
double precision, dimension(:,:), allocatable fiz
Definition: ghdtheory_mod.f:45
geometry::ody_n
double precision, dimension(:), allocatable ody_n
Definition: geometry_mod.f:123
fun_avg
Definition: fun_avg_mod.f:1
ghdtheory::dij
double precision, dimension(:,:,:), allocatable dij
Definition: ghdtheory_mod.f:24
fldvar::d_p
double precision, dimension(:,:), allocatable d_p
Definition: fldvar_mod.f:57
ghdtheory::fixvel
double precision, dimension(:,:), allocatable fixvel
Definition: ghdtheory_mod.f:48
physprop::mmax
integer mmax
Definition: physprop_mod.f:19
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
physprop::ro_g0
double precision ro_g0
Definition: physprop_mod.f:59
ghdtheory::dij_harme
logical, dimension(:,:), allocatable dij_harme
Definition: ghdtheory_mod.f:135
toleranc::zero_ep_s
double precision, parameter zero_ep_s
Definition: toleranc_mod.f:15
fldvar::theta_m
double precision, dimension(:,:), allocatable theta_m
Definition: fldvar_mod.f:149
fldvar::v_g
double precision, dimension(:), allocatable v_g
Definition: fldvar_mod.f:99
ghdtheory::dit_harmn
logical, dimension(:), allocatable dit_harmn
Definition: ghdtheory_mod.f:123
ghdtheory::joiy
double precision, dimension(:,:), allocatable joiy
Definition: ghdtheory_mod.f:33
ghdtheory::beta_cell_y
double precision, dimension(:,:), allocatable beta_cell_y
Definition: ghdtheory_mod.f:107
fldvar::w_g
double precision, dimension(:), allocatable w_g
Definition: fldvar_mod.f:111
run
Definition: run_mod.f:13
param
Definition: param_mod.f:2
ludcmp
subroutine ludcmp(a, n, np, indx, d, calledFrom)
Definition: cooling_rate.f:287
toleranc::dil_ep_s
double precision, parameter dil_ep_s
Definition: toleranc_mod.f:24
fldvar::ro_s
double precision, dimension(:,:), allocatable ro_s
Definition: fldvar_mod.f:45
geometry::no_k
logical no_k
Definition: geometry_mod.f:28
ghdtheory::dit
double precision, dimension(:,:), allocatable dit
Definition: ghdtheory_mod.f:15
physprop::mu_g
double precision, dimension(:), allocatable mu_g
Definition: physprop_mod.f:68
urnewt
subroutine urnewt(ntrial, x, s, ijk, tolx, tolf, epgN, rogN, mugN, Vg, vrelSq, rosN, dp, DijN, JoiM)
Definition: ghdmassflux.f:773
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
fldvar::u_g
double precision, dimension(:), allocatable u_g
Definition: fldvar_mod.f:87
ghdtheory::dijf_harmn
logical, dimension(:,:), allocatable dijf_harmn
Definition: ghdtheory_mod.f:130
ghdtheory::deltau
double precision, dimension(:,:), allocatable deltau
Definition: ghdtheory_mod.f:76
ghdtheory::fiminusdragx
double precision, dimension(:,:), allocatable fiminusdragx
Definition: ghdtheory_mod.f:57
fldvar::rop_s
double precision, dimension(:,:), allocatable rop_s
Definition: fldvar_mod.f:51
physprop::smax
integer smax
Definition: physprop_mod.f:22
ghdtheory::fiy
double precision, dimension(:,:), allocatable fiy
Definition: ghdtheory_mod.f:42
ghdtheory::fiyvel
double precision, dimension(:,:), allocatable fiyvel
Definition: ghdtheory_mod.f:51
ghdtheory::joiminusdragy
double precision, dimension(:,:), allocatable joiminusdragy
Definition: ghdtheory_mod.f:66
lubksb
subroutine lubksb(a, n, np, indx, b)
Definition: cooling_rate.f:383
ghdtheory::dijf_harme
logical, dimension(:,:), allocatable dijf_harme
Definition: ghdtheory_mod.f:128
ghdtheory::dijf
double precision, dimension(:,:,:), allocatable dijf
Definition: ghdtheory_mod.f:18
ghdtheory::beta_cell_z
double precision, dimension(:,:), allocatable beta_cell_z
Definition: ghdtheory_mod.f:110
ghdtheory::fizvel
double precision, dimension(:,:), allocatable fizvel
Definition: ghdtheory_mod.f:54
constant::pi
double precision, parameter pi
Definition: constant_mod.f:158
ghdtheory::joiminusdragz
double precision, dimension(:,:), allocatable joiminusdragz
Definition: ghdtheory_mod.f:72
geometry::odz_t
double precision, dimension(:), allocatable odz_t
Definition: geometry_mod.f:125
ghdtheory::deltav
double precision, dimension(:,:), allocatable deltav
Definition: ghdtheory_mod.f:79
fldvar::rop_g
double precision, dimension(:), allocatable rop_g
Definition: fldvar_mod.f:38
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
ghdtheory::joiz
double precision, dimension(:,:), allocatable joiz
Definition: ghdtheory_mod.f:36
ghdtheory::dij_harmt
logical, dimension(:,:), allocatable dij_harmt
Definition: ghdtheory_mod.f:139
ghdtheory::fiminusdragz
double precision, dimension(:,:), allocatable fiminusdragz
Definition: ghdtheory_mod.f:69