1!-------------------------------------------------------------------------------
2
3! This file is part of Code_Saturne, a general-purpose CFD tool.
4!
5! Copyright (C) 1998-2021 EDF S.A.
6!
7! This program is free software; you can redistribute it and/or modify it under
8! the terms of the GNU General Public License as published by the Free Software
9! Foundation; either version 2 of the License, or (at your option) any later
10! version.
11!
12! This program is distributed in the hope that it will be useful, but WITHOUT
13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
15! details.
16!
17! You should have received a copy of the GNU General Public License along with
18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
19! Street, Fifth Floor, Boston, MA 02110-1301, USA.
20
21!-------------------------------------------------------------------------------
22
23subroutine cs_coal_noxst &
24 ( ncel   ,                                                               &
25   indpdf ,                                                               &
26   pdfm1  , pdfm2  , doxyd  , dfuel  , hrec ,                             &
27   f3m    , f4m    , f5m    , f6m    , f7m  , f8m , f9m ,                 &
28   fs3no  , fs4no  , yfs4no , enthox )
29
30!===============================================================================
31! FONCTION :
32! --------
33! CALCUL DE :
34!
35!           K1 exp(-E1/RT)   (conversion HCN en N2)
36!           K2 exp(-E2/RT)   (conversion HCN en NO)
37!           K3 exp(-E3/RT)   (NO thermique)
38!
39!----------------
40! Arguments
41!__________________.____._____.________________________________________________.
42! name             !type!mode ! role                                           !
43!__________________!____!_____!________________________________________________!
44! ncel             ! i  ! <-- ! number of cells                                !
45! indpdf           ! te ! <-- ! passage par les pdf                            !
46!__________________!____!_____!________________________________________________!
47
48!     TYPE : E (ENTIER), R (REEL), A (ALPHAMNUMERIQUE), T (TABLEAU)
49!            L (LOGIQUE)   .. ET TYPES COMPOSES (EX : TR TABLEAU REEL)
50!     MODE : <-- donnee, --> resultat, <-> Donnee modifiee
51!            --- tableau de travail
52!===============================================================================
53
54!==============================================================================
55! Module files
56!==============================================================================
57
58use paramx
59use numvar
60use optcal
61use cstphy
62use cstnum
63use entsor
64use parall
65use pointe
66use ppppar
67use ppthch
68use coincl
69use cpincl
70use ppincl
71use ppcpfu
72use field
73use cs_c_bindings
74
75!===============================================================================
76
77implicit none
78
79! Arguments
80
81integer          ncel
82integer          indpdf(ncel)
83
84double precision pdfm1(ncel) , pdfm2(ncel) , dfuel(ncel)
85double precision f3m(ncel)   , f4m(ncel)   , f5m(ncel)
86double precision f6m(ncel)   , f7m(ncel)   , f8m(ncel) ,f9m(ncel)
87double precision doxyd(ncel) , hrec(ncel)
88double precision fs3no(ncel) , fs4no(ncel) , yfs4no(ncel,ngazg)
89double precision enthox(ncel)
90!
91! VARIABLES LOCALES
92!
93integer          iel , icla , i
94integer          ipdf1 , ipdf2 , ipdf3
95integer          npart
96parameter        (npart = 200)
97
98double precision ee1,ee2,ee3,kk1,kk2,kk3
99double precision wmo2,tg,xo2,bb
100double precision xash,xch,xck,xmx2
101double precision bb1 , bb2 , bb3 , bb4
102double precision dirac , tfuel , tmpgaz , lrf ,lro
103double precision qqq , rrr  , sss , ttt , uuu
104double precision gt1,gt2,gt3,gt10,gt20,gt30
105double precision ts4,ts4num,ts4den
106double precision dgs
107double precision val(npart+1),tt(npart+1) , gs(npart+1) , yyo2(npart+1)
108!
109integer          icha,ige,numcha,mode
110double precision xxf , hhf , hfuel,tfs4ad,hfs4ad
111double precision xhf1,xhf2 , aa , den
112double precision coefe(ngazem),f1mc(ncharm),f2mc(ncharm)
113double precision yo2moy,yo2ox,yo2cb,yo24den,yo24num,yo2s4,yo2
114double precision toxyd,hoxyd , som , react , deltamol
115!
116integer          i300 ,i000 ,imini,i2500,i2600,i2700,i2800,i3000,i3500
117integer          i4000,i5000,imaxi,inok
118integer          nbpt , nbclip1 , nbclip2
119integer          nbclip30,nbclip31
120double precision ts4min,ts4max
121double precision somm  , sommax , sommin , ts4admin,ts4admax
122double precision yo2min,yo2max,yo2min1,yo2max1
123double precision yo2oxmin,yo2oxmax
124double precision toxmin,toxmax
125
126!LOCAL VARIABLES
127!===============
128!
129! Kinetic constants
130! -----------------
131double precision kk4,ee4,kk5,ee5,kkrb,eerb
132!
133! Indicators for integration on the pdf
134! -------------------------------------
135integer ipdf4,ipdf5
136
137type(pmapper_double_r1), dimension(:), allocatable :: cvar_xckcl, cvar_xchcl
138type(pmapper_double_r1), dimension(:), allocatable :: cvar_xnpcl, cvar_xwtcl
139type(pmapper_double_r1), dimension(:), allocatable :: cvar_f1m, cvar_f2m
140type(pmapper_double_r1), dimension(:), allocatable :: cpro_temp2
141double precision, dimension(:), pointer :: cpro_exp1, cpro_exp2, cpro_exp3
142double precision, dimension(:), pointer :: cpro_exp4, cpro_exp5, cpro_exprb
143double precision, dimension(:), pointer :: cpro_temp, cpro_yo2, cpro_mmel
144
145logical(kind=c_bool) :: log_active
146
147!===============================================================================
148! 1. Preliminary computations
149!===============================================================================
150
151! Molar masses
152!
153wmo2 = wmole(io2)
154!
155! pointers
156
157call field_get_val_s(ighcn1, cpro_exp1)
158call field_get_val_s(ighcn2, cpro_exp2)
159call field_get_val_s(ignoth, cpro_exp3)
160call field_get_val_s(ignh31, cpro_exp4)
161call field_get_val_s(ignh32, cpro_exp5)
162call field_get_val_s(igrb, cpro_exprb)
163
164! Parameters of Arrhenius laws
165
166kk1 = 3.0e12
167ee1 = 3.0e4
168kk2 = 1.2e10
169ee2 = 3.35e4
170kk3 = 3.4e12
171ee3 = 6.69e4
172
173kk4 = 4.1e6
174ee4 = 1.6e4
175kk5 = 1.8e8
176ee5 = 1.35e4
177kkrb= 2.7e12
178eerb= 9.467e3
179
180! Pour les termes, indicateur de calcul par PDF ou non
181!       = 1 --> passage par pdf
182!       = 0 --> on ne passe pas par les pdf
183
184ipdf1 = 0
185ipdf2 = 0
186ipdf3 = 1
187
188ipdf4 = 0
189ipdf5 = 0
190
191! Initialisation
192
193do iel = 1, ncel
194  cpro_exp1(iel)  = zero
195  cpro_exp2(iel)  = zero
196  cpro_exp3(iel)  = zero
197
198  cpro_exp4(iel)  = zero
199  cpro_exp5(iel)  = zero
200  cpro_exprb(iel) = zero
201enddo
202
203!===============================================================================
204! 2. CALCUL SANS LES PDF
205!===============================================================================
206
207call field_get_val_s(itemp, cpro_temp)
208call field_get_val_s(iym1(io2), cpro_yo2)
209call field_get_val_s(immel, cpro_mmel)
210
211do iel = 1, ncel
212
213  tg  = cpro_temp(iel)
214  yo2 = cpro_yo2(iel)
215  xo2 = yo2*cpro_mmel(iel)/wmo2
216
217  ! Reaction HCN + NO + 1/4 O2 ---> N2 + 1/2 H2O + CO
218
219  cpro_exp1(iel)  = kk1*exp(-ee1/tg)
220
221  ! Reaction HCN + 5/4 O2 --> NO + 1/2 H2O  + CO
222
223  if (xo2 .gt. 0.018d0) then
224    bb = 0.d0
225  else if (xo2 .lt. 0.0025d0) then
226    bb= 1.d0
227  else
228    bb=(0.018d0-xo2)/(0.018d0-0.0025d0)
229  endif
230  cpro_exp2(iel)  = kk2*exp(-ee2/tg)*(xo2**bb)
231
232  ! No thermique : Zeldovich
233
234  cpro_exp3(iel)  = kk3*exp(-ee3/tg)*(xo2**0.5d0)
235
236  ! Reaction NH3 + O2 --> NO + ...
237  cpro_exp4(iel)  = kk4*exp(-ee4/tg)*(xo2**bb)
238
239  ! Reaction NH3 + NO --> N2 + ...
240  cpro_exp5(iel)  = kk5*exp(-ee5/tg)
241
242  ! Reburning (Model de Chen)
243  cpro_exprb(iel) = kkrb*exp(-eerb/tg)
244
245enddo
246
247!===============================================================================
248! 3. Computation with the PDFs
249!===============================================================================
250!
251if (     ipdf1.eq.1 .or. ipdf2.eq.1 .or. ipdf3.eq.1 &
252    .or. ipdf4.eq.1 .or. ipdf5.eq.1) then
253
254  ! Arrays of pointers containing the field values for each class
255  ! (loop on cells outside loop on classes)
256  allocate(cvar_xckcl(nclacp), cvar_xchcl(nclacp))
257  allocate(cvar_xnpcl(nclacp), cvar_xwtcl(nclacp))
258  allocate(cvar_f1m(ncharb), cvar_f2m(ncharb))
259  allocate(cpro_temp2(nclacp))
260
261  do icla = 1, nclacp
262    call field_get_val_s(ivarfl(isca(ixck(icla))), cvar_xckcl(icla)%p)
263    call field_get_val_s(ivarfl(isca(ixch(icla))), cvar_xchcl(icla)%p)
264    call field_get_val_s(ivarfl(isca(inp(icla))), cvar_xnpcl(icla)%p)
265    call field_get_val_s(itemp2(icla), cpro_temp2(icla)%p)
266    if (ippmod(iccoal) .eq. 1) then
267      call field_get_val_s(ivarfl(isca(ixwt(icla))), cvar_xwtcl(icla)%p)
268    endif
269  enddo
270
271  do numcha = 1, ncharb
272    call field_get_val_s(ivarfl(isca(if1m(numcha))), cvar_f1m(numcha)%p)
273    call field_get_val_s(ivarfl(isca(if2m(numcha))), cvar_f2m(numcha)%p)
274  enddo
275
276  inok = 0
277  i300 = 0
278  i000 = 0
279  imini= 0
280  i2500= 0
281  i2600= 0
282  i2700= 0
283  i2800= 0
284  i3000= 0
285  i3500= 0
286  i4000= 0
287  i5000= 0
288  imaxi= 0
289  nbpt = 0
290  nbclip1 =0
291  nbclip2 = 0
292  ts4min= 1.d+20
293  ts4max=-1.d+20
294  sommin= 1.d+20
295  sommax=-1.d+20
296  ts4admin= 1.d+20
297  ts4admax=-1.d+20
298  toxmin = 1.d+20
299  toxmax =-1.d+20
300  yo2oxmin= 1.d+20
301  yo2oxmax=-1.d20
302
303  nbclip30 = 0
304  nbclip31 = 0
305  yo2min  = 1.d+20
306  yo2max  =-1.d+20
307  yo2min1 = 1.d+20
308  yo2max1 =-1.d+20
309
310  do iel=1, ncel
311
312    if (indpdf(iel).ne.1) cycle
313
314    if ((fs3no (iel).gt.fs4no(iel)) .and. (fs4no (iel).lt.1.d0)) then
315
316      ! Calcul de Yo2 dans l'oxydant
317      !           Yo2 en fs4
318
319      bb1 = max(0.d0      ,pdfm1(iel))
320      bb2 = min(fs3no(iel),pdfm2(iel))
321      if (bb2 .gt. bb1) then
322        lro = hrec(iel)
323      else
324        lro = 0.d0
325      endif
326      qqq = bb2**2 - bb1**2
327      rrr = bb2    - bb1
328      gt1 = lro*qqq/(2.d0*fs3no(iel))
329      gt2 = lro*rrr
330      gt3 = doxyd(iel)
331      yo2cb = 0.d0
332
333      if (cpro_yo2(iel) .gt. 0.d0) then
334        yo2ox = cpro_yo2(iel)/(-gt1+gt2+gt3)
335
336        yo2oxmin = min(yo2oxmin,yo2ox)
337        yo2oxmax = max(yo2oxmax,yo2ox)
338
339        yo2moy = cpro_yo2(iel)
340        dirac  =  dfuel(iel)*yo2cb + doxyd(iel)*yo2ox
341
342        bb1 = max(0.D0      ,pdfm1(iel))
343        bb2 = min(fs4no(iel),pdfm2(iel))
344        bb3 = max(fs4no(iel),pdfm1(iel))
345        bb4 = min(fs3no(iel),pdfm2(iel))
346        if (bb2 .gt. bb1) then
347          lro = hrec(iel)
348        else
349          lro = 0.d0
350        endif
351        if (bb4 .gt. bb3) then
352          lrf = hrec(iel)
353        else
354          lrf = 0.d0
355        endif
356
357        qqq = bb2**2 - bb1**2
358        rrr = bb2    - bb1
359        sss = bb4**2 - bb3**2
360        ttt = bb4    - bb3
361        uuu = fs4no(iel)-fs3no(iel)
362
363        gt1 = lro*qqq/(2.d0*fs4no(iel))
364        gt2 = lro*rrr
365
366        gt10= lrf*sss/(2.d0*uuu)
367        gt20= lrf*ttt*fs3no(iel)/uuu
368
369        yo24num = yo2moy - dirac + yo2ox*(gt1 -gt2)
370        yo24den = gt1+gt10-gt20
371
372        yo2s4  = yo24num/yo24den
373
374      else
375        yo2ox = 0.d0
376        yo2s4 = 0.d0
377      endif
378
379      ! Affichage et clipping
380
381      yo2min = min(yo2s4,yo2min)
382      yo2max = max(yo2s4,yo2max)
383      if (yo2s4 .lt. 0.d0) then
384         nbclip30 = nbclip30+1
385         yo2s4 = 0.d0
386      endif
387      if (yo2s4 .gt. yo2ox) then
388         nbclip31 = nbclip31+1
389         yo2s4 = yo2ox
390      endif
391      yo2min1 = min(yo2s4,yo2min1)
392      yo2max1 = max(yo2s4,yo2max1)
393
394      ! Calcul de Tfioul moyen
395      ! ----------------------
396
397      tfuel = 0.d0
398      xmx2  = 0.d0
399
400      do icla = 1, nclacp
401        xck  = cvar_xckcl(icla)%p(iel)
402        xch  = cvar_xchcl(icla)%p(iel)
403        xash = cvar_xnpcl(icla)%p(iel)*xmash(icla)
404        xmx2   = xmx2 + xch + xck + xash
405
406        ! Prise en compte de l'humidite
407
408        if (ippmod(iccoal) .eq. 1) then
409          xmx2 = xmx2+cvar_xwtcl(icla)%p(iel)
410        endif
411      enddo
412      if (xmx2 .gt. 0.d0) then
413        do icla=1,nclacp
414          tfuel = tfuel +(cvar_xckcl(icla)%p(iel)                            &
415                        + cvar_xchcl(icla)%p(iel)                            &
416                        + cvar_xnpcl(icla)%p(iel)*xmash(icla))*              &
417                          cpro_temp2(icla)%p(iel)
418
419          ! Prise en compte de l'humidite
420
421          if (ippmod(iccoal) .eq. 1) then
422            tfuel = tfuel +(cvar_xwtcl(icla)%p(iel))*cpro_temp2(icla)%p(iel)
423          endif
424        enddo
425
426        tfuel = tfuel/xmx2
427
428      else
429        tfuel = cpro_temp(iel)
430      endif
431
432      ! On recupere la valeur de Toxyd a partir de hoxyd
433
434      hoxyd = enthox(iel)
435
436      do ige = 1, ngazem
437        coefe(ige) = zero
438      enddo
439      coefe(io2) =( af3(io2)*f3m(iel)+af4(io2)*f4m(iel)    &
440                   +af5(io2)*f5m(iel)+af6(io2)*f6m(iel)    &
441                   +af7(io2)*f7m(iel)+af8(io2)*f8m(iel)    &
442                   +af9(io2)*f9m(iel)) * wmole(io2)
443
444      coefe(in2) =( af3(in2)*f3m(iel)+af4(in2)*f4m(iel)    &
445                   +af5(in2)*f5m(iel)+af6(in2)*f6m(iel)    &
446                   +af7(in2)*f7m(iel)+af8(in2)*f8m(iel)    &
447                   +af9(in2)*f9m(iel)) * wmole(in2)
448
449      coefe(ico2)=( af3(ico2)*f3m(iel)+af4(ico2)*f4m(iel)    &
450                   +af5(ico2)*f5m(iel)+af6(ico2)*f6m(iel)    &
451                   +af7(ico2)*f7m(iel)+af8(ico2)*f8m(iel)    &
452                   +af9(ico2)*f9m(iel)) * wmole(ico2)
453
454      coefe(ih2o)=( af3(ih2o)*f3m(iel)+af4(ih2o)*f4m(iel)    &
455                   +af5(ih2o)*f5m(iel)+af6(ih2o)*f6m(iel)    &
456                   +af7(ih2o)*f7m(iel)+af8(ih2o)*f8m(iel)    &
457                   +af9(ih2o)*f9m(iel)) * wmole(ih2o)
458
459      coefe(ico)=(  af3(ico)*f3m(iel)+af4(ico)*f4m(iel)      &
460                   +af5(ico)*f5m(iel)+af6(ico)*f6m(iel)      &
461                   +af7(ico)*f7m(iel)+af8(ico)*f8m(iel)      &
462                   +af9(ico)*f9m(iel)) * wmole(ico)
463
464      coefe(ihy)=(  af3(ihy)*f3m(iel)+af4(ihy)*f4m(iel)      &
465                   +af5(ihy)*f5m(iel)+af6(ihy)*f6m(iel)      &
466                   +af7(ihy)*f7m(iel)+af8(ihy)*f8m(iel)      &
467                   +af9(ihy)*f9m(iel)) * wmole(ihy)
468
469      ! Nbre de mole qui a reagis avec CO pour faire du CO2
470
471      deltamol   = (coefe(io2)-yo2ox)/wmole(io2)
472      react      = min(2.d0*deltamol,coefe(ico)/wmole(ico))
473      coefe(io2) = coefe(io2)-react/2.d0*wmole(io2)
474      coefe(ico) = coefe(ico)-react     *wmole(ico)
475      coefe(ico2) = coefe(ico2)+react   *wmole(ico2)
476
477      som = coefe(io2) +coefe(in2)+coefe(ico2) &
478           +coefe(ih2o)+coefe(ico)+coefe(ihy)
479      coefe(io2) = coefe(io2) /som
480      coefe(in2) = coefe(in2) /som
481      coefe(ico2) = coefe(ico2) /som
482      coefe(ih2o) = coefe(ih2o) /som
483      coefe(ico) = coefe(ico) /som
484      coefe(ihy) = coefe(ihy) /som
485
486      do icha = 1, ncharm
487        f1mc(icha) = zero
488        f2mc(icha) = zero
489      enddo
490
491      mode = 1
492      call cs_coal_htconvers1(mode, hoxyd, coefe, f1mc, f2mc, toxyd)
493
494      toxmin = min(toxmin,toxyd)
495      toxmax = max(toxmax,toxyd)
496
497      if (toxyd .gt. cpro_temp(iel)) then
498        toxyd = cpro_temp(iel)
499      endif
500
501      ! On initialise par les temperatures Toxy et Tfuel aux extremites
502
503      dirac  =  dfuel(iel)*tfuel + doxyd(iel)*toxyd
504
505      ! On recupere la valeur de la temperature moyenne
506
507      tmpgaz = cpro_temp(iel)
508
509      bb1 = max(0.D0      ,pdfm1(iel))
510      bb2 = min(fs4no(iel),pdfm2(iel))
511      bb3 = max(fs4no(iel),pdfm1(iel))
512      bb4 = min(1.d0      ,pdfm2(iel))
513
514      if (bb2 .gt. bb1) then
515        lro = hrec(iel)
516      else
517        lro = 0.d0
518      endif
519      if (bb4 .gt. bb3) then
520        lrf = hrec(iel)
521      else
522        lrf = 0.d0
523      endif
524
525      qqq = bb2**2 - bb1**2
526      rrr = bb2    - bb1
527      sss = bb4**2 - bb3**2
528      ttt = bb4    - bb3
529      uuu = 1.d0   - fs4no(iel)
530
531      gt1 = lro*qqq/(2.d0*fs4no(iel))
532      gt2 = lro*rrr
533
534      gt10= lrf*sss/(2.d0*uuu)
535      gt20= lrf*ttt
536      gt30= lrf*ttt/uuu
537
538      ts4num = tmpgaz - dirac + toxyd*(gt1 -gt2)        &
539                              - tfuel*(gt10+gt20-gt30)
540      ts4den = gt1-gt10+gt30
541
542      ts4 = ts4num/ts4den
543
544      ! Calcul de l'enthalpie du charbon a Tfuel
545
546     xxf = 0.d0
547     hhf = 0.d0
548
549     do numcha=1,ncharb
550
551       ! H(mv1, TFUEL)
552
553        do ige = 1, ngazem
554          coefe(ige) = zero
555        enddo
556        den = a1(numcha)*wmole(ichx1c(numcha))+b1(numcha)*wmole(ico)    &
557             +c1(numcha)*wmole(ih2o)          +d1(numcha)*wmole(ih2s)   &
558             +e1(numcha)*wmole(ihcn)          +f1(numcha)*wmole(inh3)
559        coefe(ichx1) = a1(numcha)*wmole(ichx1c(numcha)) / den
560        coefe(ico ) = b1(numcha)*wmole(ico)            / den
561        coefe(ih2o) = c1(numcha)*wmole(ih2o)           / den
562        coefe(ih2s) = d1(numcha)*wmole(ih2s)           / den
563        coefe(ihcn) = e1(numcha)*wmole(ihcn)           / den
564        coefe(inh3) = f1(numcha)*wmole(inh3)           / den
565
566        do icha = 1, ncharm
567          f1mc(icha) = zero
568          f2mc(icha) = zero
569        enddo
570        f1mc(numcha) = 1.d0
571
572        mode      = -1
573        call cs_coal_htconvers1(mode, xhf1, coefe, f1mc, f2mc, tfuel)
574
575        ! H(mv2, TFUEL)
576
577        do ige = 1, ngazem
578          coefe(ige) = zero
579        enddo
580        den = a2(numcha)*wmole(ichx2c(numcha)) +b2(numcha)*wmole(ico)   &
581             +c2(numcha)*wmole(ih2o)           +d2(numcha)*wmole(ih2s)  &
582             +e2(numcha)*wmole(ihcn)           +f2(numcha)*wmole(inh3)
583        coefe(ichx2) = a2(numcha)*wmole(ichx2c(numcha)) /den
584        coefe(ico ) = b2(numcha)*wmole(ico)            /den
585        coefe(ih2o) = c2(numcha)*wmole(ih2o)           /den
586        coefe(ih2s) = d2(numcha)*wmole(ih2s)           /den
587        coefe(ihcn) = e2(numcha)*wmole(ihcn)           /den
588        coefe(inh3) = f2(numcha)*wmole(ihcn)           /den
589
590        do icha = 1, ncharm
591          f1mc(icha) = zero
592          f2mc(icha) = zero
593        enddo
594        f2mc(numcha) = 1.d0
595
596        mode      = -1
597        call cs_coal_htconvers1(mode, xhf2, coefe, f1mc, f2mc, tfuel)
598
599        xxf = xxf + cvar_f1m(numcha)%p(iel)                   &
600                  + cvar_f2m(numcha)%p(iel)
601        hhf = hhf + cvar_f1m(numcha)%p(iel)*xhf1              &
602                  + cvar_f2m(numcha)%p(iel)*xhf2
603
604      enddo
605
606      if (xxf .gt. epsicp) then
607
608        hfuel = hhf / xxf
609
610        hfs4ad = fs4no(iel)*hfuel+(1.d0-fs4no(iel))*hoxyd
611
612        do ige = 1, ngazem
613          coefe(ige) = zero
614        enddo
615        do ige = 1, ngazg
616          coefe(ige) = yfs4no(iel,ige)
617        enddo
618
619        do icha = 1, ncharm
620          f1mc(icha) = zero
621          f2mc(icha) = zero
622        enddo
623
624        mode = 1
625        call cs_coal_htconvers1(mode, hfs4ad, coefe, f1mc, f2mc, tfs4ad)
626
627        ! Calcul pour affichage
628
629        ts4min = min(ts4min,ts4)
630        ts4max = max(ts4max,ts4)
631
632        ts4admin= min(ts4admin,tfs4ad)
633        ts4admax= max(ts4admax,tfs4ad)
634
635        somm = 0.d0
636        do ige = 1, ngazg
637          somm = somm + yfs4no(iel,ige)
638        enddo
639        sommin=min(sommin,somm)
640        sommax=max(sommax,somm)
641
642        if ((ts4 .gt. min(toxyd,tmpgaz,tfuel)) .and. &
643            (ts4 .lt. 2400.d0)) inok = inok + 1
644        if (ts4 .lt. min(toxyd,tmpgaz,tfuel)) then
645          if (ts4 .ge. 300.d0) then
646            i300  = i300  + 1
647          else if (ts4 .gt. 0) then
648            i000  = i000  + 1
649          else
650            imini = imini + 1
651          endif
652        endif
653
654        if (ts4 .gt. 2400.d0) then
655          if (ts4 .lt. 2500) then
656            i2500 = i2500 +1
657          else if (ts4 .lt. 2600) then
658            i2600 = i2600 +1
659          else if (ts4 .lt. 2700) then
660            i2700 = i2700 +1
661          else if (ts4 .lt. 2800) then
662            i2800 = i2800 +1
663          else if (ts4 .lt. 3000) then
664            i3000 = i3000 +1
665          else if (ts4 .lt. 3500) then
666            i3500 = i3500 +1
667          else if (ts4 .lt. 4000) then
668            i4000 = i4000 +1
669          else if (ts4 .lt. 5000) then
670            i5000 = i5000 +1
671          else
672            imaxi = imaxi + 1
673          endif
674        endif
675
676        ! Fin calcul pour affichage
677
678        ! Clipping de Ts4 : a min(toxyd,tmpgaz,tfuel) en min
679        !                   a ts4ad                   en max
680
681        nbpt = nbpt + 1
682        if (ts4 .lt. min(toxyd,tmpgaz,tfuel)) then
683           nbclip1 = nbclip1 + 1
684           ts4 = min(toxyd,tmpgaz,tfuel)
685        endif
686
687        if (ts4 .gt. tfs4ad) then
688           nbclip2 = nbclip2 + 1
689           ts4 = tfs4ad
690        endif
691
692        ! Concentration oxygene
693
694        xo2 = cpro_yo2(iel)             &
695             *cpro_mmel(iel)/wmole(io2)
696
697        ! Integration
698
699        do i = 1, npart+1
700          gs(i) = pdfm1(iel)+dble(i-1)/dble(npart)*(pdfm2(iel)-pdfm1(iel))
701          ! compute T
702          if(gs(i) .lt. fs4no(iel)) then
703            tt(i) = (ts4-toxyd)/fs4no(iel)* gs(i) + toxyd
704          else
705            tt(i) = (tfuel-ts4)/(1.d0-fs4no(iel))*gs(i)           &
706                   + tfuel - (tfuel-ts4)/(1.d0-fs4no(iel))
707          endif
708          ! compute yo2
709          if (gs(i) .lt. fs4no(iel))  then
710            yyo2(i) = (yo2s4-yo2ox)/fs4no(iel) * gs(i) + yo2ox
711          else if (gs(i) .lt. fs3no(iel)) then
712            aa = yo2s4/(fs4no(iel)-fs3no(iel))
713            yyo2(i) = aa * (gs(i) -fs3no(iel))
714          else
715            yyo2(i) = 0.d0
716          endif
717        enddo
718
719        ! no integration
720
721        dgs = (pdfm2(iel) - pdfm1(iel)) / dble(npart)
722
723        ! Compute K1*EXP(-E1/T)
724
725        if (ipdf1 .eq. 1) then
726
727          cpro_exp1(iel) =   kk1*exp(-ee1/toxyd)*doxyd(iel)  &
728                           + kk1*exp(-ee1/tfuel)*dfuel(iel)
729
730          do i = 1, npart+1
731            val(i) = kk1*exp(-ee1/tt(i))*hrec(iel)
732          enddo
733
734          do i = 1, npart
735            cpro_exp1(iel) =   cpro_exp1(iel)               &
736                             + 0.5d0*dgs*(val(i)+val(i+1))
737          enddo
738
739        endif
740
741        ! Compute K2*EXP(-E2/T)
742
743          if (ipdf2 .eq. 1) then
744
745            if (xo2 .gt. 0.d0) then
746
747              if (xo2.gt.0.018d0) then
748                bb=0.d0
749              else if (xo2 .lt. 0.0025d0) then
750                bb=1.d0
751              else
752                bb=(0.018d0-xo2)/(0.018d0-0.0025d0)
753              endif
754
755              cpro_exp2(iel) =   kk2*exp(-ee2/toxyd)*doxyd(iel)     &
756                                    *(xo2**bb)                      &
757                               + kk2*exp(-ee2/tfuel)*dfuel(iel)     &
758                                    *(xo2**bb)
759
760              do i = 1, npart+1
761                val(i) = kk2*exp(-ee2/tt(i))*hrec(iel)
762              enddo
763
764              do i = 1, npart
765                cpro_exp2(iel) =   cpro_exp2(iel)                   &
766                                 + 0.5d0*dgs*(val(i)+val(i+1))*(xo2**bb)
767              enddo
768            else
769              cpro_exp2(iel) = 0.d0
770            endif
771
772          endif
773
774          ! Compute K3*EXP(-E3/T)
775
776          if (ipdf3 .eq. 1) then
777
778            if (xo2 .gt. 0.d0) then
779
780              cpro_exp3(iel) =  kk3*exp(-ee3/toxyd)                &
781                                   *doxyd(iel)*(yo2ox**0.5d0)      &
782                               +kk3*exp(-ee3/tfuel)                &
783                                   *dfuel(iel)*(yo2cb**0.5d0)
784
785              do i = 1, npart+1
786                if (yyo2(i).gt.0.d0) then
787                  if (gs(i).le.fs3no(iel)) then
788                    val(i) = kk3*exp(-ee3/tt(i))*hrec(iel)*(yyo2(i)**0.5d0)
789                  else
790                    val(i) = 0.d0
791                  endif
792                else
793                  val(i) = 0.d0
794                endif
795              enddo
796
797              do i = 1, npart
798                cpro_exp3(iel) =   cpro_exp3(iel)                &
799                                 + 0.5d0*dgs*(val(i)+val(i+1))
800              enddo
801
802            else
803              cpro_exp3(iel) = 0.d0
804            endif
805          endif
806
807          ! Compute K4*EXP(-E4/T)
808
809          if (ipdf4 .eq. 1) then
810
811            if (xo2 .gt. 0.d0) then
812
813              if(xo2.gt.0.018d0) then
814                bb=0.d0
815              else if(xo2 .lt. 0.0025d0) then
816                bb=1.d0
817              else
818                bb=(0.018d0-xo2)/(0.018d0-0.0025d0)
819              endif
820
821              cpro_exp4(iel) =   kk4*exp(-ee4/toxyd)*doxyd(iel)     &
822                                    *(xo2**bb)                      &
823                               + kk4*exp(-ee4/tfuel)*dfuel(iel)     &
824                                    *(xo2**bb)
825
826              do i = 1, npart+1
827                val(i) = kk4*exp(-ee4/tt(i))*hrec(iel)
828              enddo
829
830              do i = 1, npart
831                cpro_exp4(iel) =   cpro_exp4(iel)                &
832                                 + 0.5d0*dgs*(val(i)+val(i+1))*(xo2**bb)
833              enddo
834            else
835              cpro_exp4(iel) = 0.d0
836            endif
837
838          endif
839
840        ! Compute K5*EXP(-E5/T)
841
842        if (ipdf5 .eq. 1) then
843
844          cpro_exp5(iel) =   kk5*exp(-ee5/toxyd)*doxyd(iel)   &
845                           + kk5*exp(-ee5/tfuel)*dfuel(iel)
846
847          do i = 1, npart+1
848            val(i) = kk5*exp(-ee5/tt(i))*hrec(iel)
849          enddo
850
851          do i = 1, npart
852            cpro_exp5(iel) =   cpro_exp5(iel)                 &
853                             + 0.5d0*dgs*(val(i)+val(i+1))
854          enddo
855
856        endif
857
858      endif
859
860    endif
861  enddo
862
863  deallocate(cvar_xckcl, cvar_xchcl)
864  deallocate(cvar_xnpcl, cvar_xwtcl)
865  deallocate(cvar_f1m, cvar_f2m)
866  deallocate(cpro_temp2)
867
868  log_active = cs_log_default_is_active()
869  if (log_active .eqv. .true.) then
870
871    if (irangp .ge. 0) then
872      call parcpt(inok)
873      call parcpt(i300)
874      call parcpt(i000)
875      call parcpt(imini)
876      call parcpt(i2500)
877      call parcpt(i2600)
878      call parcpt(i2700)
879      call parcpt(i2800)
880      call parcpt(i3000)
881      call parcpt(i3500)
882      call parcpt(i4000)
883      call parcpt(i5000)
884      call parcpt(imaxi)
885      call parcpt(nbpt)
886      call parcpt(nbclip1)
887      call parcpt(nbclip2)
888      call parmin(ts4min)
889      call parmax(ts4max)
890      call parmin(ts4admin)
891      call parmax(ts4admax)
892      call parmin(sommin)
893      call parmax(sommax)
894
895      call parcpt(nbclip30)
896      call parcpt(nbclip31)
897      call parmin(yo2min)
898      call parmax(yo2max)
899      call parmin(yo2min1)
900      call parmax(yo2max1)
901
902      call parmin(toxmin)
903      call parmax(toxmax)
904      call parmin(yo2oxmin)
905      call parmax(yo2oxmax)
906    endif
907
908    !---------------------------------------------------------------------------
909
910    write(nfecra,*) ' '
911    write(nfecra,*) ' Min max of TSox ', toxmin, toxmax
912    write(nfecra,*) ' Min max of TS4 ', ts4min, ts4max
913    write(nfecra,*) '    number of points at Tmini < ts4 < 2400  ',inok
914    write(nfecra,*) '    number of points at         ts4 < 0     ',imini
915    write(nfecra,*) '    number of points at 0     < ts4 < 300   ',i000
916    write(nfecra,*) '    number of points at 300   < ts4 < Tmini ',i300
917    write(nfecra,*) '    number of points at 2400  < ts4 < 2500  ',i2500
918    write(nfecra,*) '    number of points at 2500  < ts4 < 2600  ',i2600
919    write(nfecra,*) '    number of points at 2600  < ts4 < 2700  ',i2700
920    write(nfecra,*) '    number of points at 2700  < ts4 < 2800  ',i2800
921    write(nfecra,*) '    number of points at 2800  < ts4 < 3000  ',i3000
922    write(nfecra,*) '    number of points at 3000  < ts4 < 3500  ',i3500
923    write(nfecra,*) '    number of points at 3500  < ts4 < 4000  ',i4000
924    write(nfecra,*) '    number of points at 4000  < ts4 < 5000  ',i5000
925    write(nfecra,*) '    number of points at 5000  < ts4         ',imaxi
926    write(nfecra,*) ' Min max of TS4ad ',ts4admin,ts4admax
927    write(nfecra,*) ' Min max concentration in fs4 ', sommin, sommax
928    write(nfecra,*) ' Clip to min at ', nbclip1,' of ', nbpt,' points'
929    write(nfecra,*) ' Clip to max at ', nbclip2,' of ', nbpt,' points'
930
931    write(nfecra,*) ' '
932    write(nfecra,*) ' Min max of Yo2ox at 0 ',yo2oxmin, yo2oxmax
933    write(nfecra,*) ' Min max of Yo2 in fs4 before clipping ', yo2min, yo2max
934    write(nfecra,*) ' Clip to min over Yo2 in fs4       ', nbclip30
935    write(nfecra,*) ' Clip to max over Yo2 in fs4       ', nbclip31
936    write(nfecra,*) ' Min max of Yo2 in fs4 after clipping ', yo2min1, yo2max1
937
938    !---------------------------------------------------------------------------
939
940  endif
941
942endif
943
944!----
945! End
946!----
947
948return
949end subroutine
950