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!> \file modini.f90 24!> \brief Modify calculation parameters after user changes (module variables) 25!> 26!------------------------------------------------------------------------------ 27 28subroutine modini 29 30!=============================================================================== 31! Module files 32!=============================================================================== 33 34use paramx 35use cstnum 36use dimens 37use field 38use numvar 39use optcal 40use cstphy 41use entsor 42use albase 43use alstru 44use cplsat 45use post 46use ppincl 47use rotation 48use darcy_module 49use turbomachinery 50use vof 51use cs_c_bindings 52 53!=============================================================================== 54 55implicit none 56 57! Arguments 58 59 60! Local variables 61 62integer f_id 63integer ii, jj, iok, ikw 64integer nbccou 65integer nscacp, iscal, ivar 66integer imrgrp 67integer iscacp, kcpsyr, icpsyr 68integer nfld, f_type 69integer key_t_ext_id, icpext, kscmin, kscmax 70integer iviext 71integer kturt, turb_flux_model, turb_flux_model_type 72 73double precision relxsp, clvfmn, clvfmx, visls_0, visls_cmp 74double precision scminp 75 76character(len=80) :: name 77 78type(var_cal_opt) :: vcopt , vcopt1 79 80!=============================================================================== 81 82! Indicateur erreur (0 : pas d'erreur) 83iok = 0 84 85call field_get_key_id("syrthes_coupling", kcpsyr) 86 87call field_get_key_id("time_extrapolated", key_t_ext_id) 88 89call field_get_key_id("min_scalar_clipping", kscmin) 90call field_get_key_id("max_scalar_clipping", kscmax) 91call field_get_key_id('turbulent_flux_model', kturt) 92 93!=============================================================================== 94! 1. ENTREES SORTIES entsor 95!=============================================================================== 96 97!---> sorties chrono? 98 99if (idtvar.lt.0) then 100 call hide_property(icour) 101 call hide_property(ifour) 102endif 103 104!---> sorties historiques ? 105! Si une valeur non modifiee par l'utilisateur (=-999) 106! on la met a sa valeur par defaut 107! On sort toutes les variables a tous les pas de temps par defaut 108! NTHIST = -1 : on ne sort pas d'historiques 109! NTHIST = n : on sort des historiques tous les n pas de temps 110 111! Adapt the output frequency parameters according to the time scheme. 112if (idtvar.lt.0.or.idtvar.eq.2) then 113 frhist = -1.d0 114else 115 if (frhist > 0.d0) then 116 nthist = -1 117 endif 118endif 119 120! Logging and postprocessing output 121 122if (irovar.eq.0) then 123 call hide_property(icrom) 124 call hide_property(ibrom) 125endif 126 127if (ivivar.eq.0) then 128 call hide_property(iviscl) 129endif 130 131if (idtvar.lt.0) then 132 call hide_property(icour) 133 call hide_property(ifour) 134endif 135 136!=============================================================================== 137! 2. POSITION DES VARIABLES DE numvar 138!=============================================================================== 139 140! ---> Reperage des variables qui disposeront de deux types de CL 141 142! Fait dans varpos. 143! Si l'utilisateur y a touche ensuite, on risque l'incident. 144 145!=============================================================================== 146! 3. OPTIONS DU CALCUL : TABLEAUX DE optcal 147!=============================================================================== 148 149! time scheme 150 151if (ntmabs.eq.-1 .and. ttmabs.lt.-0.5) then 152 ntmabs = 10 153endif 154 155! restart 156 157call indsui(isuite) 158 159if (isuit1.eq.-1) isuit1 = isuite 160 161! -- Proprietes physiques 162call field_get_key_int(iviscl, key_t_ext_id, iviext) 163if (abs(thetvi+999.d0).gt.epzero) then 164 write(nfecra,1011) 'IVIEXT',iviext,'THETVI' 165 iok = iok + 1 166elseif (iviext.eq.0) then 167 thetvi = 0.0d0 168elseif (iviext.eq.1) then 169 thetvi = 0.5d0 170elseif (iviext.eq.2) then 171 thetvi = 1.d0 172endif 173 174if (icp.ge.0) then 175 call field_get_key_int(icp, key_t_ext_id, icpext) 176 if (abs(thetcp+999.d0).gt.epzero) then 177 write(nfecra,1011) 'ICPEXT',icpext,'THETCP' 178 iok = iok + 1 179 elseif (icpext.eq.0) then 180 thetcp = 0.0d0 181 elseif (icpext.eq.1) then 182 thetcp = 0.5d0 183 elseif (icpext.eq.2) then 184 thetcp = 1.d0 185 endif 186endif 187 188! -- Termes sources NS 189if (abs(thetsn+999.d0).gt.epzero) then 190 write(nfecra,1011) 'ISNO2T',isno2t,'THETSN' 191 iok = iok + 1 192elseif (isno2t.eq.1) then 193 thetsn = 0.5d0 194elseif (isno2t.eq.2) then 195 thetsn = 1.d0 196elseif (isno2t.eq.0) then 197 thetsn = 0.d0 198endif 199 200! -- Termes sources grandeurs turbulentes 201if (abs(thetst+999.d0).gt.epzero) then 202 write(nfecra,1011) 'ISTO2T',isto2t,'THETST' 203 iok = iok + 1 204elseif (isto2t.eq.1) then 205 thetst = 0.5d0 206elseif (isto2t.eq.2) then 207 thetst = 1.d0 208elseif (isto2t.eq.0) then 209 thetst = 0.d0 210endif 211 212do iscal = 1, nscal 213! -- Termes sources des scalaires 214 if (abs(thetss(iscal)+999.d0).gt.epzero) then 215 write(nfecra,1021) iscal,'ISSO2T',isso2t(iscal),'THETSS' 216 iok = iok + 1 217 elseif (isso2t(iscal).eq.1) then 218 thetss(iscal) = 0.5d0 219 elseif (isso2t(iscal).eq.2) then 220 thetss(iscal) = 1.d0 221 elseif (isso2t(iscal).eq.0) then 222 thetss(iscal) = 0.d0 223 endif 224 ! Scalars diffusivity 225 call field_get_key_int(ivarfl(isca(iscal)), kivisl, f_id) 226 if (f_id.ge.0) then 227 call field_get_key_int(f_id, key_t_ext_id, iviext) 228 else 229 iviext = 0 230 endif 231 if (abs(thetvs(iscal)+999.d0).gt.epzero) then 232 write(nfecra,1021) iscal,'IVSEXT',iviext,'THETVS' 233 iok = iok + 1 234 elseif (iviext.eq.0) then 235 thetvs(iscal) = 0.d0 236 elseif (iviext.eq.1) then 237 thetvs(iscal) = 0.5d0 238 elseif (iviext.eq.2) then 239 thetvs(iscal) = 1.d0 240 endif 241enddo 242 243! Loop on on field variables 244call field_get_n_fields(nfld) 245 246do f_id = 0, nfld - 1 247 call field_get_type(f_id, f_type) 248 ! Is the field of type FIELD_VARIABLE? 249 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 250 call field_get_key_struct_var_cal_opt(f_id, vcopt) 251 if (abs(vcopt%thetav+1.d0).gt.epzero) then 252 call field_get_name(f_id, name) 253 write(nfecra,1131) trim(name),'THETAV' 254 else 255 if (vcopt%istat.eq.0) then 256 vcopt%thetav = 1.d0 257 else if (ischtp.eq.1) then 258 vcopt%thetav = 1.d0 259 else if (ischtp.eq.2) then 260 vcopt%thetav = 0.5d0 261 endif 262 endif 263 call field_set_key_struct_var_cal_opt(f_id, vcopt) 264 endif 265enddo 266 267! Diffusivity model: 268! Daly Harlow (GGDH) on Rij and epsilon by default 269if (itytur.eq.3) then 270 271 call field_get_key_struct_var_cal_opt(ivarfl(irij), vcopt1) 272 call field_get_key_struct_var_cal_opt(ivarfl(iep), vcopt) 273 274 ! Diffusivity model: 275 ! Daly Harlow (GGDH) on Rij and epsilon by default 276 if (idirsm.ne.0) then 277 vcopt1%idften = ANISOTROPIC_RIGHT_DIFFUSION 278 vcopt%idften = ANISOTROPIC_RIGHT_DIFFUSION 279 ! Scalar diffusivity (Shir model) elswhere (idirsm = 0) 280 else 281 vcopt1%idften = ISOTROPIC_DIFFUSION 282 vcopt%idften = ISOTROPIC_DIFFUSION 283 endif 284 285 call field_set_key_struct_var_cal_opt(ivarfl(irij), vcopt1) 286 call field_set_key_struct_var_cal_opt(ivarfl(iep), vcopt) 287endif 288 289! ---> ISSTPC 290! Si l'utilisateur n'a rien specifie pour le test de pente (=-1), 291! On impose 1 (ie sans) pour la vitesse en LES 292! 0 (ie avec) sinon 293 294if (itytur.eq.4) then 295 ii = iu 296 call field_get_key_struct_var_cal_opt(ivarfl(ii), vcopt) 297 if (vcopt%isstpc.eq.-999) then 298 vcopt%isstpc = 1 299 call field_set_key_struct_var_cal_opt(ivarfl(ii), vcopt) 300 endif 301 do jj = 1, nscal 302 ii = isca(jj) 303 call field_get_key_struct_var_cal_opt(ivarfl(ii), vcopt) 304 if (vcopt%isstpc.eq.-999) then 305 vcopt%isstpc = 0 306 call field_set_key_struct_var_cal_opt(ivarfl(ii), vcopt) 307 endif 308 enddo 309endif 310 311do f_id = 0, nfld - 1 312 call field_get_type(f_id, f_type) 313 ! Is the field of type FIELD_VARIABLE? 314 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 315 call field_get_key_struct_var_cal_opt(f_id, vcopt) 316 if (vcopt%isstpc.eq.-999) then 317 vcopt%isstpc = 0 318 call field_set_key_struct_var_cal_opt(f_id, vcopt) 319 endif 320 endif 321enddo 322 323! ---> BLENCV 324! Si l'utilisateur n'a rien specifie pour le schema convectif 325! 1 (ie centre) pour les vitesses 326! les scalaires utilisateurs 327! le scalaire thermique 328! 0 (ie upwind pur) pour le reste 329! (en particulier, en L.E.S. toutes les variables sont donc en centre) 330 331! Pour le modele de cavitation on force dans tous les cas le taux de vide en 332! upwind et on affiche un message si l'utilisateur avait specifie autre chose 333 334ii = iu 335call field_get_key_struct_var_cal_opt(ivarfl(ii), vcopt) 336if (abs(vcopt%blencv+1.d0).lt.epzero) then 337 vcopt%blencv = 1.d0 338 call field_set_key_struct_var_cal_opt(ivarfl(ii), vcopt) 339endif 340 341if (iand(ivofmt,VOF_MERKLE_MASS_TRANSFER).ne.0) then 342 ii = ivolf2 343 call field_get_key_struct_var_cal_opt(ivarfl(ii), vcopt) 344 if (abs(vcopt%blencv+1.d0).lt.epzero) then 345 if (abs(vcopt%blencv+1.d0).gt.epzero) & 346 write(nfecra,3000) 0.d0, vcopt%blencv 347 vcopt%blencv = 0.d0 348 call field_set_key_struct_var_cal_opt(ivarfl(ii), vcopt) 349 endif 350else if (ivofmt.gt.0) then 351 ii = ivolf2 352 call field_get_key_struct_var_cal_opt(ivarfl(ii), vcopt) 353 if (abs(vcopt%blencv+1.d0).lt.epzero) then 354 vcopt%blencv = 1.d0 355 call field_set_key_struct_var_cal_opt(ivarfl(ii), vcopt) 356 endif 357endif 358 359do jj = 1, nscaus 360 ii = isca(jj) 361 call field_get_key_struct_var_cal_opt(ivarfl(ii), vcopt) 362 if (abs(vcopt%blencv+1.d0).lt.epzero) then 363 vcopt%blencv = 1.d0 364 call field_set_key_struct_var_cal_opt(ivarfl(ii), vcopt) 365 endif 366enddo 367 368if (iscalt.gt.0) then 369 ii = isca(iscalt) 370 call field_get_key_struct_var_cal_opt(ivarfl(ii), vcopt) 371 if (abs(vcopt%blencv+1.d0).lt.epzero) then 372 vcopt%blencv = 1.d0 373 call field_set_key_struct_var_cal_opt(ivarfl(ii), vcopt) 374 endif 375endif 376 377do f_id = 0, nfld - 1 378 call field_get_type(f_id, f_type) 379 ! Is the field of type FIELD_VARIABLE? 380 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 381 call field_get_key_struct_var_cal_opt(f_id, vcopt) 382 if (abs(vcopt%blencv+1.d0).lt.epzero) then 383 vcopt%blencv = 0.d0 384 call field_set_key_struct_var_cal_opt(f_id, vcopt) 385 endif 386 endif 387enddo 388 389! ---> NSWRSM, EPSRSM ET EPSILO 390! Si l'utilisateur n'a rien specifie (NSWRSM=-1), 391! On impose 392! a l'ordre 1 : 393! 2 pour la pression 394! 1 pour les autres variables 395! on initialise EPSILO a 1.d-8 396! pour la pression 397! on initialise EPSILO a 1.d-5 398! pour les autres variables 399! on initialise EPSRSM a 10*EPSILO 400! a l'ordre 2 : 401! 5 pour la pression 402! 10 pour les autres variables 403! on initialise EPSILO a 1.D-5 404! on initialise EPSRSM a 10*EPSILO 405! Attention aux tests dans verini 406 407if (ischtp.eq.2) then 408 ii = ipr 409 call field_get_key_struct_var_cal_opt(ivarfl(ii), vcopt) 410 if (vcopt%nswrsm.eq.-1) vcopt%nswrsm = 5 411 if (abs(vcopt%epsilo+1.d0).lt.epzero) vcopt%epsilo = 1.d-5 412 if (abs(vcopt%epsrsm+1.d0).lt.epzero) vcopt%epsrsm = 10.d0*vcopt%epsilo 413 call field_set_key_struct_var_cal_opt(ivarfl(ii), vcopt) 414 ii = iu 415 call field_get_key_struct_var_cal_opt(ivarfl(ii), vcopt) 416 if (vcopt%nswrsm.eq.-1) vcopt%nswrsm = 10 417 if (abs(vcopt%epsilo+1.d0).lt.epzero) vcopt%epsilo = 1.d-5 418 if (abs(vcopt%epsrsm+1.d0).lt.epzero) vcopt%epsrsm = 10.d0*vcopt%epsilo 419 call field_set_key_struct_var_cal_opt(ivarfl(ii), vcopt) 420 do jj = 1, nscal 421 ii = isca(jj) 422 call field_get_key_struct_var_cal_opt(ivarfl(ii), vcopt) 423 if (vcopt%nswrsm.eq.-1) vcopt%nswrsm = 10 424 if (abs(vcopt%epsilo+1.d0).lt.epzero) vcopt%epsilo = 1.d-5 425 if (abs(vcopt%epsrsm+1.d0).lt.epzero) vcopt%epsrsm = 10.d0*vcopt%epsilo 426 call field_set_key_struct_var_cal_opt(ivarfl(ii), vcopt) 427 enddo 428endif 429 430! For the pressure, default solver precision 1e-8 431! because the mass conservation is up to this precision 432ii = ipr 433call field_get_key_struct_var_cal_opt(ivarfl(ii), vcopt) 434if (vcopt%nswrsm.eq.-1) vcopt%nswrsm = 2 435if (abs(vcopt%epsilo+1.d0).lt.epzero) vcopt%epsilo = 1.d-8 436call field_set_key_struct_var_cal_opt(ivarfl(ii), vcopt) 437 438do f_id = 0, nfld - 1 439 call field_get_type(f_id, f_type) 440 ! Is the field of type FIELD_VARIABLE? 441 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 442 call field_get_key_struct_var_cal_opt(f_id, vcopt) 443 if (vcopt%nswrsm.eq.-1) vcopt%nswrsm = 1 444 if (abs(vcopt%epsilo+1.d0).lt.epzero) vcopt%epsilo = 1.d-5 445 if (abs(vcopt%epsrsm+1.d0).lt.epzero) vcopt%epsrsm = 10.d0*vcopt%epsilo 446 call field_set_key_struct_var_cal_opt(f_id, vcopt) 447 endif 448enddo 449 450! ---> IMLIGR 451! Si l'utilisateur n'a rien specifie pour la limitation des 452! gradients (=-999), 453! On impose -1 avec gradrc (pas de limitation) 454! et 1 avec gradmc (limitation) 455imrgrp = abs(imrgra) 456if (imrgrp.ge.10) imrgrp = imrgrp - 10 457 458if (imrgrp.eq.0.or.imrgrp.ge.4) then 459 do f_id = 0, nfld - 1 460 call field_get_type(f_id, f_type) 461 ! Is the field of type FIELD_VARIABLE? 462 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 463 call field_get_key_struct_var_cal_opt(f_id, vcopt) 464 if (vcopt%imligr.eq.-999) then 465 vcopt%imligr = -1 466 call field_set_key_struct_var_cal_opt(f_id, vcopt) 467 endif 468 endif 469 enddo 470else 471 do f_id = 0, nfld - 1 472 call field_get_type(f_id, f_type) 473 ! Is the field of type FIELD_VARIABLE? 474 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 475 call field_get_key_struct_var_cal_opt(f_id, vcopt) 476 if (vcopt%imligr.eq.-999) then 477 vcopt%imligr = 1 478 call field_set_key_struct_var_cal_opt(f_id, vcopt) 479 endif 480 endif 481 enddo 482endif 483 484! ---> DTMIN DTMAX CDTVAR 485 486if (dtmin.le.-grand) then 487 dtmin = 0.1d0*dtref 488endif 489if (dtmax.le.-grand) then 490 dtmax = 1.0d3*dtref 491endif 492 493! Init. of time step factor for velocity, pressure and turbulent variables 494! FIXME time step factor is used ONLY for additional variables (user or model) 495 496cdtvar(iv ) = cdtvar(iu) 497cdtvar(iw ) = cdtvar(iu) 498cdtvar(ipr) = cdtvar(iu) 499 500if (itytur.eq.2) then 501 cdtvar(iep ) = cdtvar(ik ) 502elseif (itytur.eq.3) then 503 cdtvar(ir22) = cdtvar(ir11) 504 cdtvar(ir33) = cdtvar(ir11) 505 cdtvar(ir12) = cdtvar(ir11) 506 cdtvar(ir13) = cdtvar(ir11) 507 cdtvar(ir23) = cdtvar(ir11) 508 cdtvar(iep ) = cdtvar(ir11) 509 ! cdtvar(ial) is useless because no time dependance in the equation of alpha. 510 if (iturb.eq.32) then 511 cdtvar(ial) = cdtvar(ir11) 512 endif 513elseif (itytur.eq.5) then 514 cdtvar(iep ) = cdtvar(ik ) 515 cdtvar(iphi) = cdtvar(ik ) 516! CDTVAR(IFB/IAL) est en fait inutile car pas de temps dans 517! l'eq de f_barre/alpha 518 if (iturb.eq.50) then 519 cdtvar(ifb ) = cdtvar(ik ) 520 elseif (iturb.eq.51) then 521 cdtvar(ial ) = cdtvar(ik ) 522 endif 523elseif (iturb.eq.60) then 524 cdtvar(iomg) = cdtvar(ik ) 525elseif (iturb.eq.70) then 526 ! cdtvar est a 1.0 par defaut dans iniini.f90 527 cdtvar(inusa)= cdtvar(inusa) 528endif 529 530! ---> IWALLF 531! For laminar cases or when using low Reynolds model: no wall function. 532! When using mixing length, Spalart-Allmaras or LES: one scale log law. 533! When using EB-RSM : all y+ wall functions 534! In all other cases: 2 scales log law. 535! Here iwallf is set automatically only if it wasn't set in the gui or 536! in a user subroutine. 537 538if (iwallf.eq.-999) then 539 if ( iturb.eq.10.or.iturb.eq.70 & 540 .or.itytur.eq.4) then 541 iwallf = 2 542 elseif (iturb.eq.0.or.itytur.eq.5) then 543 iwallf = 0 544 elseif (iturb.eq.32) then 545 iwallf = 7 546 else 547 iwallf = 3 548 endif 549endif 550 551! ---> IWALFS 552! If the wall function for the velocity is the two scales wall function using 553! Van Driest mixing length (iwallf=5), then the corresponding wall function for 554! scalar should be used (iwalfs=1). 555! For atmospheric Flows, it is by default Louis, or Monin-Obukhov 556! Here iwalfs is set automatically only if it wasn't set in a user subroutine. 557 558if (iwalfs.eq.-999) then 559 if (ippmod(iatmos).ge.0) then 560 iwalfs = 2 561 else if (iwallf.eq.5) then 562 iwalfs = 1 563 else 564 iwalfs = 0 565 endif 566endif 567 568! ---> YPLULI 569! 1/XKAPPA est la valeur qui assure la continuite de la derivee 570! entre la zone lineaire et la zone logarithmique. 571 572! Dans le cas des lois de paroi invariantes, on utilise la valeur de 573! continuite du profil de vitesse, 10.88. 574 575! Pour la LES, on remet 10.88, afin d'eviter des clic/clac quand on est a 576! la limite (en modele a une echelle en effet, YPLULI=1/XKAPPA ne permet pas 577! forcement de calculer u* de maniere totalement satisfaisante). 578! Idem en Spalart-Allmaras. 579 580if (ypluli.lt.-grand) then 581 if (iwallf.eq.4 .or. itytur.eq.4 .or. iturb.eq.70.or.iwallf.eq.6.or.iturb.eq.60 & 582 .or. iturb.eq.22 ) then 583 ypluli = 10.88d0 584 else 585 ypluli = 1.d0/xkappa 586 endif 587endif 588 589! ---> ICPSYR 590! Si l'utilisateur n'a pas modifie ICPSYR, on prend par defaut : 591! s'il n y a pas de couplage 592! 0 pour tous les scalaires 593! sinon 594! 1 pour le scalaire ISCALT s'il existe 595! 0 pour les autres 596! Les modifs adequates devront etre ajoutees pour les physiques 597! particulieres 598! Les tests de coherence seront faits dans verini. 599 600if (nscal.gt.0) then 601 602! On regarde s'il y a du couplage 603 604 nbccou = cs_syr_coupling_n_couplings() 605 606! S'il y a du couplage 607 if (nbccou .ne. 0) then 608 609! On compte le nombre de scalaires couples 610 nscacp = 0 611 do iscal = 1, nscal 612 call field_get_key_int(ivarfl(isca(iscal)), kcpsyr, icpsyr) 613 if (icpsyr.eq.1) then 614 nscacp = nscacp + 1 615 endif 616 enddo 617 618! Si l'utilisateur n'a pas couple de scalaire, 619 if (nscacp.eq.0) then 620 621! On couple le scalaire temperature de la phase 622 if (iscalt.gt.0.and.iscalt.le.nscal) then 623 icpsyr = 1 624 call field_set_key_int(ivarfl(isca(iscalt)), kcpsyr, icpsyr) 625 goto 100 626 endif 627 100 continue 628 629 endif 630 631 endif 632 633endif 634 635! Temperature scale 636 637if (itherm.ge.1 .and. itpscl.le.0) then 638 itpscl = 1 639endif 640 641! ---> "is_temperature" 642! If the user has not modified "is_temperature", we take by default: 643! passive scalar of scalars other than iscalt 644! = 0 : passive, enthalpy, or energy 645! = 1 : temperature 646 647if (nscal.gt.0) then 648 do ii = 1, nscal 649 call field_get_key_int(ivarfl(isca(ii)), kscacp, iscacp) 650 if (iscacp.eq.-1) then 651 if (ii.eq.iscalt .and. itherm.eq.1) then 652 iscacp = 1 653 else 654 iscacp = 0 655 endif 656 call field_set_key_int(ivarfl(isca(ii)), kscacp, iscacp) 657 endif 658 enddo 659endif 660 661! ---> ICALHY 662! Calcul de la pression hydrostatique en sortie pour les conditions de 663! Dirichlet sur la pression. Se deduit de IPHYDR et de la valeur de 664! la gravite (test assez arbitraire sur la norme). 665! ICALHY est initialise a -1 (l'utilisateur peut avoir force 666! 0 ou 1 et dans ce cas, on ne retouche pas) 667 668if (icalhy.ne.-1.and.icalhy.ne.0.and.icalhy.ne.1) then 669 write(nfecra,1061)icalhy 670 iok = iok + 1 671endif 672 673 674! ---> ICDPAR 675! Calcul de la distance a la paroi. En standard, on met ICDPAR a -1, au cas 676! ou les faces de bord auraient change de type d'un calcul a l'autre. En k-omega, 677! il faut la distance a la paroi pour une suite propre, donc on initialise a 1 et 678! on avertit (dans verini). 679ikw = 0 680if (iturb.eq.60) ikw = 1 681if (icdpar.eq.-999) then 682 icdpar = -1 683 if (ikw.eq.1) icdpar = 1 684 if (isuite.eq.1 .and. ikw.eq.1) write(nfecra,2000) 685endif 686if (icdpar.eq.-1 .and. ikw.eq.1 .and. isuite.eq.1) & 687 write(nfecra,2001) 688 689! ---> IKECOU 690! If the fluid_solid option is enabled, we force ikecou to 0. 691if (fluid_solid) then 692 if(ikecou .eq. 1) then 693 ikecou = 0 694 write(nfecra,5000) 695 endif 696endif 697 698 699! ---> RELAXV 700if (idtvar.lt.0) then 701 relxsp = 1.d0-relxst 702 if (relxsp.le.epzero) relxsp = relxst 703 call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt) 704 if (abs(vcopt%relaxv+1.d0).le.epzero) then 705 vcopt%relaxv = relxsp 706 call field_set_key_struct_var_cal_opt(ivarfl(ipr), vcopt) 707 endif 708 do f_id = 0, nfld - 1 709 call field_get_type(f_id, f_type) 710 ! Is the field of type FIELD_VARIABLE? 711 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 712 call field_get_key_struct_var_cal_opt(f_id, vcopt) 713 if (abs(vcopt%relaxv+1.d0).le.epzero) then 714 vcopt%relaxv = relxst 715 call field_set_key_struct_var_cal_opt(f_id, vcopt) 716 endif 717 endif 718 enddo 719else 720 721 do f_id = 0, nfld - 1 722 call field_get_type(f_id, f_type) 723 ! Is the field of type FIELD_VARIABLE? 724 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 725 call field_get_key_struct_var_cal_opt(f_id, vcopt) 726 if (abs(vcopt%relaxv+1.d0).le.epzero) then 727 vcopt%relaxv = 1.d0 728 call field_set_key_struct_var_cal_opt(f_id, vcopt) 729 endif 730 endif 731 enddo 732endif 733 734! Options specific to steady case 735if (idtvar.lt.0) then 736 ipucou = 0 737 dtref = 1.d0 738 dtmin = 1.d0 739 dtmax = 1.d0 740 do f_id = 0, nfld - 1 741 call field_get_type(f_id, f_type) 742 ! Is the field of type FIELD_VARIABLE? 743 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 744 call field_get_key_struct_var_cal_opt(f_id, vcopt) 745 vcopt%istat = 0 746 call field_set_key_struct_var_cal_opt(f_id, vcopt) 747 endif 748 enddo 749 call field_get_key_struct_var_cal_opt(ivarfl(iu), vcopt) 750 arak = arak/max(vcopt%relaxv,epzero) 751endif 752 753! With a staggered approach no Rhie and Chow correction is needed 754if (staggered.eq.1) then 755 arak = 0.d0 756endif 757 758!=============================================================================== 759! 4. TABLEAUX DE cstphy 760!=============================================================================== 761 762! ---> Constantes 763! Ca fait un calcul en double, mais si qqn a bouge cmu, apow, bpow, 764! ca servira. 765 766cpow = apow**(2.d0/(1.d0-bpow)) 767dpow = 1.d0/(1.d0+bpow) 768 769! Modified value of Cmu for V2f and Bl-v2k 770if (iturb.eq.50.or.iturb.eq.51) cmu = 0.22d0 771 772cmu025 = cmu**0.25d0 773 774if (idirsm.eq.0) then 775 csrij = 0.11d0 776else 777 if (iturb.eq.32) then 778 csrij = 0.21d0 779 else 780 csrij = 0.22d0 781 endif 782endif 783 784! Constant for the Buoyant production term of Rij 785! EBRSM 786if (iturb.eq.32) then 787 crij3 = 0.6d0 788else 789 crij3 = 0.55d0 790endif 791 792if (iturb.eq.60) then !sst-ddes 793 ! SST DDES 794 if (hybrid_turb.eq.2) then 795 cddes = 0.65d0 796 else if (hybrid_turb.eq.1) then 797 cddes = 0.61d0 798 endif 799 ! SST SAS 800 csas = 0.11d0 801 csas_eta2 = 3.51d0 802elseif (iturb.eq.51) then !phif-ddes 803 cddes = 0.60d0 804endif 805 806! ---> ICLVFL 807! Si l'utilisateur n'a pas modifie ICLVFL, on prend par defaut : 808! 0 pour les variances 809! Les modifs adequates devront etre ajoutees pour les physiques 810! particulieres 811! If the user gives a value we put iclcfl to 2. 812 813do iscal = 1, nscal 814 if (iscavr(iscal).gt.0) then 815 816 ! Get the min clipping 817 call field_get_key_double(ivarfl(isca(iscal)), kscmin, scminp) 818 ! If modified put 2 819 if (iclvfl(iscal).eq.-1 .and. abs(scminp+grand).ge.epzero) then 820 iclvfl(iscal) = 2 821 822 else if (iclvfl(iscal).eq.-1) then 823 iclvfl(iscal) = 0 824 endif 825 826 ! Min for variances is 0 or greater 827 call field_get_key_double(ivarfl(isca(iscal)), kscmin, scminp) 828 ! set min clipping to 0 829 scminp = max(0.d0, scminp) 830 call field_set_key_double(ivarfl(isca(iscal)), kscmin, scminp) 831 endif 832enddo 833 834do ii = 1, nscal 835 f_id = ivarfl(isca(ii)) 836 call field_get_key_double(f_id, kvisl0, visls_0) 837 838 ! For scalars which are not variances, define the reference diffusivity 839 if (iscavr(ii).le.0 .and. visls_0.lt.-grand) then 840 call field_get_key_int(f_id, kscacp, iscacp) 841 if (iscacp.gt.0) then 842 ! For temperature, the diffusivity factor is directly the thermal conductivity 843 ! lambda = Cp * mu / Pr 844 ! where Pr is the (molecular) Prandtl number 845 visls_0 = viscl0 * cp0 846 else 847 visls_0 = viscl0 848 endif 849 call field_set_key_double(f_id, kvisl0, visls_0) 850 endif 851 852 ! For fluctuation variances, the diffusivity is that of the associated scalar. 853 iscal = iscavr(ii) 854 if (iscal.gt.0.and.iscal.le.nscal)then 855 call field_get_key_double(ivarfl(isca(iscal)), kvisl0, visls_0) 856 call field_get_key_double(f_id, kvisl0, visls_cmp) 857 call field_set_key_double(f_id, kvisl0, visls_0) 858 if (visls_cmp.gt.-grand) then 859 write(nfecra,1071) ii, iscal, ii, iscal, visls_0 860 endif 861 endif 862enddo 863 864! xyzp0 : reference point for hydrostatic pressure 865! The user should specify the 3 coordinates, otherwise 866! it is set to (0.,0.,0.). 867 868if (xyzp0(1).gt.-0.5d0*rinfin.and. & 869 xyzp0(2).gt.-0.5d0*rinfin.and. & 870 xyzp0(3).gt.-0.5d0*rinfin ) then 871 ixyzp0 = 1 872else 873 do ii = 1, 3 874 xyzp0(ii) = 0.d0 875 enddo 876endif 877 878! Turbulent fluxes constant for GGDH, AFM and DFM 879if (nscal.gt.0) then 880 do iscal = 1, nscal 881 882 call field_get_key_int(ivarfl(isca(iscal)), kturt, turb_flux_model) 883 turb_flux_model_type = turb_flux_model / 10 884 885 ! AFM and GGDH on the scalar 886 if (turb_flux_model_type.eq.1.or.turb_flux_model_type.eq.2) then 887 call field_get_key_struct_var_cal_opt(ivarfl(isca(iscal)), vcopt) 888 vcopt%idften = ANISOTROPIC_RIGHT_DIFFUSION 889 ctheta(iscal) = cthafm 890 call field_set_key_struct_var_cal_opt(ivarfl(isca(iscal)), vcopt) 891 ! DFM on the scalar 892 elseif (turb_flux_model_type.eq.3) then 893 call field_get_key_struct_var_cal_opt(ivarfl(isca(iscal)), vcopt) 894 vcopt%idifft = 0 895 vcopt%idften = ISOTROPIC_DIFFUSION 896 if (turb_flux_model.eq.31) then 897 ctheta(iscal) = cthebdfm 898 c2trit = 0.3d0 899 else 900 ctheta(iscal) = cthdfm 901 end if 902 call field_set_key_struct_var_cal_opt(ivarfl(isca(iscal)), vcopt) 903 ! GGDH on the thermal fluxes is automatically done 904 905 ! GGDH on the variance of the thermal scalar 906 do ii = 1, nscal 907 if (iscavr(ii).eq.iscal) then 908 call field_get_key_struct_var_cal_opt(ivarfl(isca(ii)), vcopt) 909 vcopt%idften = ANISOTROPIC_RIGHT_DIFFUSION 910 ctheta(ii) = csrij 911 call field_set_key_struct_var_cal_opt(ivarfl(isca(ii)), vcopt) 912 endif 913 enddo 914 else 915 ctheta(iscal) = csrij 916 endif 917 enddo 918endif 919 920! harmonic face viscosity interpolation 921if (imvisf.eq.1) then 922 do ivar = 1, nvar 923 call field_get_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 924 vcopt%imvisf = 1 925 call field_set_key_struct_var_cal_opt(ivarfl(ivar), vcopt) 926 enddo 927endif 928 929! VoF model enabled 930if (ivofmt.gt.0) then 931 ro0 = rho2 932 viscl0 = mu2 933 934 ! VOF algorithm: continuity of the flux across internal faces 935 call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt) 936 vcopt%imvisf = 1 937 call field_set_key_struct_var_cal_opt(ivarfl(ipr), vcopt) 938endif 939 940! Anisotropic diffusion/permeability for Darcy module 941if (ippmod(idarcy).eq.1) then 942 943 if (darcy_anisotropic_permeability.eq.1) then 944 call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt) 945 vcopt%idften = ANISOTROPIC_LEFT_DIFFUSION 946 call field_set_key_struct_var_cal_opt(ivarfl(ipr), vcopt) 947 endif 948 949 if (darcy_anisotropic_dispersion.eq.1) then 950 do iscal = 1, nscal 951 call field_get_key_struct_var_cal_opt(ivarfl(isca(iscal)), vcopt) 952 vcopt%idften = ANISOTROPIC_LEFT_DIFFUSION 953 call field_set_key_struct_var_cal_opt(ivarfl(isca(iscal)), vcopt) 954 enddo 955 endif 956 957 ! csrij = 1 and ctheta(iscal) = 1 for Darcy module 958 csrij = 1.d0 959 do iscal = 1, nscal 960 ctheta(iscal) = 1.d0 961 enddo 962 963 ! reference values for pressure and density 964 p0 = 0.d0 965 ro0 = 1.d0 966 967 ! be careful: if iturb was not initialized iturb is set to 0 to pass verini 968 if (iturb.eq.-999) iturb = 0 969 if (iturb.gt.0) then 970 write(nfecra,4001) 971 call csexit (1) 972 endif 973 974endif 975 976! Set iswdyn to 2 by default if not modified for pure diffusion equations 977do f_id = 0, nfld - 1 978 call field_get_type(f_id, f_type) 979 ! Is the field of type FIELD_VARIABLE? 980 if (iand(f_type, FIELD_VARIABLE).eq.FIELD_VARIABLE) then 981 call field_get_key_struct_var_cal_opt(f_id, vcopt) 982 if (vcopt%iswdyn.eq.-1.and. vcopt%iconv.eq.0) then 983 vcopt%iswdyn = 2 984 call field_set_key_struct_var_cal_opt(f_id, vcopt) 985 endif 986 endif 987enddo 988 989!=============================================================================== 990! 5. ELEMENTS DE albase 991!=============================================================================== 992 993if (iale.ge.1) then 994 if (isuite.eq.0 .and. italin.eq.-999 ) italin = 1 995else 996 italin = 0 997endif 998 999!=============================================================================== 1000! 6. COEFFICIENTS DE alstru 1001!=============================================================================== 1002 1003if (betnmk.lt.-0.5d0*grand) betnmk = (1.d0-alpnmk)**2/4.d0 1004if (gamnmk.lt.-0.5d0*grand) gamnmk = (1.d0-2.d0*alpnmk)/2.d0 1005if (aexxst.lt.-0.5d0*grand) aexxst = 0.5d0 1006if (bexxst.lt.-0.5d0*grand) bexxst = 0.0d0 1007if (cfopre.lt.-0.5d0*grand) cfopre = 2.0d0 1008 1009!=============================================================================== 1010! 7. PARAMETRES DE cplsat 1011!=============================================================================== 1012 1013! Get number of couplings 1014 1015call nbccpl(nbrcpl) 1016!========== 1017 1018if (nbrcpl.ge.1.and.iturbo.ne.0) then 1019 ifaccp = 1 1020endif 1021 1022!=============================================================================== 1023! 8. Define Min/Max clipping values of void fraction of VOF model 1024!=============================================================================== 1025 1026if (ivofmt.gt.0) then 1027 call field_get_key_double(ivarfl(ivolf2), kscmin, clvfmn) 1028 call field_get_key_double(ivarfl(ivolf2), kscmax, clvfmx) 1029 1030 if (clvfmn.lt.-0.5d0*grand) then 1031 clvfmn = 0.d0 1032 if (iand(ivofmt,VOF_MERKLE_MASS_TRANSFER).ne.0) clvfmn = epzero 1033 endif 1034 if (clvfmx.gt.0.5d0*grand) then 1035 clvfmx = 1.d0 1036 if (iand(ivofmt,VOF_MERKLE_MASS_TRANSFER).ne.0) clvfmx = 1.d0-epzero 1037 endif 1038 1039 call field_set_key_double(ivarfl(ivolf2), kscmin, clvfmn) 1040 call field_set_key_double(ivarfl(ivolf2), kscmax, clvfmx) 1041endif 1042 1043!=============================================================================== 1044! 9. STOP SI PB 1045!=============================================================================== 1046 1047if (iok.ne.0) then 1048 call csexit (1) 1049endif 1050 1051 1011 format( & 1052'@', /,& 1053'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1054'@', /,& 1055'@ @@ WARNING: ABORT IN THE DATA SPECIFICATION', /,& 1056'@ ========', /,& 1057'@ ',a6,' = ', i10, /,& 1058'@ ',a6,' WILL BE INITIALIZED AUTOMATICALLY', /,& 1059'@ DO NOT MODIFY IT.,' /,& 1060'@', /,& 1061'@ The calculation will not be run.', /,& 1062'@', /,& 1063'@ Check cs_user_parameters.f90', /,& 1064'@', /,& 1065'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1066'@', /) 1067 1021 format( & 1068'@', /,& 1069'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1070'@', /,& 1071'@ @@ WARNING: ABORT IN THE DATA SPECIFICATION', /,& 1072'@ ========', /,& 1073'@ SCALAR ', i10,' ',a6,' = ', i10, /,& 1074'@ ',a6,' WILL BE INITIALIZED AUTOMATICALLY', /,& 1075'@ DO NOT MODIFY IT.,' /,& 1076'@', /,& 1077'@ The calculation will not be run.', /,& 1078'@', /,& 1079'@ Check cs_user_parameters.f90', /,& 1080'@', /,& 1081'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1082'@', /) 1083 1131 format( & 1084'@', /,& 1085'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1086'@', /,& 1087'@ @@ WARNING: ADVANCED MODIFICATION FOR', /,& 1088'@ ========,' /,& 1089'@ ',a17,' OF THE VARIABLE' /,& 1090'@ ',a6,'.' /,& 1091'@', /,& 1092'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1093'@', /) 1094 1095 1061 format( & 1096'@', /,& 1097'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1098'@', /,& 1099'@ @@ WARNING: ABORT IN THE DATA SPECIFICATION', /,& 1100'@ ========', /,& 1101'@ ICALHY must be an integer equal to 0 or 1', /,& 1102'@', /,& 1103'@ Its value is ',i10, /,& 1104'@', /,& 1105'@ The calculation will not be run.', /,& 1106'@', /,& 1107'@ Check cs_user_parameters.f90', /,& 1108'@', /,& 1109'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1110'@', /) 1111 1071 format( & 1112'@', /,& 1113'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1114'@', /,& 1115'@ @@ WARNING: IN THE DATA SPECIFICATION', /,& 1116'@ ========', /,& 1117'@', /,& 1118'@ The scalar ',i10, ' is the fluctuations variance', /,& 1119'@ of the scalar ',i10, /,& 1120'@', /,& 1121'@ The diffusivity_ref value of the scalar ', i10, /,& 1122'@ must not be set:', /,& 1123'@ it is automatically set equal to the scalar', /,& 1124'@ diffusivity ', i10, ' i.e. ',e14.5, /,& 1125'@', /,& 1126'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1127'@', /) 1128 2000 format( & 1129'@', /,& 1130'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1131'@', /,& 1132'@ @@ WARNING: IN THE DATA SPECIFICATION', /,& 1133'@ ========', /,& 1134'@', /,& 1135'@ The k-omega turbulence model has been chosen. In order to', /,& 1136'@ have a correct calculation restart, the ICDPAR indicator',/,& 1137'@ has been set to 1 (read the wall distance in the restart',/,& 1138'@ file).', /,& 1139'@ If this initialization raises any issue (modification of,' /,& 1140'@ the number and position of the wall faces since the', /,& 1141'@ previous calcuation), force ICDPAR=1 (there might be,' /,& 1142'@ a small shift in the turbulent viscosity at the,' /,& 1143'@ first time-step).,' /,& 1144'@', /,& 1145'@ The calculation will be run.', /,& 1146'@', /,& 1147'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1148'@', /) 1149 2001 format( & 1150'@', /,& 1151'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1152'@', /,& 1153'@ @@ WARNING: IN THE DATA SPECIFICATION', /,& 1154'@ ========', /,& 1155'@', /,& 1156'@ The k-omega turbulence model has been chosen, with the,' /,& 1157'@ option for a re-calculation of the wall distance', /,& 1158'@ (ICDPAR=-1). There might be a small shift in the', /,& 1159'@ turbulent viscosity at the first time-step.', /,& 1160'@', /,& 1161'@ The calculation will be run.', /,& 1162'@', /,& 1163'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1164'@', /) 1165 3000 format( & 1166'@', /,& 1167'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1168'@', /,& 1169'@ @@ WARNING: IN THE DATA SPECIFICATION', /,& 1170'@ ========', /,& 1171'@', /,& 1172'@ The cavitation model requires an upwind convection scheme' ,/,& 1173'@ for the void fraction (BLENCV(IVOLF2)=',e14.5,').', /,& 1174'@ The user has set BLENCV(IVOLF2)=',e14.5, /,& 1175'@', /,& 1176'@ The upwind scheme for the void fraction is forced.', /,& 1177'@', /,& 1178'@ The calculation will be run.', /,& 1179'@', /,& 1180'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1181'@', /) 1182 4001 format( & 1183'@', /,& 1184'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1185'@', /,& 1186'@ @@ WARNING: IN THE DATA SPECIFICATION', /,& 1187'@ ========', /,& 1188'@', /,& 1189'@ A turbulence model can not be used with the' /,& 1190'@ gound water flows modeling.', /,& 1191'@', /,& 1192'@ The calculation will not be run.', /,& 1193'@', /,& 1194'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1195'@', /) 1196 5000 format( & 1197'@', /,& 1198'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1199'@', /,& 1200'@ @@ WARNING: IN THE DATA SPECIFICATION', /,& 1201'@ ========', /,& 1202'@', /,& 1203'@ The pseudo coupling of turbulent dissipation and turbulent',/,& 1204'@ kinetic energy (ikecou = 1) is not compatible with the use',/,& 1205'@ of fluid/solid option to disable the dynamic in the solid ',/,& 1206'@ cells (fluid_solid =1). ', /,& 1207'@', /,& 1208'@ The parameter ikecou is forced to 0 (no coupling)', /,& 1209'@', /,& 1210'@ The calculation will be run.', /,& 1211'@', /,& 1212'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',/,& 1213'@', /) 1214return 1215end subroutine 1216