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_fuel_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_fuel_incl
74use cs_c_bindings
75
76!===============================================================================
77
78implicit none
79
80! Arguments
81
82integer          ncel
83integer          indpdf(ncel)
84
85double precision pdfm1(ncel) , pdfm2(ncel) , dfuel(ncel)
86double precision f3m(ncel)   , f4m(ncel)   , f5m(ncel)
87double precision f6m(ncel)   , f7m(ncel)   , f8m(ncel) ,f9m(ncel)
88double precision doxyd(ncel) , hrec(ncel)
89double precision fs3no(ncel) , fs4no(ncel) , yfs4no(ncel,ngazg)
90double precision enthox(ncel)
91!
92! VARIABLES LOCALES
93!
94integer          iel , icla , i
95integer          ipdf1 , ipdf2 , ipdf3
96integer          npart
97parameter        (npart = 200)
98
99double precision ee1,ee2,ee3,kk1,kk2,kk3
100double precision wmo2,tg,xo2,bb,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,mode
110double precision hhf , hfuel,tfs4ad,hfs4ad
111double precision aa
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
126double precision, dimension(:), pointer :: cpro_exp1, cpro_exp2, cpro_exp3
127double precision, dimension(:), pointer :: cpro_temp, cpro_yo2, cpro_mmel
128type(pmapper_double_r1), dimension(:), allocatable :: cvar_yfolcl, cvar_h2cl
129type(pmapper_double_r1), dimension(:), allocatable :: cpro_temp2
130
131logical(kind=c_bool) :: log_active
132
133!===============================================================================
134! 1. Preliminary computations
135!===============================================================================
136
137! Molar masses
138!
139wmo2 = wmole(io2)
140!
141! pointers
142
143call field_get_val_s(ighcn1,cpro_exp1)
144call field_get_val_s(ighcn2,cpro_exp2)
145call field_get_val_s(ignoth,cpro_exp3)
146
147! Parameters of Arrhenius laws
148
149kk1 = 3.0e12
150ee1 = 3.0e4
151kk2 = 1.2e10
152ee2 = 3.35e4
153kk3 = 3.4e12
154ee3 = 6.69e4
155
156! Pour les termes, indicateur de calcul par PDF ou non
157!       = 1 --> passage par pdf
158!       = 0 --> on ne passe pas par les pdf
159
160ipdf1 = 0
161ipdf2 = 0
162ipdf3 = 1
163
164! Initialisation
165
166do iel = 1, ncel
167  cpro_exp1(iel)  = zero
168  cpro_exp2(iel)  = zero
169  cpro_exp3(iel)  = zero
170enddo
171
172!===============================================================================
173! 2. CALCUL SANS LES PDF
174!===============================================================================
175
176call field_get_val_s(itemp, cpro_temp)
177call field_get_val_s(iym1(io2), cpro_yo2)
178call field_get_val_s(immel, cpro_mmel)
179
180do iel = 1, ncel
181
182  tg  = cpro_temp(iel)
183  yo2 = cpro_yo2(iel)
184  xo2 = yo2*cpro_mmel(iel)/wmo2
185
186  ! Reaction HCN + NO + 1/4 O2 ---> N2 + 1/2 H2O + CO
187
188  cpro_exp1(iel)  = kk1*exp(-ee1/tg)
189
190  ! Reaction HCN + 5/4 O2 --> NO + 1/2 H2O  + CO
191
192  if (xo2 .gt. 0.018d0) then
193    bb = 0.d0
194  else if (xo2 .lt. 0.0025d0) then
195    bb= 1.d0
196  else
197    bb=(0.018d0-xo2)/(0.018d0-0.0025d0)
198  endif
199  cpro_exp2(iel)  = kk2*exp(-ee2/tg)*(xo2**bb)
200
201  ! No thermique : Zeldovich
202
203  cpro_exp3(iel)  = kk3*exp(-ee3/tg)*(xo2**0.5d0)
204
205enddo
206
207!===============================================================================
208! 3. Computation with the PDFs
209!===============================================================================
210!
211if (ipdf1.eq.1 .or. ipdf2.eq.1 .or. ipdf3.eq.1) then
212
213  ! Arrays of pointers containing the field values for each class
214  ! (loop on cells outside loop on classes)
215  allocate(cvar_yfolcl(nclafu), cvar_h2cl(nclafu))
216  allocate(cpro_temp2(nclafu))
217  do icla = 1, nclafu
218    call field_get_val_s(ivarfl(isca(iyfol(icla))), cvar_yfolcl(icla)%p)
219    call field_get_val_s(ivarfl(isca(ih2(icla))), cvar_h2cl(icla)%p)
220    call field_get_val_s(itemp2(icla), cpro_temp2(icla)%p)
221  enddo
222
223  inok = 0
224  i300 = 0
225  i000 = 0
226  imini= 0
227  i2500= 0
228  i2600= 0
229  i2700= 0
230  i2800= 0
231  i3000= 0
232  i3500= 0
233  i4000= 0
234  i5000= 0
235  imaxi= 0
236  nbpt = 0
237  nbclip1 =0
238  nbclip2 = 0
239  ts4min= 1.d+20
240  ts4max=-1.d+20
241  sommin= 1.d+20
242  sommax=-1.d+20
243  ts4admin= 1.d+20
244  ts4admax=-1.d+20
245  toxmin = 1.d+20
246  toxmax =-1.d+20
247  yo2oxmin= 1.d+20
248  yo2oxmax=-1.d20
249
250  nbclip30 = 0
251  nbclip31 = 0
252  yo2min  = 1.d+20
253  yo2max  =-1.d+20
254  yo2min1 = 1.d+20
255  yo2max1 =-1.d+20
256
257  do iel=1,ncel
258
259    if (indpdf(iel) .eq. 1) then
260
261      ! Calcul de Yo2 dans l'oxydant
262      !           Yo2 en fs4
263
264      bb1 = max(0.d0      ,pdfm1(iel))
265      bb2 = min(fs3no(iel),pdfm2(iel))
266      if (bb2 .gt. bb1) then
267        lro = hrec(iel)
268      else
269        lro = 0.d0
270      endif
271      qqq = bb2**2 - bb1**2
272      rrr = bb2    - bb1
273      gt1 = lro*qqq/(2.d0*fs3no(iel))
274      gt2 = lro*rrr
275      gt3 = doxyd(iel)
276      yo2cb = 0.d0
277
278      if (cpro_yo2(iel) .gt. 0.d0) then
279        yo2ox = cpro_yo2(iel)/(-gt1+gt2+gt3)
280
281        yo2oxmin = min(yo2oxmin,yo2ox)
282        yo2oxmax = max(yo2oxmax,yo2ox)
283
284        yo2moy = cpro_yo2(iel)
285        dirac  =  dfuel(iel)*yo2cb + doxyd(iel)*yo2ox
286
287        bb1 = max(0.D0      ,pdfm1(iel))
288        bb2 = min(fs4no(iel),pdfm2(iel))
289        bb3 = max(fs4no(iel),pdfm1(iel))
290        bb4 = min(fs3no(iel),pdfm2(iel))
291        if (bb2 .gt. bb1) then
292          lro = hrec(iel)
293        else
294          lro = 0.d0
295        endif
296        if (bb4 .gt. bb3) then
297          lrf = hrec(iel)
298        else
299          lrf = 0.d0
300        endif
301
302        qqq = bb2**2 - bb1**2
303        rrr = bb2    - bb1
304        sss = bb4**2 - bb3**2
305        ttt = bb4    - bb3
306        uuu = fs4no(iel)-fs3no(iel)
307
308        gt1 = lro*qqq/(2.d0*fs4no(iel))
309        gt2 = lro*rrr
310
311        gt10= lrf*sss/(2.d0*uuu)
312        gt20= lrf*ttt*fs3no(iel)/uuu
313
314        yo24num = yo2moy - dirac + yo2ox*(gt1 -gt2)
315        yo24den = gt1+gt10-gt20
316
317        yo2s4  = yo24num/yo24den
318
319      else
320        yo2ox = 0.d0
321        yo2s4 = 0.d0
322      endif
323
324      ! Affichage et clipping
325
326      yo2min = min(yo2s4,yo2min)
327      yo2max = max(yo2s4,yo2max)
328      if (yo2s4 .lt. 0.d0) then
329         nbclip30 = nbclip30+1
330         yo2s4 = 0.d0
331      endif
332      if (yo2s4 .gt. yo2ox) then
333         nbclip31 = nbclip31+1
334         yo2s4 = yo2ox
335      endif
336      yo2min1 = min(yo2s4,yo2min1)
337      yo2max1 = max(yo2s4,yo2max1)
338
339      ! Calcul de Tfioul moyen
340      ! ----------------------
341
342      tfuel = 0.d0
343      xmx2  = 0.d0
344
345      do icla = 1, nclafu
346        xmx2 = xmx2 + cvar_yfolcl(icla)%p(iel)
347      enddo
348      if (xmx2 .gt. 0.d0) then
349        do icla=1,nclafu
350          tfuel = tfuel +cvar_yfolcl(icla)%p(iel)*cpro_temp2(icla)%p(iel)
351        enddo
352
353        tfuel = tfuel/xmx2
354
355      else
356        tfuel = cpro_temp(iel)
357      endif
358
359      ! On recupere la valeur de Toxyd a partir de hoxyd
360
361      hoxyd = enthox(iel)
362
363      do ige = 1, ngazem
364        coefe(ige) = zero
365      enddo
366      coefe(io2) =( af3(io2)*f3m(iel)+af4(io2)*f4m(iel)+af5(io2)*f5m(iel) &
367                   +af6(io2)*f6m(iel)+af7(io2)*f7m(iel)+af8(io2)*f8m(iel) &
368                   +af9(io2)*f9m(iel)) * wmole(io2)
369
370      coefe(in2) =( af3(in2)*f3m(iel)+af4(in2)*f4m(iel)+af5(in2)*f5m(iel) &
371                   +af6(in2)*f6m(iel)+af7(in2)*f7m(iel)+af8(in2)*f8m(iel) &
372                   +af9(in2)*f9m(iel)) * wmole(in2)
373
374      coefe(ico2)=( af3(ico2)*f3m(iel)+af4(ico2)*f4m(iel)+af5(ico2)*f5m(iel) &
375                   +af6(ico2)*f6m(iel)+af7(ico2)*f7m(iel)+af8(ico2)*f8m(iel) &
376                   +af9(ico2)*f9m(iel)) * wmole(ico2)
377
378      coefe(ih2o)=( af3(ih2o)*f3m(iel)+af4(ih2o)*f4m(iel)+af5(ih2o)*f5m(iel) &
379                   +af6(ih2o)*f6m(iel)+af7(ih2o)*f7m(iel)+af8(ih2o)*f8m(iel) &
380                   +af9(ih2o)*f9m(iel)) * wmole(ih2o)
381
382      coefe(ico)=(  af3(ico)*f3m(iel)+af4(ico)*f4m(iel)+af5(ico)*f5m(iel)    &
383                   +af6(ico)*f6m(iel)+af7(ico)*f7m(iel)+af8(ico)*f8m(iel)    &
384                   +af9(ico)*f9m(iel)) * wmole(ico)
385
386      coefe(ihy)=(  af3(ihy)*f3m(iel)+af4(ihy)*f4m(iel)+af5(ihy)*f5m(iel)    &
387                   +af6(ihy)*f6m(iel)+af7(ihy)*f7m(iel)+af8(ihy)*f8m(iel)    &
388                   +af9(ihy)*f9m(iel)) * wmole(ihy)
389
390      ! Nbre de mole qui a reagis avec CO pour faire du CO2
391
392      deltamol   = (coefe(io2)-yo2ox)/wmole(io2)
393      react      = min(2.d0*deltamol,coefe(ico)/wmole(ico))
394      coefe(io2) = coefe(io2)-react/2.d0*wmole(io2)
395      coefe(ico) = coefe(ico)-react     *wmole(ico)
396      coefe(ico2) = coefe(ico2)+react   *wmole(ico2)
397
398      som = coefe(io2) +coefe(in2)+coefe(ico2) &
399           +coefe(ih2o)+coefe(ico)+coefe(ihy)
400      coefe(io2) = coefe(io2) /som
401      coefe(in2) = coefe(in2) /som
402      coefe(ico2) = coefe(ico2) /som
403      coefe(ih2o) = coefe(ih2o) /som
404      coefe(ico) = coefe(ico) /som
405      coefe(ihy) = coefe(ihy) /som
406
407      do icha = 1, ncharm
408        f1mc(icha) = zero
409        f2mc(icha) = zero
410      enddo
411
412      mode = 1
413      call cs_fuel_htconvers1(mode, hoxyd, coefe, toxyd)
414
415      toxmin = min(toxmin,toxyd)
416      toxmax = max(toxmax,toxyd)
417
418      if (toxyd .gt. cpro_temp(iel)) then
419        toxyd = cpro_temp(iel)
420      endif
421
422      ! On initialise par les temperatures Toxy et Tfuel aux extremites
423
424      dirac  =  dfuel(iel)*tfuel + doxyd(iel)*toxyd
425
426      ! On recupere la valeur de la temperature moyenne
427
428      tmpgaz = cpro_temp(iel)
429
430      bb1 = max(0.D0      ,pdfm1(iel))
431      bb2 = min(fs4no(iel),pdfm2(iel))
432      bb3 = max(fs4no(iel),pdfm1(iel))
433      bb4 = min(1.d0      ,pdfm2(iel))
434
435      if (bb2 .gt. bb1) then
436        lro = hrec(iel)
437      else
438        lro = 0.d0
439      endif
440      if (bb4 .gt. bb3) then
441        lrf = hrec(iel)
442      else
443        lrf = 0.d0
444      endif
445
446      qqq = bb2**2 - bb1**2
447      rrr = bb2    - bb1
448      sss = bb4**2 - bb3**2
449      ttt = bb4    - bb3
450      uuu = 1.d0   - fs4no(iel)
451
452      gt1 = lro*qqq/(2.d0*fs4no(iel))
453      gt2 = lro*rrr
454
455      gt10= lrf*sss/(2.d0*uuu)
456      gt20= lrf*ttt
457      gt30= lrf*ttt/uuu
458
459      ts4num = tmpgaz - dirac + toxyd*(gt1 -gt2)        &
460                              - tfuel*(gt10+gt20-gt30)
461      ts4den = gt1-gt10+gt30
462
463      ts4 = ts4num/ts4den
464
465      ! Calcul de l'enthalpie des vapeurs a Tfuel
466
467      hhf = 0.D0
468      do icla=1,nclafu
469        hhf = hhf + cvar_h2cl(icla)%p(iel)
470      enddo
471
472      if (xmx2 .gt. epsifl) then
473
474        hfuel = hhf / xmx2
475
476        hfs4ad = fs4no(iel)*hfuel+(1.d0-fs4no(iel))*hoxyd
477
478        do ige = 1, ngazem
479          coefe(ige) = zero
480        enddo
481        do ige = 1, ngazg
482          coefe(ige) = yfs4no(iel,ige)
483        enddo
484
485        mode = 1
486        call cs_fuel_htconvers1(mode, hfs4ad, coefe, tfs4ad)
487
488        ! Calcul pour affichage
489
490        ts4min = min(ts4min,ts4)
491        ts4max = max(ts4max,ts4)
492
493        ts4admin= min(ts4admin,tfs4ad)
494        ts4admax= max(ts4admax,tfs4ad)
495
496        somm = 0.d0
497        do ige = 1, ngazg
498          somm = somm + yfs4no(iel,ige)
499        enddo
500        sommin=min(sommin,somm)
501        sommax=max(sommax,somm)
502
503        if ((ts4 .gt. min(toxyd,tmpgaz,tfuel)) .and. &
504            (ts4 .lt. 2400.d0)) inok = inok + 1
505        if (ts4 .lt. min(toxyd,tmpgaz,tfuel)) then
506          if (ts4 .ge. 300.d0) then
507            i300  = i300  + 1
508          else if (ts4 .gt. 0) then
509            i000  = i000  + 1
510          else
511            imini = imini + 1
512          endif
513        endif
514
515        if (ts4 .gt. 2400.d0) then
516          if (ts4 .lt. 2500) then
517            i2500 = i2500 +1
518          else if (ts4 .lt. 2600) then
519            i2600 = i2600 +1
520          else if (ts4 .lt. 2700) then
521            i2700 = i2700 +1
522          else if (ts4 .lt. 2800) then
523            i2800 = i2800 +1
524          else if (ts4 .lt. 3000) then
525            i3000 = i3000 +1
526          else if (ts4 .lt. 3500) then
527            i3500 = i3500 +1
528          else if (ts4 .lt. 4000) then
529            i4000 = i4000 +1
530          else if (ts4 .lt. 5000) then
531            i5000 = i5000 +1
532          else
533            imaxi = imaxi + 1
534          endif
535        endif
536
537        ! Fin calcul pour affichage
538
539        ! Clipping de Ts4 : a min(toxyd,tmpgaz,tfuel) en min
540        !                   a ts4ad                   en max
541
542        nbpt = nbpt + 1
543        if (ts4 .lt. min(toxyd,tmpgaz,tfuel)) then
544           nbclip1 = nbclip1 + 1
545           ts4 = min(toxyd,tmpgaz,tfuel)
546        endif
547
548        if (ts4 .gt. tfs4ad) then
549           nbclip2 = nbclip2 + 1
550           ts4 = tfs4ad
551        endif
552
553        ! Concentration oxygene
554
555        xo2 = cpro_yo2(iel)             &
556             *cpro_mmel(iel)/wmole(io2)
557
558        ! Integration
559
560        do i = 1, npart+1
561          gs(i) = pdfm1(iel)+dble(i-1)/dble(npart)*(pdfm2(iel)-pdfm1(iel))
562          ! compute T
563          if(gs(i) .lt. fs4no(iel)) then
564            tt(i) = (ts4-toxyd)/fs4no(iel)* gs(i) + toxyd
565          else
566            tt(i) = (tfuel-ts4)/(1.d0-fs4no(iel))*gs(i)           &
567                   + tfuel - (tfuel-ts4)/(1.d0-fs4no(iel))
568          endif
569          ! compute yo2
570          if (gs(i) .lt. fs4no(iel))  then
571            yyo2(i) = (yo2s4-yo2ox)/fs4no(iel) * gs(i) + yo2ox
572          else if (gs(i) .lt. fs3no(iel)) then
573            aa = yo2s4/(fs4no(iel)-fs3no(iel))
574            yyo2(i) = aa * (gs(i) -fs3no(iel))
575          else
576            yyo2(i) = 0.d0
577          endif
578        enddo
579
580        ! no integration
581
582        dgs = (pdfm2(iel) - pdfm1(iel)) / dble(npart)
583
584        ! Compute K1*EXP(-E1/T)
585
586        if (ipdf1 .eq. 1) then
587
588          cpro_exp1(iel) =   kk1*exp(-ee1/toxyd)*doxyd(iel)  &
589                           + kk1*exp(-ee1/tfuel)*dfuel(iel)
590
591          do i = 1, npart+1
592            val(i) = kk1*exp(-ee1/tt(i))*hrec(iel)
593          enddo
594
595          do i = 1, npart
596            cpro_exp1(iel) =   cpro_exp1(iel)               &
597                             + 0.5d0*dgs*(val(i)+val(i+1))
598          enddo
599
600        endif
601
602        ! Compute K2*EXP(-E2/T)
603
604          if (ipdf2 .eq. 1) then
605
606            if (xo2 .gt. 0.d0) then
607
608              if (xo2.gt.0.018d0) then
609                bb=0.d0
610              else if (xo2 .lt. 0.0025d0) then
611                bb=1.d0
612              else
613                bb=(0.018d0-xo2)/(0.018d0-0.0025d0)
614              endif
615
616              cpro_exp2(iel) =   kk2*exp(-ee2/toxyd)*doxyd(iel)*(xo2**bb)  &
617                               + kk2*exp(-ee2/tfuel)*dfuel(iel)*(xo2**bb)
618
619              do i = 1, npart+1
620                val(i) = kk2*exp(-ee2/tt(i))*hrec(iel)
621              enddo
622
623              do i = 1, npart
624                cpro_exp2(iel) =   cpro_exp2(iel)                   &
625                                 + 0.5d0*dgs*(val(i)+val(i+1))*(xo2**bb)
626              enddo
627            else
628              cpro_exp2(iel) = 0.d0
629            endif
630
631          endif
632
633          ! Compute K3*EXP(-E3/T)
634
635          if (ipdf3 .eq. 1) then
636
637            if (xo2 .gt. 0.d0) then
638
639              cpro_exp3(iel) =  kk3*exp(-ee3/toxyd)                &
640                                   *doxyd(iel)*(yo2ox**0.5d0)      &
641                               +kk3*exp(-ee3/tfuel)                &
642                                   *dfuel(iel)*(yo2cb**0.5d0)
643
644              do i = 1, npart+1
645                if (yyo2(i).gt.0.d0) then
646                  if (gs(i).le.fs3no(iel)) then
647                    val(i) = kk3*exp(-ee3/tt(i))*hrec(iel)*(yyo2(i)**0.5d0)
648                  else
649                    val(i) = 0.d0
650                  endif
651                else
652                  val(i) = 0.d0
653                endif
654              enddo
655
656              do i = 1, npart
657                cpro_exp3(iel) =   cpro_exp3(iel)                &
658                                 + 0.5d0*dgs*(val(i)+val(i+1))
659              enddo
660
661            else
662              cpro_exp3(iel) = 0.d0
663            endif
664          endif
665
666      endif
667
668    endif
669  enddo
670
671  deallocate(cvar_yfolcl, cvar_h2cl)
672  deallocate(cpro_temp2)
673
674  log_active = cs_log_default_is_active()
675  if (log_active .eqv. .true.) then
676
677    if (irangp .ge. 0) then
678      call parcpt(inok)
679      call parcpt(i300)
680      call parcpt(i000)
681      call parcpt(imini)
682      call parcpt(i2500)
683      call parcpt(i2600)
684      call parcpt(i2700)
685      call parcpt(i2800)
686      call parcpt(i3000)
687      call parcpt(i3500)
688      call parcpt(i4000)
689      call parcpt(i5000)
690      call parcpt(imaxi)
691      call parcpt(nbpt)
692      call parcpt(nbclip1)
693      call parcpt(nbclip2)
694      call parmin(ts4min)
695      call parmax(ts4max)
696      call parmin(ts4admin)
697      call parmax(ts4admax)
698      call parmin(sommin)
699      call parmax(sommax)
700
701      call parcpt(nbclip30)
702      call parcpt(nbclip31)
703      call parmin(yo2min)
704      call parmax(yo2max)
705      call parmin(yo2min1)
706      call parmax(yo2max1)
707
708      call parmin(toxmin)
709      call parmax(toxmax)
710      call parmin(yo2oxmin)
711      call parmax(yo2oxmax)
712    endif
713
714    !---------------------------------------------------------------------------
715
716    write(nfecra,*) ' '
717    write(nfecra,*) ' Min max of TSox ', toxmin, toxmax
718    write(nfecra,*) ' Min max of TS4 ', ts4min, ts4max
719    write(nfecra,*) '    number of points at Tmini < ts4 < 2400  ',inok
720    write(nfecra,*) '    number of points at         ts4 < 0     ',imini
721    write(nfecra,*) '    number of points at 0     < ts4 < 300   ',i000
722    write(nfecra,*) '    number of points at 300   < ts4 < Tmini ',i300
723    write(nfecra,*) '    number of points at 2400  < ts4 < 2500  ',i2500
724    write(nfecra,*) '    number of points at 2500  < ts4 < 2600  ',i2600
725    write(nfecra,*) '    number of points at 2600  < ts4 < 2700  ',i2700
726    write(nfecra,*) '    number of points at 2700  < ts4 < 2800  ',i2800
727    write(nfecra,*) '    number of points at 2800  < ts4 < 3000  ',i3000
728    write(nfecra,*) '    number of points at 3000  < ts4 < 3500  ',i3500
729    write(nfecra,*) '    number of points at 3500  < ts4 < 4000  ',i4000
730    write(nfecra,*) '    number of points at 4000  < ts4 < 5000  ',i5000
731    write(nfecra,*) '    number of points at 5000  < ts4         ',imaxi
732    write(nfecra,*) ' Min max of TS4ad ',ts4admin,ts4admax
733    write(nfecra,*) ' Min max concentration in fs4 ', sommin, sommax
734    write(nfecra,*) ' Clip to min at ', nbclip1,' of ', nbpt,' points'
735    write(nfecra,*) ' Clip to max at ', nbclip2,' of ', nbpt,' points'
736
737    write(nfecra,*) ' '
738    write(nfecra,*) ' Min max of Yo2ox at 0 ',yo2oxmin, yo2oxmax
739    write(nfecra,*) ' Min max of Yo2 in fs4 before clipping ', yo2min, yo2max
740    write(nfecra,*) ' Clip to min over Yo2 in fs4       ', nbclip30
741    write(nfecra,*) ' Clip to max over Yo2 in fs4       ', nbclip31
742    write(nfecra,*) ' Min max of Yo2 in fs4 after clipping ', yo2min1, yo2max1
743
744    !---------------------------------------------------------------------------
745
746  endif
747
748endif
749
750!----
751! End
752!----
753
754return
755end subroutine
756