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