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 23!=============================================================================== 24! Function: 25! --------- 26 27!> \file covofv.f90 28!> 29!> \brief This subroutine performs the solving the convection/diffusion 30!> equation (with eventually source terms and/or drift) for a vectorial quantity 31!> over a time step. 32!> 33!------------------------------------------------------------------------------- 34 35!------------------------------------------------------------------------------- 36! Arguments 37!______________________________________________________________________________. 38! mode name role ! 39!______________________________________________________________________________! 40!> \param[in] nvar total number of variables 41!> \param[in] nscal total number of scalars 42!> \param[in] ncepdp number of cells with head loss 43!> \param[in] ncesmp number of cells with mass source term 44!> \param[in] iterns Navier-Stokes iteration number 45!> \param[in] iscal scalar number 46!> \param[in] icepdc index of cells with head loss 47!> \param[in] icetsm index of cells with mass source term 48!> \param[in] itypsm type of mass source term for the variables 49!> \param[in] dt time step (per cell) 50!> \param[in] ckupdc work array for the head loss 51!> \param[in] smacel variable value associated to the mass source 52!> term (for ivar=ipr, smacel is the mass flux 53!> \f$ \Gamma^n \f$) 54!> \param[in] viscf visc*surface/dist at internal faces 55!> \param[in] viscb visc*surface/dist at boundary faces 56!_______________________________________________________________________________ 57 58subroutine covofv & 59 ( nvar , nscal , ncepdp , ncesmp , & 60 iterns , iscal , & 61 icepdc , icetsm , & 62 itypsm , & 63 dt , & 64 ckupdc , smacel , & 65 viscf , viscb ) 66 67!=============================================================================== 68 69!=============================================================================== 70! Module files 71!=============================================================================== 72 73use, intrinsic :: iso_c_binding 74 75use paramx 76use numvar 77use entsor 78use optcal 79use cstphy 80use cstnum 81use ppppar 82use ppthch 83use coincl 84use cpincl 85use cs_fuel_incl 86use ppincl 87use ppcpfu 88use lagran 89use radiat 90use field 91use field_operator 92use mesh 93use parall 94use period 95use cs_f_interfaces 96use atchem 97use darcy_module 98use cs_c_bindings 99 100!=============================================================================== 101 102implicit none 103 104! Arguments 105 106integer nvar , nscal 107integer ncepdp , ncesmp 108integer iterns , iscal 109 110integer icepdc(ncepdp) 111integer icetsm(ncesmp), itypsm(ncesmp,nvar) 112 113double precision dt(ncelet) 114double precision ckupdc(6,ncepdp), smacel(ncesmp,nvar) 115double precision viscf(nfac), viscb(nfabor) 116 117! Local variables 118 119logical lprev 120character(len=80) :: chaine 121integer ivar 122integer ifac , iel, isou, jsou 123integer iiun, ibcl 124integer ivarsc 125integer iiscav 126integer ifcvsl, iflmas, iflmab 127integer iwarnp 128integer iconvp, imvisp 129integer iescap, ivissv 130integer idftnp 131integer iflid , st_prv_id, st_id, keydri, iscdri 132integer icvflb, f_dim, iflwgr 133integer key_buoyant_id, is_buoyant_fld 134 135integer ivoid(1) 136 137double precision sclnor 138double precision thetv , thets , thetap, thetp1 139double precision smbexp(3) 140double precision temp, idifftp 141double precision turb_schmidt, visls_0 142 143double precision rvoid(1) 144 145double precision, allocatable, dimension(:) :: w1 146double precision, allocatable, dimension(:,:) :: smbrv, gavinj 147double precision, allocatable, dimension(:,:,:) :: fimp 148double precision, allocatable, dimension(:,:) :: viscce 149double precision, allocatable, dimension(:,:) :: weighf 150double precision, allocatable, dimension(:) :: weighb 151 152double precision, dimension(:,:), pointer :: visten 153double precision, dimension(:,:), pointer :: cpro_wgrec_v 154double precision, dimension(:), pointer :: cpro_wgrec_s 155double precision, dimension(:), pointer :: imasfl, bmasfl 156double precision, dimension(:), pointer :: crom, croma, pcrom 157double precision, dimension(:,:), pointer :: coefap, cofafp 158double precision, dimension(:,:,:), pointer :: coefbp, cofbfp 159double precision, dimension(:), pointer :: visct 160double precision, dimension(:,:), pointer :: cpro_vect_st, cproa_vect_st 161double precision, dimension(:), pointer :: cpro_viscls 162! Darcy arrays 163double precision, allocatable, dimension(:) :: diverg, divflu 164double precision, dimension(:), pointer :: cpro_delay, cpro_sat 165double precision, dimension(:), pointer :: cproa_delay, cproa_sat 166double precision, dimension(:,:), pointer :: cvar_var, cvara_var 167 168type(var_cal_opt) :: vcopt 169type(var_cal_opt), target :: vcopt_loc 170type(var_cal_opt), pointer :: p_k_value 171type(c_ptr) :: c_k_value 172 173type(gwf_soilwater_partition) :: sorption_scal 174 175!=============================================================================== 176 177!=============================================================================== 178! 1. Initialization 179!=============================================================================== 180 181! --- Variable number 182ivar = isca(iscal) 183! Index of the field 184iflid = ivarfl(ivar) 185 186call field_get_val_v(iflid, cvar_var) 187call field_get_val_prev_v(iflid, cvara_var) 188 189! Key id for buoyant field (inside the Navier Stokes loop) 190call field_get_key_id("is_buoyant", key_buoyant_id) 191call field_get_key_int(iflid, key_buoyant_id, is_buoyant_fld) 192 193! If the vector is buoyant, it is inside the Navier Stokes loop, and so iterns >=1 194! otherwise it is outside of the loop and iterns = -1. 195if ( (is_buoyant_fld.eq. 1 .and. iterns.eq.-1) & 196 .or.(is_buoyant_fld.eq. 0 .and. iterns.ne.-1)) return 197 198! Key id for drift scalar 199call field_get_key_id("drift_scalar_model", keydri) 200 201! Id of the mass flux 202call field_get_key_int(iflid, kimasf, iflmas) ! interior mass flux 203! Pointer to the internal mass flux 204call field_get_val_s(iflmas, imasfl) 205 206! Id of the mass flux 207call field_get_key_int(iflid, kbmasf, iflmab) ! boundary mass flux 208! Pointer to the Boundary mass flux 209call field_get_val_s(iflmab, bmasfl) 210 211call field_get_key_struct_var_cal_opt(iflid, vcopt) 212 213if (vcopt%iwgrec.eq.1) then 214 ! Id weighting field for gradient 215 call field_get_key_int(iflid, kwgrec, iflwgr) 216 call field_get_dim(iflwgr, f_dim) 217 if (f_dim.gt.1) then 218 call field_get_val_v(iflwgr, cpro_wgrec_v) 219 else 220 call field_get_val_s(iflwgr, cpro_wgrec_s) 221 endif 222endif 223 224! Allocate temporary arrays 225allocate(smbrv(3,ncelet), fimp(3,3,ncelet)) 226allocate(weighf(2,nfac)) 227allocate(weighb(nfabor)) 228 229if (ippmod(idarcy).eq.1) then 230 allocate(diverg(ncelet)) 231endif 232 233! --- Numero du scalaire eventuel associe dans le cas fluctuation 234! et numero de variable de calcul 235iiscav = iscavr(iscal) 236if (iiscav.gt.0.and.iiscav.le.nscal) then 237 ivarsc = isca(iiscav) 238else 239 ivarsc = 0 240endif 241 242! --- Numero des grandeurs physiques 243call field_get_val_s(icrom, crom) 244call field_have_previous(icrom, lprev) 245if (lprev) then 246 call field_get_val_prev_s(icrom, croma) 247endif 248call field_get_val_s(ivisct, visct) 249 250call field_get_key_int (ivarfl(isca(iscal)), kivisl, ifcvsl) 251if (ifcvsl.ge.0) then 252 call field_get_val_s(ifcvsl, cpro_viscls) 253endif 254 255! --- Numero de propriété du terme source si extrapolation 256call field_get_key_int(iflid, kstprv, st_prv_id) 257if (st_prv_id .ge.0) then 258 call field_get_val_v(st_prv_id, cproa_vect_st) 259else 260 cproa_vect_st => null() 261endif 262 263! S pour Source, V pour Variable 264thets = thetss(iscal) 265thetv = vcopt%thetav 266 267call field_get_name(iflid, chaine) 268 269if (vcopt%iwarni.ge.1) then 270 write(nfecra,1000) chaine(1:16) 271endif 272 273! Retrieve turbulent Schmidt value for current vector 274call field_get_key_double(ivarfl(isca(iscal)), ksigmas, turb_schmidt) 275 276! Reference diffusivity 277call field_get_key_double(ivarfl(isca(iscal)), kvisl0, visls_0) 278 279!=============================================================================== 280! 2. Source terms 281!=============================================================================== 282 283! --> Initialization 284 285do iel = 1, ncel 286 do isou = 1, 3 287 do jsou = 1, 3 288 fimp(isou,jsou,iel) = 0.d0 289 enddo 290 smbrv(isou,iel) = 0.d0 291 enddo 292enddo 293 294call ustsvv & 295( nvar , nscal , ncepdp , ncesmp , & 296 iscal , & 297 icepdc , icetsm , itypsm , & 298 dt , & 299 ckupdc , smacel , smbrv , fimp ) 300 301! C version 302call user_source_terms(ivarfl(isca(iscal)), smbrv, fimp) 303 304! Store the source terms for convective limiter 305call field_get_key_int(iflid, kst, st_id) 306if (st_id .ge.0) then 307 call field_get_dim(st_id, f_dim) 308 if (f_dim.ne.3) then 309 call csexit(1) 310 endif 311 call field_get_val_v(st_id, cpro_vect_st) 312 313 do iel = 1, ncel 314 !Fill the scalar source term field 315 do isou = 1, 3 316 cpro_vect_st(isou,iel) = smbrv(isou,iel) 317 enddo 318 end do 319 ! Handle parallelism and periodicity 320 if (irangp.ge.0.or.iperio.eq.1) then 321 call synvin(cpro_vect_st) 322 endif 323end if 324 325! Si on extrapole les TS : 326! SMBRV recoit -theta TS du pas de temps precedent 327! (on aurait pu le faire avant ustsvv, mais avec le risque que 328! l'utilisateur l'ecrase) 329! SMBRV recoit la partie du terme source qui depend de la variable 330! A l'ordre 2, on suppose que le ROVSDT fourni par l'utilisateur est <0 331! on implicite le terme (donc ROVSDT*RTPA va dans SMBRV) 332! En std, on adapte le traitement au signe de ROVSDT, mais ROVSDT*RTPA va 333! quand meme dans SMBRV (pas d'autre choix) 334if (st_prv_id .ge. 0) then 335 do iel = 1, ncel 336 ! Stockage temporaire pour economiser un tableau 337 do isou = 1, 3 338 smbexp(isou) = cproa_vect_st(isou,iel) 339 340 ! Terme source utilisateur explicite 341 cproa_vect_st(isou,iel) = smbrv(isou,iel) 342 ! Terme source du pas de temps precedent et 343 ! On suppose -fimp(isou,isou) > 0 : on implicite 344 ! le terme source utilisateur (le reste) 345 smbrv(isou,iel) = fimp(isou,isou,iel)*cvara_var(isou,iel) & 346 - thets*smbexp(isou) 347 ! Diagonale 348 fimp(isou,isou,iel) = - thetv*fimp(isou,isou,iel) 349 enddo 350 enddo 351 352! Si on n'extrapole pas les TS : 353else 354 do iel = 1, ncel 355 do isou = 1, 3 356 ! Terme source utilisateur 357 smbrv(isou,iel) = smbrv(isou,iel) & 358 + fimp(isou,isou,iel)*cvara_var(isou,iel) 359 ! Diagonale 360 fimp(isou,isou,iel) = max(-fimp(isou,isou,iel),zero) 361 enddo 362 enddo 363endif 364 365! --> Physique particulieres 366! Ordre 2 non pris en compte 367 368if (ippmod(iphpar).ge.1) then 369 call pptsvv(iscal, smbrv) 370endif 371 372! Mass source term 373if (ncesmp.gt.0) then 374 375 ! Entier egal a 1 (pour navsto : nb de sur-iter) 376 iiun = 1 377 378 ! On incremente SMBRV par -Gamma RTPA et ROVSDT par Gamma 379 allocate(gavinj(3,ncelet)) 380 call catsmv(ncesmp, 1, icetsm, itypsm(:,ivar), & 381 cell_f_vol, cvara_var, smacel(:,ivar), smacel(:,ipr), & 382 smbrv, fimp, gavinj) 383 384 ! Si on extrapole les TS on met Gamma Pinj dans cproa_vect_st 385 if (st_prv_id .ge. 0) then 386 do iel = 1, ncel 387 do isou = 1, 3 388 cproa_vect_st(isou,iel) = cproa_vect_st(isou,iel) + gavinj(isou,iel) 389 enddo 390 enddo 391 ! Sinon on le met directement dans SMBRV 392 else 393 do iel = 1, ncel 394 do isou = 1, 3 395 smbrv(isou,iel) = smbrv(isou,iel) + gavinj(isou,iel) 396 enddo 397 enddo 398 endif 399 400endif 401 402if (st_prv_id .ge. 0) then 403 thetp1 = 1.d0 + thets 404 do iel = 1, ncel 405 do isou = 1, 3 406 smbrv(isou,iel) = smbrv(isou,iel) + thetp1 * cproa_vect_st(isou,iel) 407 enddo 408 enddo 409endif 410 411! Compressible algorithm 412! or Low Mach compressible algos with mass flux prediction 413if (ippmod(icompf).ge.0.or.(idilat.gt.1.and.ipredfl.eq.1.and.irovar.eq.1)) then 414 pcrom => croma 415 416! Low Mach compressible algos (conservative in time). 417! Same algo for Volume of Fluid method. 418else if ((idilat.gt.1.or.ivofmt.gt.0).and.irovar.eq.1) then 419 if (iterns.eq.1) then 420 call field_get_val_prev2_s(icrom, pcrom) 421 else 422 call field_get_val_prev_s(icrom, pcrom) 423 endif 424 425! Deprecated algo or constant density 426else 427 call field_get_val_s(icrom, pcrom) 428endif 429 430! Get the the order of the diffusivity tensor of the variable 431call field_get_key_int_by_name(iflid, "diffusivity_tensor", idftnp) 432 433! "VITESSE" DE DIFFUSION FACETTE 434 435! On prend le MAX(mu_t,0) car en LES dynamique mu_t peut etre negatif 436! (clipping sur (mu + mu_t)). On aurait pu prendre 437! MAX(K + K_t,0) mais cela autoriserait des K_t negatif, ce qui est 438! considere ici comme non physique. 439if (vcopt%idiff.ge.1) then 440 ! Scalar diffusivity 441 if (iand(vcopt%idften, ISOTROPIC_DIFFUSION).ne.0) then 442 443 allocate(w1(ncelet)) 444 idifftp = vcopt%idifft 445 if (ifcvsl.lt.0) then 446 do iel = 1, ncel 447 w1(iel) = visls_0 & 448 + idifftp*max(visct(iel),zero)/turb_schmidt 449 enddo 450 else 451 do iel = 1, ncel 452 w1(iel) = cpro_viscls(iel) & 453 + idifftp*max(visct(iel),zero)/turb_schmidt 454 enddo 455 endif 456 457 if (vcopt%iwgrec.eq.1) then 458 ! Weighting for gradient 459 do iel = 1, ncel 460 cpro_wgrec_s(iel) = w1(iel) 461 enddo 462 call synsca(cpro_wgrec_s) 463 endif 464 465 imvisp = vcopt%imvisf 466 467 call viscfa & 468 ( imvisp , & 469 w1 , & 470 viscf , viscb ) 471 472 deallocate(w1) 473 474 ! Symmetric tensor diffusivity (GGDH) 475 elseif (iand(vcopt%idften, ANISOTROPIC_DIFFUSION).ne.0) then 476 477 ! Allocate temporary arrays 478 allocate(viscce(6,ncelet)) 479 480 if (iturb.ne.32) then 481 call field_get_val_v(ivsten, visten) 482 else ! EBRSM and (GGDH or AFM) 483 call field_get_val_v(ivstes, visten) 484 endif 485 486 if (ifcvsl.lt.0) then 487 do iel = 1, ncel 488 489 temp = vcopt%idifft*ctheta(iscal)/csrij 490 viscce(1,iel) = temp*visten(1,iel) + visls_0 491 viscce(2,iel) = temp*visten(2,iel) + visls_0 492 viscce(3,iel) = temp*visten(3,iel) + visls_0 493 viscce(4,iel) = temp*visten(4,iel) 494 viscce(5,iel) = temp*visten(5,iel) 495 viscce(6,iel) = temp*visten(6,iel) 496 497 enddo 498 else 499 do iel = 1, ncel 500 501 temp = vcopt%idifft*ctheta(iscal)/csrij 502 viscce(1,iel) = temp*visten(1,iel) + cpro_viscls(iel) 503 viscce(2,iel) = temp*visten(2,iel) + cpro_viscls(iel) 504 viscce(3,iel) = temp*visten(3,iel) + cpro_viscls(iel) 505 viscce(4,iel) = temp*visten(4,iel) 506 viscce(5,iel) = temp*visten(5,iel) 507 viscce(6,iel) = temp*visten(6,iel) 508 509 enddo 510 endif 511 512 if (vcopt%iwgrec.eq.1) then 513 ! Weighting for gradient 514 do iel = 1, ncel 515 do isou = 1, 6 516 cpro_wgrec_v(isou,iel) = viscce(isou,iel) 517 enddo 518 enddo 519 call syntis(cpro_wgrec_v) 520 endif 521 522 iwarnp = vcopt%iwarni 523 524 call vitens & 525 ( viscce , iwarnp , & 526 weighf , weighb , & 527 viscf , viscb ) 528 529 endif 530 531else 532 533 do ifac = 1, nfac 534 viscf(ifac) = 0.d0 535 enddo 536 do ifac = 1, nfabor 537 viscb(ifac) = 0.d0 538 enddo 539 540endif 541 542if (ippmod(idarcy).eq.1) then 543 call field_get_key_struct_gwf_soilwater_partition(iflid, sorption_scal) 544 call field_get_val_s(sorption_scal%idel, cpro_delay) 545 call field_get_val_prev_s(sorption_scal%idel, cproa_delay) 546 call field_get_val_prev_s_by_name('saturation', cproa_sat) 547 call field_get_val_s_by_name('saturation', cpro_sat) 548endif 549 550! Not Darcy 551if (ippmod(idarcy).eq.-1) then 552 ! --> Unsteady term and mass aggregation term 553 do iel = 1, ncel 554 do isou = 1, 3 555 fimp(isou,isou,iel) = fimp(isou,isou,iel) & 556 + vcopt%istat*pcrom(iel)*cell_f_vol(iel)/dt(iel) 557 enddo 558 enddo 559! Darcy : we take into account the porosity and delay for underground transport 560else 561 do iel = 1, ncel 562 do isou = 1, 3 563 smbrv(isou,iel) = smbrv(isou,iel)*cpro_delay(iel)*cpro_sat(iel) 564 fimp(isou,isou,iel) = (fimp(isou,isou,iel) & 565 + vcopt%istat*pcrom(iel)*volume(iel)/dt(iel) ) & 566 * cpro_delay(iel)*cpro_sat(iel) 567 enddo 568 enddo 569endif 570 571! Scalar with a Drift: 572! compute the convective flux 573!---------------------------- 574 575call field_get_key_int(iflid, keydri, iscdri) 576 577if (iscdri.ge.1) then 578 allocate(divflu(ncelet)) 579 580 call driflu & 581 ( iflid , & 582 dt , & 583 imasfl , bmasfl , & 584 divflu ) 585 586 iconvp = vcopt%iconv 587 thetap = vcopt%thetav 588 589 ! NB: if the porosity module is swiched on, the the porosity is already 590 ! taken into account in divflu 591 592 ! --> mass aggregation term 593 do iel = 1, ncel 594 do isou = 1, 3 595 fimp(isou,isou,iel) = fimp(isou,isou,iel) + iconvp*thetap*divflu(iel) 596 smbrv(isou,iel) = smbrv(isou,iel) - iconvp*divflu(iel)*cvara_var(isou,iel) 597 enddo 598 enddo 599 600 deallocate(divflu) 601endif 602 603! Darcy: 604! This step is necessary because the divergence of the 605! hydraulic head (or pressure), which is taken into account as the 606! 'aggregation term' via 607! cs_equation_iterative_solve_scalar -> cs_balance_scalar, 608! is not exactly equal to the loss of mass of water in the unsteady case 609! (just as far as the precision of the Newton scheme is good), 610! and does not take into account the sorbption of the tracer 611! (represented by the 'delay' coefficient). 612! We choose to 'cancel' the aggregation term here, and to add to smbr 613! (right hand side of the transport equation) 614! the necessary correction to take sorption into account and get the exact 615! conservation of the mass of tracer. 616 617if (ippmod(idarcy).eq.1) then 618 if (darcy_unsteady.eq.1) then 619 call divmas(1, imasfl , bmasfl , diverg) 620 do iel = 1, ncel 621 do isou = 1, 3 622 smbrv(isou,iel) = smbrv(isou,iel) - diverg(iel)*cvar_var(isou,iel) & 623 + cell_f_vol(iel)/dt(iel)*cvar_var(isou,iel) & 624 *( cproa_delay(iel)*cproa_sat(iel) & 625 - cpro_delay(iel)*cpro_sat(iel) ) 626 enddo 627 enddo 628 do iel = 1, ncel 629 do isou = 1, 3 630 fimp(isou,isou,iel) = fimp(isou,isou,iel) + thetv*diverg(iel) 631 enddo 632 enddo 633 endif 634endif 635 636!=============================================================================== 637! 3. Solving 638!=============================================================================== 639 640iescap = 0 641! all boundary convective flux with upwind 642icvflb = 0 643! transposed gradient term only for NS 644ivissv = 0 645 646call field_get_coefa_v(iflid, coefap) 647call field_get_coefb_v(iflid, coefbp) 648call field_get_coefaf_v(iflid, cofafp) 649call field_get_coefbf_v(iflid, cofbfp) 650 651vcopt_loc = vcopt 652 653vcopt_loc%istat = -1 654vcopt_loc%idifft = -1 655vcopt_loc%iwgrec = 0 656vcopt_loc%thetav = thetv 657vcopt_loc%blend_st = 0 ! Warning, may be overwritten if a field 658 659p_k_value => vcopt_loc 660c_k_value = equation_param_from_vcopt(c_loc(p_k_value)) 661 662call cs_equation_iterative_solve_vector & 663 ( idtvar , iterns , & 664 iflid , c_null_char , & 665 ivissv , iescap , c_k_value , & 666 cvara_var , cvara_var , & 667 coefap , coefbp , cofafp , cofbfp , & 668 imasfl , bmasfl , viscf , & 669 viscb , viscf , viscb , rvoid , & 670 rvoid , viscce , weighf , weighb , & 671 icvflb , ivoid , & 672 fimp , smbrv , cvar_var , rvoid ) 673 674!=============================================================================== 675! 4. Writing 676!=============================================================================== 677 678! BILAN EXPLICITE (VOIR CODITS : ON ENLEVE L'INCREMENT) 679! Ceci devrait etre valable avec le theta schema sur les Termes source 680 681if (vcopt%iwarni.ge.2) then 682 if (vcopt%nswrsm.gt.1) then 683 ibcl = 1 684 else 685 ibcl = 0 686 endif 687 do iel = 1, ncel 688 do isou = 1, 3 689 smbrv(isou,iel) = smbrv(isou,iel) & 690 - vcopt%istat*(pcrom(iel)/dt(iel))*cell_f_vol(iel) & 691 *(cvar_var(isou,iel)-cvara_var(isou,iel))*ibcl 692 enddo 693 enddo 694 sclnor = sqrt(cs_gdot(3*ncel,smbrv,smbrv)) 695 write(nfecra,1200)chaine(1:16) ,sclnor 696endif 697 698! Free memory 699deallocate(smbrv, fimp) 700deallocate(weighf, weighb) 701if (allocated(viscce)) deallocate(viscce) 702if (allocated(diverg)) deallocate(diverg) 703 704!-------- 705! Formats 706!-------- 707 708 1000 format(/, & 709' ** SOLVING VARIABLE ',A16 ,/,& 710' ----------------' ,/) 711 1200 format(1X,A16,' : EXPLICIT BALANCE = ',E14.5) 712! 9000 format( & 713!'@' ,/,& 714!'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 715!'@' ,/,& 716!'@ @@ WARNING: ERROR IN COVOFV' ,/,& 717!'@ ========' ,/,& 718!'@ IVARSC MUST BE A STRICTLY POSITIVE INTEGER' ,/,& 719!'@ ITS VALUE IS ',I10 ,/,& 720!'@' ,/,& 721!'@ The calculation will not be run.' ,/,& 722!'@' ,/,& 723!'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 724!'@ ',/) 725 726!---- 727! End 728!---- 729 730return 731 732end subroutine 733