File: N:\mfix\model\kintheory_u_s.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_EXPLICIT_MOM_SOURCE_S                              C
4     !  Purpose: Determine any additional momentum source terms that are    C
5     !  calculated explicitly here                                          C
6     !                                                                      C
7     !                                                                      C
8     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
9     
10           SUBROUTINE CALC_EXPLICIT_MOM_SOURCE_S()
11     
12     !-----------------------------------------------
13     ! Modules
14     !-----------------------------------------------
15     
16           USE param1, only: zero
17     
18     ! kinetic theories
19           USE run, only: kt_type_enum
20           USE run, only: ia_2005
21     
22     ! number of solids phases
23           USE physprop, only: smax
24     
25     ! solids source term
26           USE kintheory, only: ktmom_u_s, ktmom_v_s, ktmom_w_s
27     
28           IMPLICIT NONE
29     !-----------------------------------------------
30     ! Local variables
31     !-----------------------------------------------
32     ! Solids phase index
33           INTEGER :: M
34     !-----------------------------------------------
35     ! Additional interphase interaction terms that arise from kinetic theory
36           DO M = 1, SMAX
37              KTMOM_U_s(:,M) = ZERO
38              KTMOM_V_S(:,M) = ZERO
39              KTMOM_W_S(:,M) = ZERO
40     
41              IF (KT_TYPE_ENUM == IA_2005) THEN
42                 CALL CALC_IA_MOMSOURCE_U_S(M)
43                 CALL CALC_IA_MOMSOURCE_V_S(M)
44                 CALL CALC_IA_MOMSOURCE_W_S(M)
45              ENDIF
46           ENDDO
47     
48           RETURN
49           END SUBROUTINE CALC_EXPLICIT_MOM_SOURCE_S
50     
51     
52     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
53     !                                                                      C
54     !  Subroutine: CALC_IA_MOMSOURCE_U_S                                   C
55     !  Purpose: Determine source terms for U_S momentum equation arising   C
56     !           from IA kinetic theory constitutive relations for stress   C
57     !           and solids-solids drag (collisional momentum source)       C
58     !                                                                      C
59     !  Literature/Document References:                                     C
60     !    Idir, Y.H., "Modeling of the multiphase mixture of particles      C
61     !      using the kinetic theory approach," PhD Thesis, Illinois        C
62     !      Institute of Technology, Chicago, Illinois, 2004                C
63     !    Iddir, Y.H., & H. Arastoopour, "Modeling of multitype particle    C
64     !      flow using the kinetic theory approach," AIChE J., Vol 51,      C
65     !      No 6, June 2005                                                 C
66     !                                                                      C
67     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
68     
69           SUBROUTINE CALC_IA_MOMSOURCE_U_S(M)
70     
71     !-----------------------------------------------
72     ! Modules
73     !-----------------------------------------------
74           USE param1, only: half, zero
75           USE constant, only: pi
76     
77     ! trace of D_s at i, j, k
78           USE visc_s, only: trD_s
79     
80     ! number of solids phases
81           USE physprop, only: mmax
82     
83     ! x,y,z-components of solids velocity
84           USE fldvar, only: u_s, v_s, w_s
85     ! particle diameter, bulk density, material density
86           USE fldvar, only: d_p, rop_s, ro_s
87     ! granular temperature
88           USE fldvar, only: theta_m, ep_s
89     ! dilute threshold
90           USE toleranc, only: dil_ep_s
91     
92           Use kintheory
93     
94           USE geometry
95           USE indices
96           USE compar
97           USE fun_avg
98           USE functions
99     
100           IMPLICIT NONE
101     !-----------------------------------------------
102     ! Dummy arguments
103     !-----------------------------------------------
104     ! solids phase index
105           INTEGER, INTENT(IN) :: M
106     !-----------------------------------------------
107     ! Local variables
108     !-----------------------------------------------
109     ! Temporary variable
110           DOUBLE PRECISION :: STRESS_TERMS, DRAG_TERMS
111     ! Indices
112           INTEGER :: I, J, JM, K, KM, IJK, IJKE, IPJK, IP, IMJK, IJKN,&
113                      IJKNE, IJKS, IJKSE, IPJMK, IJMK, IJKT, IJKTE,&
114                      IJKB, IJKBE, IJKM, IPJKM, &
115                      IJPK, IJKP, IM, IJKW
116     ! solids phase index
117           INTEGER :: L
118     ! Viscosity values for stress calculations
119           DOUBLE PRECISION :: MU_sL_pE, MU_sL_pW, MU_sL_pN, MU_sL_pS, MU_sL_pT,&
120                               MU_sL_pB, Mu_sL_p
121     ! Bulk viscosity values for calculations
122           DOUBLE PRECISION :: XI_sL_pE, XI_sL_pW, LAMBDA_sL_pE, LAMBDA_sL_pW
123     ! Variables for drag calculations
124           DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL, NU_PM_pE, NU_PM_pW, NU_PM_p, &
125                               NU_PL_pE, NU_PL_pW, NU_PL_p, T_PM_pE, T_PM_pW, &
126                               T_PL_pE, T_PL_pW, Fnu_s_p, FT_sM_p, FT_sL_p
127     ! Average volume fraction
128           DOUBLE PRECISION :: EPSA
129     ! Source terms (Surface)
130           DOUBLE PRECISION :: ssux, ssuy, ssuz, ssx, ssy, ssz, ssbv
131     ! Source terms (Volumetric)
132           DOUBLE PRECISION :: tauzz, DS1, DS2, DS3, DS4, DS1plusDS2
133     !-----------------------------------------------
134     
135     ! section largely based on tau_u_g:
136     
137           DO IJK = IJKSTART3, IJKEND3
138     
139     ! Skip walls where some values are undefined.
140              IF(WALL_AT(IJK)) cycle
141     
142              D_PM = D_P(IJK,M)
143              M_PM = (Pi/6d0)*D_PM**3 * RO_S(IJK,M)
144     
145              I = I_OF(IJK)
146              IJKE = EAST_OF(IJK)
147              EPSA = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
148              IF ( .NOT.SIP_AT_E(IJK) .AND. EPSA>DIL_EP_S) THEN
149     
150                 IP = IP1(I)
151                 IM = Im1(I)
152                 IJKW  = WEST_OF(IJK)
153                 J = J_OF(IJK)
154                 JM = JM1(J)
155                 K = K_OF(IJK)
156                 KM = KM1(K)
157                 IPJK = IP_OF(IJK)
158                 IMJK = IM_OF(IJK)
159                 IJMK = JM_OF(IJK)
160                 IJKM = KM_OF(IJK)
161                 IPJMK = JM_OF(IPJK)
162                 IPJKM = IP_OF(IJKM)
163     
164                 IJKN = NORTH_OF(IJK)
165                 IJKS = SOUTH_OF(IJK)
166                 IJKNE = EAST_OF(IJKN)
167                 IJKSE = EAST_OF(IJKS)
168     
169                 IJKT = TOP_OF(IJK)
170                 IJKB = BOTTOM_OF(IJK)
171                 IJKTE = EAST_OF(IJKT)
172                 IJKBE = EAST_OF(IJKB)
173     
174     ! additional required quantities:
175                 IJPK = JP_OF(IJK)
176                 IJKP = KP_OF(IJK)
177     
178     ! initialize variable
179                 STRESS_TERMS = ZERO
180                 DRAG_TERMS = ZERO
181     
182                 DO L = 1,MMAX
183                    IF (M .ne. L) THEN
184     
185     !--------------------- Sources from Stress Terms ---------------------
186     ! Surface Forces
187     ! standard shear stress terms (i.e. ~diffusion)
188                       MU_sL_pE = MU_sL_ip(IJKE,M,L)
189                       MU_sL_pW = MU_sL_ip(IJK,M,L)
190                       SSUX = MU_sL_pE*(U_S(IPJK,L)-U_S(IJK,L))*AYZ_U(IJK)*ODX(IP)&
191                            -MU_sL_PW*(U_S(IJK,L)-U_S(IMJK,L))*AYZ_U(IMJK)*ODX(I)
192     
193                       MU_sL_pN = AVG_X_H(AVG_Y_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKN,M,L), J),&
194                            AVG_Y_H(MU_sL_ip(IJKE,M,L),MU_sL_ip(IJKNE,M,L), J), I)
195                       MU_sL_pS = AVG_X_H(AVG_Y_H(MU_sL_ip(IJKS,M,L),MU_sL_ip(IJK,M,L), JM),&
196                            AVG_Y_H(MU_sL_ip(IJKSE,M,L),MU_sL_ip(IJKE,M,L), JM), I)
197                       SSUY = MU_sL_pN*(U_S(IJPK,L)-U_S(IJK,L))*AXZ_U(IJK)*ODY_N(J)&
198                            -MU_sL_pS*(U_S(IJK,L)-U_S(IJMK,L))*AXZ_U(IJMK)*ODY_N(JM)
199     
200                       MU_sL_pT = AVG_X_H(AVG_Z_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKT,M,L),K),&
201                            AVG_Z_H(MU_sL_ip(IJKE,M,L),MU_sL_ip(IJKTE,M,L),K),I)
202                       MU_sL_pB = AVG_X_H(AVG_Z_H(MU_sL_ip(IJKB,M,L),MU_sL_ip(IJK,M,L),KM),&
203                            AVG_Z_H(MU_sL_ip(IJKBE,M,L),MU_sL_ip(IJKE,M,L),KM),I)
204                       SSUZ = MU_sL_pT*(U_S(IJKP,L)-U_S(IJK,L))*AXY_U(IJK)*ODZ_T(K)*OX_E(I)&
205                            -MU_sL_pB*(U_S(IJK,L)-U_S(IJKM,L))*AXY_U(IJKM)*ODZ_T(KM)*OX_E(I)
206     
207     ! bulk viscosity term
208                       XI_sL_pE = XI_sL_ip(IJKE,M,L)
209                       XI_sL_pW = XI_sL_ip(IJK,M,L)
210                       LAMBDA_sL_pE = -(2.d0/3.d0)*Mu_sL_pE + XI_sL_pE
211                       LAMBDA_sL_pW = -(2.d0/3.d0)*Mu_sL_pW + XI_sL_pW
212                       SSBV = (LAMBDA_sL_pE*TRD_S(IJKE,L)-LAMBDA_sL_pW*TRD_S(IJK,L))*AYZ(IJK)
213     
214     ! off diagonal shear stress terms
215                       SSX = SSUX
216                       SSY = MU_sL_pN*(V_S(IPJK,L)-V_S(IJK,L))*AXZ_U(IJK)*ODX_E(I)&
217                            -MU_sL_pS*(V_S(IPJMK,L)-V_S(IJMK,L))*AXZ_U(IJMK)*ODX_E(I)
218                       SSZ = MU_sL_pT*(W_S(IPJK,L)-W_S(IJK,L))*AXY_U(IJK)*ODX_E(I)&
219                            -MU_sL_pB*(W_S(IPJKM,L)-W_S(IJKM,L))*AXY_U(IJKM)*ODX_E(I)
220     
221     ! special terms for cylindrical coordinates
222                       IF (CYLINDRICAL) THEN
223     
224     ! modify Ssz: (1/x) (d/dz) (tau_zz)
225     !                    integral of (1/x) (d/dz) (mu*(-w/x))
226     !                    (normally part of tau_u_s) - explicit
227                          SSZ = SSZ - (MU_sL_pT*(W_S(IPJK,L)+W_S(IJK,L))*HALF*OX_E(I)*AXY_U(IJK)&
228                               -MU_sL_pB*(W_S(IPJKM,L)+W_S(IJKM,L))*HALF*OX_E(I)*AXY_U(IJKM))
229     
230     ! tau_zz/x terms: (tau_zz/x)
231     !                    integral of -(2mu/x)*(1/x)*dw/dz
232     !                    (normally part of tau_u_s) - explicit
233                          MU_sL_p = AVG_X(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKE,M,L),I)
234                          tauzz = -2.d0*MU_sL_p*OX_E(I)*HALF*(&
235                               ((W_S(IPJK,L)-W_S(IPJKM,L))*OX(IP)*ODZ(K))+&
236                               ((W_S(IJK,L)-W_S(IJKM,L))*OX(I)*ODZ(K)) &
237                               ) * VOL_U(IJK)
238     
239     !                    integral of -(2mu/x)*(1/x)*u
240     !                    (normally part of source_u_s)
241                          tauzz = tauzz + (-2.d0*MU_sL_p*OX_E(I)*OX_E(I)*&
242                               U_S(IJK,L)*VOL_U(IJK))
243                       ELSE
244                          tauzz = ZERO
245                       ENDIF
246     !--------------------- End Sources from Stress Term ---------------------
247     
248     
249     !--------------------- Sources from Momentum Source Term ---------------------
250     ! Momentum source associated with the difference in the gradients in
251     ! number density of solids phase m and all other solids phases
252                       D_PL = D_P(IJK,L)
253                       M_PL = (PI/6.d0)*D_PL**3 * RO_S(IJK,L)
254     
255                       NU_PM_pE = ROP_S(IJKE,M)/M_PM
256                       NU_PM_pW = ROP_S(IJK,M)/M_PM
257                       NU_PM_p = AVG_X(NU_PM_pW,NU_PM_pE,I)
258     
259                       NU_PL_pE = ROP_S(IJKE,L)/M_PL
260                       NU_PL_pW = ROP_S(IJK,L)/M_PL
261                       NU_PL_p = AVG_X(NU_PL_pW,NU_PL_pE,I)
262     
263                       Fnu_s_p = AVG_X(Fnu_s_ip(IJK,M,L),Fnu_s_ip(IJKE,M,L),I)
264                       DS1 = Fnu_s_p*NU_PL_p*(NU_PM_pE-NU_PM_pW)*ODX_E(I)
265                       DS2 = -Fnu_s_p*NU_PM_p*(NU_PL_pE-NU_PL_pW)*ODX_E(I)
266                       DS1plusDS2 = DS1 + DS2
267     
268     ! Momentum source associated with the gradient in granular
269     ! temperature of species M
270                       T_PM_pE = Theta_M(IJKE,M)
271                       T_PM_pW = Theta_M(IJK,M)
272     
273                       FT_sM_p = AVG_X(FT_sM_ip(IJK,M,L),FT_sM_ip(IJKE,M,L),I)
274                       DS3 = FT_sM_p*(T_PM_pE-T_PM_pW)*ODX_E(I)
275     
276     ! Momentum source associated with the gradient in granular
277     ! temperature of species L
278                       T_PL_pE = Theta_M(IJKE,L)
279                       T_PL_pW = Theta_M(IJK,L)
280     
281                       FT_sL_p = AVG_X(FT_sL_ip(IJK,M,L),FT_sL_ip(IJKE,M,L),I)
282                       DS4 = FT_sL_p*(T_PL_pE-T_PL_pW)*ODX_E(I)
283     !--------------------- End Sources from Momentum Source Term ---------------------
284     
285     
286     ! Add the terms
287                       STRESS_TERMS = STRESS_TERMS + SSUX + SSUY + SSUZ + &
288                            SSBV + SSX + SSY + SSZ + tauzz
289                       DRAG_TERMS = DRAG_TERMS + (DS1plusDS2+DS3+DS4)*VOL_U(IJK)
290     
291                    ELSE ! if m .ne. L
292     ! for m = l all stress terms should already be handled in existing routines
293     ! for m = l all drag terms should become zero
294                       STRESS_TERMS = STRESS_TERMS + ZERO
295                       DRAG_TERMS = DRAG_TERMS + ZERO
296     
297                    ENDIF ! if m .ne. L
298                 ENDDO     ! over L
299     
300                 KTMOM_U_S(IJK,M) = STRESS_TERMS + DRAG_TERMS
301              ELSE
302                 KTMOM_U_S(IJK,M) = ZERO
303              ENDIF     ! dilute
304           ENDDO        ! over ijk
305     
306           RETURN
307           END SUBROUTINE CALC_IA_MOMSOURCE_U_S
308     
309     
310     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
311     !                                                                      !
312     !  Subroutine: COLL_MOMENTUM_COEFF_IA                                  !
313     !  Purpose: Compute collisional momentum source terms betweem solids   !
314     !  phase M and solids phase L using Iddir Arastoopour (2005) kinetic   !
315     !  theory model that are not proportional to the relative velocity     !
316     !  between the two phases.  Specifically, terms proportional to the    !
317     !  gradient in number density and gradient in temperature              !
318     !                                                                      !
319     !  Literature/Document References:                                     !
320     !    Iddir, Y.H., "Modeling of the multiphase mixture of particles     !
321     !       using the kinetic theory approach," PhD Thesis, Illinois       !
322     !       Institute of Technology, Chicago, Illinois, 2004               !
323     !    Iddir, Y.H., & H. Arastoopour, "Modeling of Multitype particle    !
324     !      flow using the kinetic theory approach," AIChE J., Vol 51,      !
325     !      no. 6, June 2005                                                !
326     !                                                                      !
327     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
328     
329           SUBROUTINE COLL_MOMENTUM_COEFF_IA(L, M)
330     
331     !-----------------------------------------------
332     ! Modules
333     !-----------------------------------------------
334           USE compar
335           USE constant
336           USE drag
337           USE fldvar
338           USE functions
339           USE geometry
340           USE indices
341           USE kintheory
342           USE param1
343           USE physprop
344           USE rdf
345           USE sendrecv
346     
347           IMPLICIT NONE
348     !-----------------------------------------------
349     ! Dummy arguments
350     !-----------------------------------------------
351     ! Solids phase index
352           INTEGER, INTENT(IN) :: L
353           INTEGER, INTENT(IN) :: M
354     !-----------------------------------------------
355     ! Local variables
356     !-----------------------------------------------
357     ! Indices
358           INTEGER :: IJK
359     ! Particle diameters
360           DOUBLE PRECISION :: D_PM, D_PL
361     ! Sum of particle diameters
362           DOUBLE PRECISION :: DPSUM
363     !
364           DOUBLE PRECISION :: M_PM, M_PL, MPSUM, DPSUMo2, NU_PL, NU_PM
365           DOUBLE PRECISION :: Ap_lm, Dp_lm, Bp_lm
366           DOUBLE PRECISION :: R0p_lm, R3p_lm, R4p_lm, R10p_lm
367           DOUBLE PRECISION :: Fnus_ip, FTsM_ip, FTsL_ip, F_common_term
368     !-----------------------------------------------
369     
370           DO IJK = ijkstart3, ijkend3
371              IF (.NOT.WALL_AT(IJK)) THEN
372     
373                 IF (M == L) THEN
374                    Fnu_s_ip(IJK,M,L) = ZERO
375                    FT_sM_ip(IJK,M,L) = ZERO
376                    FT_sL_ip(IJK,M,L) = ZERO
377     
378                 ELSE
379                    D_PM = D_P(IJK,M)
380                    D_PL = D_P(IJK,L)
381                    DPSUM = D_PL + D_PM
382                    M_PM = (Pi/6.d0) * D_PM**3 *RO_S(IJK,M)
383                    M_PL = (Pi/6.d0) * D_PL**3 *RO_S(IJK,L)
384                    MPSUM = M_PM + M_PL
385                    DPSUMo2 = DPSUM/2.d0
386                    NU_PM = ROP_S(IJK,M)/M_PM
387                    NU_PL = ROP_S(IJK,L)/M_PL
388     
389                    IF(Theta_m(IJK,M) > ZERO .AND. Theta_m(IJK,L) > ZERO) THEN
390     
391                       Ap_lm = (M_PM*Theta_m(IJK,L)+M_PL*Theta_m(IJK,M))/&
392                             2.d0
393                       Bp_lm = (M_PM*M_PL*(Theta_m(IJK,L)-Theta_m(IJK,M) ))/&
394                            (2.d0*MPSUM)
395                       Dp_lm = (M_PL*M_PM*(M_PM*Theta_m(IJK,M)+M_PL*Theta_m(IJK,L) ))/&
396                            (2.d0*MPSUM*MPSUM)
397     
398                       R0p_lm = ( 1.d0/( Ap_lm**1.5 * Dp_lm**2.5 ) )+ &
399                                 ( (15.d0*Bp_lm*Bp_lm)/( 2.d0* Ap_lm**2.5 * Dp_lm**3.5 ) )+&
400                                 ( (175.d0*(Bp_lm**4))/( 8.d0*Ap_lm**3.5 * Dp_lm**4.5 ) )
401     
402                       R3p_lm = ( 1.d0/( (Ap_lm**1.5)*(Dp_lm**3.5) ) )+&
403                                 ( (21.d0*Bp_lm*Bp_lm)/( 2.d0 * Ap_lm**2.5 * Dp_lm**4.5 ) )+&
404                                 ( (315.d0*Bp_lm**4)/( 8.d0 * Ap_lm**3.5 *Dp_lm**5.5 ) )
405     
406                       R4p_lm = ( 3.d0/( Ap_lm**2.5 * Dp_lm**3.5 ) )+&
407                                 ( (35.d0*Bp_lm*Bp_lm)/( 2.d0 * Ap_lm**3.5 * Dp_lm**4.5 ) )+&
408                                 ( (441.d0*Bp_lm**4)/( 8.d0 * Ap_lm**4.5 * Dp_lm**5.5 ) )
409     
410                       R10p_lm = ( 1.d0/( Ap_lm**2.5 * Dp_lm**2.5 ) )+&
411                                 ( (25.d0*Bp_lm*Bp_lm)/( 2.d0* Ap_lm**3.5 * Dp_lm**3.5 ) )+&
412                                 ( (1225.d0*Bp_lm**4)/( 24.d0* Ap_lm**4.5 * Dp_lm**4.5 ) )
413     
414                       F_common_term = (DPSUMo2*DPSUMo2/4.d0)*(M_PM*M_PL/MPSUM)*&
415                              G_0(IJK,M,L)*(1.d0+C_E)*(M_PM*M_PL)**1.5
416     
417     ! Momentum source associated with the difference in the gradients in
418     ! number density of solids phase m and all other solids phases
419                       Fnus_ip = F_common_term*(PI*DPSUMo2/12.d0)*R0p_lm*&
420                              (Theta_m(IJK,M)*Theta_m(IJK,L))**2.5
421     
422     ! Momentum source associated with the gradient in granular temperature
423     ! of solid phase M
424                       FTsM_ip = F_common_term*NU_PM*NU_PL*DPSUMo2*PI*&
425                                 (Theta_m(IJK,M)**1.5 * Theta_m(IJK,L)**2.5) *&
426                                 (  (-1.5d0/12.d0*R0p_lm)+&
427                                 Theta_m(IJK,L)/16.d0*(  (-M_PM*R10p_lm) - &
428                                 ((5.d0*M_PL*M_PL*M_PM/(192.d0*MPSUM*MPSUM))*R3p_lm)+&
429                                 ((5.d0*M_PM*M_PL)/(96.d0*MPSUM)*R4p_lm*Bp_lm)  )  )
430     
431     ! Momentum source associated with the gradient in granular temperature
432     ! of solid phase L ! no need to recompute (sof Aug 30 2006)
433                       FTsL_ip = F_common_term*NU_PM*NU_PL*DPSUMo2*PI*&
434                                (Theta_m(IJK,L)**1.5 * Theta_m(IJK,M)**2.5) *&
435                                 (  (1.5d0/12.d0*R0p_lm)+&
436                                 Theta_m(IJK,M)/16.d0*(  (M_PL*R10p_lm)+&
437                                 (5.d0*M_PM*M_PM*M_PL/(192.d0*MPSUM*MPSUM)*R3p_lm)+&
438                                 (5.d0*M_PM*M_PL/(96.d0*MPSUM) *R4p_lm*Bp_lm)  )  )
439     
440     
441                       Fnu_s_ip(IJK,M,L) = Fnus_ip
442     
443     ! WARNING: the following two terms have caused some convergence problems
444     ! earlier. Set them to ZERO for debugging in case of converegence
445     ! issues. (sof)
446                       FT_sM_ip(IJK,M,L) = FTsM_ip  ! ZERO
447                       FT_sL_ip(IJK,M,L) = FTsL_ip  ! ZERO
448                    ELSE
449                       Fnu_s_ip(IJK,M,L) = ZERO
450                       FT_sM_ip(IJK,M,L) = ZERO
451                       FT_sL_ip(IJK,M,L) = ZERO
452                    ENDIF
453                 ENDIF
454              ENDIF
455           ENDDO
456     
457           RETURN
458           END SUBROUTINE COLL_MOMENTUM_COEFF_IA
459     
460