$$$$ KEPSILON
* KEPSILON  PROCEDUR  MAGN      09/09/01    21:15:17     6478           
'DEBPROC' KEPSILON ;
ARGU RX*TABLE ;
************************************************************************
*      Ro UN Mu DT (GB T)
*
*         EN
*         KN
*       -> CHPOINT générés dans la table inco
*         TKTE      teta=k/epsilon
*         TETK      i   =epsilon/k     alias TKTE
*         NUTI      valeur intermédiaire de NUT
*         FI        inconnue Fi
*         PRODT     Prodution turbulente
*         TKTI      teta=k/epsilon intermédiaire
*         Ksi       facteur de déséquilibre : nut P / epsilon
*         MUF       viscosité dynamique effective
*                             (tourbillonnaire+moléculaire)
************************************************************************
rv=rx.'EQEX' ;
iarg=rx.'IARG' ;
*NASTOK = rv.'NAVISTOK' ;
 $mod=rx.'DOMZ' ;
 Dg=doma $mod 'XXDIAGSI' ;
*mess ' DEBUT KEPSILON';

* Lecture du 1er Argument la densité
Si(ega ('TYPE' rx.'ARG1') 'MOT     ');
Ro = rv.inco.(rx.'ARG1');
Sinon ;
Ro = rx.'ARG1';
Finsi ;

* Lecture du 3ème Argument la viscosité cinématique
Si(ega ('TYPE' rx.'ARG3') 'MOT     ');
Mu = rv.inco.(rx.'ARG3');
Sinon ;
Mu = rx.'ARG3';
Finsi ;
Mus2= Mu*(0.5);

Si (EGA (TYPE Ro) 'FLOTTANT');
 iRo = 1. / Ro ;
Sinon ;
 iRo = inve Ro;
Finsi ;

Si (EGA (TYPE Mu) 'FLOTTANT');
 Mum = Mu ;
 iMu = 1./MU ;
Sinon ;
 Mum = (Maxi Mu) + (Mini Mu) * 0.5;
 iMu = inve MU ;
Finsi ;

 Nu  = Mu * iRo;
 iNu = Ro * iMu;

Si (EGA (TYPE Nu) 'FLOTTANT');
 Num = Nu*0.5 ;
Sinon ;
 Num = (Maxi Nu) + (Mini Nu) * 0.25;
Finsi ;

Si( (non(ega rx.'KOPT'.'KIMPL' 1)) et
    (non(ega rx.'KOPT'.'KFORM' 1)) );
 mess ' KEPSILON KIMPL' (rx.'KOPT'.'KIMPL');
 mess ' KEPSILON KFORM' (rx.'KOPT'.'KFORM');
 mess 'Options non prevues IMPL obligatoire ' ;
 QUITTER KEPSILON ;
Finsi ;

******* Options du K-epsilon *******************************************
 EDPNUT=VRAI;Tmin=FAUX;
 EDPFI =FAUX;
 Kbw=FAUX ;KCnu=FAUX;RNG=FAUX;Filtre=FAUX;CSTE=FAUX;M2M=FAUX;Kimpr=FAUX;
 KRet=FAUX;KLbr=FAUX;Kchn=FAUX;Brls=FAUX;Kkl=FAUX;Kbrey=FAUX;
 Brjl=FAUX;Brlb=FAUX;
 V0=FAUX;TETA=VRAI;PERI=FAUX;

 Si ('EXIST' rv 'ALGO_KEPSILON');
  Si('EGA' (TYPE rv.'ALGO_KEPSILON') LISTMOTS);
    ko=0;
    lko=dime  rv.'ALGO_KEPSILON';
   Si ('EXIST' rv.'ALGO_KEPSILON' 'IMPR')   ; Kimpr=VRAI;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'V0')     ; V0=VRAI   ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'RNG')    ; RNG=VRAI  ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'CSTE')   ; CSTE=VRAI ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'M2M')    ; M2M=VRAI  ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'Bw')     ; Kbw=VRAI  ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'Cnu')    ; KCnu=VRAI ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'Filt') ;Filtre=VRAI  ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'Ret')    ; KRet=VRAI ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'KLbr')   ; KLbr=VRAI ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'Chie')   ; Kchn=VRAI ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'Shar')   ; Brls=VRAI ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'Jone')   ; Brjl=VRAI ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'Lam' )   ; Brlb=VRAI ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'Perio')  ; PERI=VRAI ;ko=ko+1; Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'Nut');EDPNUT=VRAI;EDPFI=FAUX;ko=ko+1;
   Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'Fi') ;EDPNUT=FAUX;EDPFI=VRAI;ko=ko+1;
   Finsi;
   Si ('EXIST' rv.'ALGO_KEPSILON' 'KL') ;EDPNUT=FAUX;EDPFI=FAUX;ko=ko+1;
   Kkl=VRAI; Finsi;

   Si (NON (EGA ko lko));
    Mess ' ' ;Mess ' ' ;
    Mess '************************************************************'
    '********';
    Mess '******* Il y a une option invalide pour le  modèle K-Epsilon'
    ' *******';
    Mess '******* ou bien elle apparait plusieurs fois                '
    ' *******';
    Mess '******* Liste de rv. ALGO_KEPSILON pour contrôle            '
    ' *******';Mess ' ';
    list rv.'ALGO_KEPSILON';
    Mess '************************************************************'
    '********';
    QUITTER KEPSILON ;
   Finsi ;
  Finsi ;
 Sinon ;
  Si Kimpr;
 mess '*********** Utilisation Standard du modèle K-Epsilon **********';
  Finsi ;
 Finsi ;


   KFORM=rx.'KOPT'.'KFORM' ;
   IDCEN=rx.'KOPT'.'IDCEN' ;
   CMD  =rx.'KOPT'.'CMD'   ;
   IDIV =rx.'KOPT'.'IDIV'  ;

   rx.'KOPT'.'IDIV'  =0  ;
   rx.'KOPT'.'IDCEN' =2  ;
*  rx.'KOPT'.'CMD'   =0.5;
   rx.'KOPT'.'CMD'   =0.2;
   NBTETA=10;

 Si V0 ;
   rx.'KOPT'.'CMD' =0.2 ; EDPNUT=VRAI;EDPFI=FAUX;
   Tmin=FAUX;Kbw=FAUX;KCnu=FAUX;M2M=FAUX;CSTE=FAUX;
   Filtre=FAUX;Kret=FAUX;KLbr=FAUX;Kchn=FAUX;Brls=FAUX;Brjl=FAUX;
   Brlb=FAUX;PERI =FAUX;
   TETA=VRAI
 Finsi ;

 Si(Kkl ou KLbr ou Kchn ou Brls ou Brjl ou Brlb);
 Kbrey= VRAI;
 Tmin = FAUX;
 Finsi;

Si Kimpr ;
   lopt = chai ' ' ;
   Si V0     ; lopt = chai lopt ' V0 '    ; Finsi ;
   Si EDPFI  ; lopt = chai lopt ' FI '    ; Finsi ;
   Si EDPNUT ; lopt = chai lopt ' NUT '   ; Finsi ;
   Si RNG    ; lopt = chai lopt ' RNG '   ; Finsi ;
   Si CSTE   ; lopt = chai lopt ' CSTE '  ; Finsi ;
   Si M2M    ; lopt = chai lopt ' M2M  '  ; Finsi ;
   Si Kbw    ; lopt = chai lopt ' Bw '    ; Finsi ;
   Si KCnu   ; lopt = chai lopt ' Cnu '   ; Finsi ;
   Si KRet   ; lopt = chai lopt ' Ret '   ; Finsi ;
   Si Kkl    ; lopt = chai lopt ' KL    ' ; Finsi ;
   Si KLbr   ; lopt = chai lopt ' KLbr'   ; Finsi ;
   Si Kchn   ; lopt = chai lopt ' Chien ' ; Finsi ;
   Si Brls   ; lopt = chai lopt ' Sharma' ; Finsi ;
   Si Brjl   ; lopt = chai lopt ' Jones ' ; Finsi ;
   Si Brlb   ; lopt = chai lopt ' Lam   ' ; Finsi ;
   Si Peri   ; lopt = chai lopt ' Peri  ' ; Finsi ;
   Si Filtre ;
        lopt = chai lopt ' Filtre (URANS) '      ;
                                            Finsi ;
 Si ('EXIST' rv 'ALGO_KEPSILON');
 mess '* K-epsilon : Options : ' lopt ;
 Finsi ;
Finsi ;

   IPT=rv.'PASDETPS'.'NUPASDT';

******* Options du K-epsilon * Fin *************************************

************** Constantes diverses **********
* Constante de Bradshaw Bw  de Smagorinski Cs
Bw=0.3 ;
Cs=0.1 ;
Ka=0.41;
* Paramètre de Menter Nut P < Ksim * epsilon (ksim=10)
ksim=10.;

************** Constantes du K epsilon ****************
cnuo=0.09;
c1=1.44;
c2=1.92;
sgk=1. ;
sge=1.3;

Si (Kkl ou KLbr);
****** Constantes du modèle K-L ou K-L Bas Reynolds ***
cnuo=0.09;
c1=1.44;
cl=1./Ka;
cd=0.9;
Aetmu=63.;
Aeteps=3.8;
sgk=0.9;
finsi;

Si Kchn;
****** Constantes K epsilon Bas Reynolds Chien ********
cnuo=0.09;
c1=1.35;
c2=1.8 ;
sgk=1. ;
sge=1.3;
finsi;

Si RNG ;
************** Constantes du RNG K epsilon ************
cnuo=0.0845;
c1=1.42;
c2=1.68;
sgk=0.7179;
sge=0.7179;
etao=4.377;
beta=0.012;
etai=1./etao;
finsi;

Si M2M ;
************** Constantes de Mohammadi Medic **********
cnuo=0.09;
c1=0.1296;
c2=11./6.;
sgk=1.;
sge=Ka*Ka*(cnuo**0.5)/(c2*cnuo - c1);
finsi;

Si CSTE;
************** Constantes lues dans rv.inco ***********
cnuo=rv.inco.'cnu';
c2=rv.inco.'c2';
sgk=rv.inco.'sgk';
sge=rv.inco.'sge';
c1=c2 - (Ka*Ka/sge/(cnuo**0.5)) ;
finsi;

Si Kimpr;
 mess '* K-epsilon, Constantes: cnuo=' cnuo ' c1=' c1 ' c2=' c2;
 mess '* K-epsilon, Constantes:  sgk=' sgk ' sge=' sge;
Finsi;

rcnu = cnuo**0.5 ;

*    :  mp=-9. ; mq = 5. ;
*    : mp= 3. ; mq =-2. ;
*    : mp=-17.; mq = 9. ;
*    : mp=(-1.)*c1; mq = 1. ;
*    : mp=2.      ; mq =-1. ;
*    : mp=-2.     ; mq = 1. ;
*    : mp=(-1.)*c2; mq = 1. ;
* Fi : mp=-3. ; mq = 2. ;

**************** INITIALISATIONs ***********************
*TYPINV = 4 ;  cette option est prematureeé
TYPINV = 3 ;

Si(non (exist rx 'rxk'));

rxt1=table 'KIZX';
rxt2=table 'KIZX';
rxk=table 'KIZX';
rxe=table 'KIZX';
rx.'rxt1'=rxt1;
rx.'rxti'=rxti;
rx.'rxt2'=rxt2;
rx.'rxk'=rxk;
rx.'rxe'=rxe;

mtinv = (eqex) . 'METHINV' ; 
rv.'METH_KEP'=mtinv;
 mtinv.'TYPINV'=TYPINV;
mtinv.'IMPINV'=1;
mtinv.'TYRENU'='SLOA';
mtinv.'PCMLAG'='APR2';
mtinv.'OUBMAT' = 0 ;
mtinv.'SCALING' = 0 ;
mtinv.'NITMAX'=500;
mtinv.'RESID'=1.e-12;
mtinv.'PRECOND'=3;
mtinv.'FCPRECT'=1;
mtinv.'FCPRECI'=1;
mtinv.'BCGSBTOL'=1.D-40;

*********** étape de diffusion de K *********************
rxk.'EQEX'=rx.'EQEX';
rxk.'NOMZONE'=' ';
rxk.'DOMZ'=rx.'DOMZ';
rxk.'TDOMZ'=0;
rxk.'NOMOPER'=mot 'LAPN';
rxk.'KOPT'=rx.'KOPT';
rxk.'IARG'=1 ;
rxk.'ARG1'=1.e-10;

*********** étape de diffusion de Epsilon ***************
rxe.'EQEX'=rx.'EQEX';
rxe.'NOMZONE'=' ';
rxe.'DOMZ'=rx.'DOMZ';
rxe.'TDOMZ'=0;
rxe.'NOMOPER'=mot 'LAPN';
rxe.'KOPT'=rx.'KOPT';
rxe.'IARG'=1 ;
rxe.'ARG1'=1.e-10;

Si (NON (EGA (dime rx.'LISTINCO') 2));
  mess ' Erreur dans la procédure Kepsilon !';
  mess ' Il doit y avoir deux inconnues !';
  quitter Kepsilon;
Finsi;
MI1=EXTR rx.'LISTINCO' 1 ;
MI2=EXTR rx.'LISTINCO' 2 ;
rxk.'LISTINCO'=MOTS MI1 ;
rxe.'LISTINCO'=MOTS MI2 ;
Ko=(abs rv.inco.MI1) + 1.e-10;
Eo=(abs rv.inco.MI2) + 1.e-10;
To=Ko*(inve Eo);
rv.inco.'TKTE' =To ;

Fio=inve (To*To*Ko);
rv.inco.'Fio'  =Fio;

************************ dT/dt **************************

rxt1.'EQEX'=rx.'EQEX';
rxt1.'NOMZONE'=' ';
rxt1.'DOMZ'=rx.'DOMZ';
rxt1.'TDOMZ'=0;
rxt1.'NOMOPER'=mot 'TSCA_TK';
rxt1.'KOPT'=rx.'KOPT';
rxt1.'LISTINCO'=MOTS 'TKTE' ;
rxt1.'IARG'=3;
*rxt1.'ARG1'=Nu ;
rxt1.'ARG2'=rx.'ARG2';
rxt1.'ARG3'=0. ;

rxt2.'EQEX'=rx.'EQEX';
rxt2.'NOMZONE'=' ';
rxt2.'DOMZ'=rx.'DOMZ';
rxt2.'TDOMZ'=0;
rxt2.'NOMOPER'=mot 'TSCA_TK';
rxt2.'KOPT'=rx.'KOPT';
*rxt2.'LISTINCO'=MOTS 'NUTI' ;
rxt2.'LISTINCO'=MOTS 'NUTI' ;
rxt2.'IARG'=3;
*rxt2.'ARG1'=Nu;
rxt2.'ARG2'=rx.'ARG2';
rxt2.'ARG3'=0. ;
 Si (NON ('EXIST' (rv.'INCO') 'NUTI'));
  rv.inco.'NUTI'=kcht $mod scal sommet Nu;
 Finsi ;
 Si (NON ('EXIST' (rv.'INCO') 'FI'));
  rv.inco.'FI'=kcht $mod scal sommet 0.;
 Finsi ;
Sinon;

rxk=rx.'rxk';
mtinv=rv.'METH_KEP';
rxe=rx.'rxe';
rxt1=rx.'rxt1';
rxt2=rx.'rxt2';

******** On crée les matrices de périodicité ************
Si Peri;
*st1 mat1 = KOPS 'MATRIK' ;
*----------------------------------------
 mr1=rv.'INCO'.'M1PERIODIC';
 mr2=rv.'INCO'.'M2PERIODIC';
 mati=rela 'UX' mr1 - 'UX' mr2;
 sti =depi mati 0.;
 mati=kops 'RIMA' mati;
 rxk.'PERIODIC_mat'=mati;
 rxk.'PERIODIC_st' =sti;
*----------------------------------------
Finsi;
**** Fin On crée les matrices de périodicité ************
Finsi ;

********************** Fin Initialisations *******************

 Si('EXIST' rxk 'PERIODIC_mat');
  matipc=rxk.'PERIODIC_mat';
  stipc =rxk.'PERIODIC_st' ;
 Finsi;

 MI1=EXTR rxk.'LISTINCO' 1 ;
 MI2=EXTR rxe.'LISTINCO' 1 ;

MEN=FAUX ;
 Si (exist rv.'CLIM' MI1);
  climk=nomc MI1 (exco MI1 rv.clim) ;
  climk=((abs climk) + 1.e-10);
    Si(exist rv 'climk');
    xo=nomc MI1 (exco MI1 rv.'CLIM');
    x1=redu rv.'climk' (extr xo MAILLAGE) ;
    x2=rv.'climk' - x1;
    climk=x2 + xo;
    Finsi ;
  MEN=VRAI;
 Sinon ;
    Si(exist rv 'climk');
    climk=abs (rv.'climk');
    Sinon;
  mess ' Erreur dans la procédure Kepsilon !';
  mess ' Une condition limite sur K est nécessaire !';
  quitter Kepsilon;
    climk a = kops 'MATRIK' ;
    Finsi ;
 Finsi ;

 Si (exist rv.'CLIM' MI2);
  clime=nomc MI2 (exco MI2 rv.clim) ;
  clime=((abs clime) + 1.e-10);
  MEN=VRAI;
    Si(exist rv 'clime');
    xo=nomc MI2 (exco MI2 rv.'CLIM');
    x1=redu rv.'clime' (extr xo MAILLAGE) ;
    x2=rv.'clime' - x1;
    clime=x2 + xo;
    Finsi ;
 Sinon ;
    Si(exist rv 'clime');
    clime=abs (rv.'clime');
    Sinon;
  mess ' Erreur dans la procédure Kepsilon !';
  mess ' Une condition limite sur Epsilon est nécessaire !';
  quitter Kepsilon;
    clime a = kops 'MATRIK' ;
    Finsi ;
 Finsi ;

  Si MEN ;
 lic=extr rv.clim 'COMP' ;
 nbic=dime lic ;
 climr a = kops 'MATRIK' ;

 repeter bloc1 nbic ;
 mic=extr &bloc1 lic ;
 Si((non (ega mic MI1)) et (non (ega mic MI2)));
 climr = climr et (nomc mic (exco mic (rv.clim)));
 Finsi ;
 Fin Bloc1 ;
 rv.'CLIM'=climr ;
 lic=extr climr 'COMP' ;
  Finsi;

rv.'climk'=climk;
rv.'clime'=clime;

 dtn=rx.'ARG4';
 tdt=type dtn;
 alfa=rv.'ALFA';
 si (ega tdt 'MOT     ');
* dtn = 'MOT' ('TEXTE' ('CHAINE' dtn)) ;
  si(ega ('MOT' dtn) ('MOT' 'DELTAT'));
  dtn=rv.'PASDETPS'.'DELTAT';
  sinon ;
  dtn=rv.'INCO'.dtn;
  finsi;
 finsi;
  dtn=dtn*alfa;
* mess ' Kepsilon : Pas de temps ' dtn ;

  dt=dtn ;
  dti=1./dt;
  mdt=(-1.)*dt;

*mess 'On calcule les grandeurs initiales Ko Eo To ';
 Ko=(abs rv.inco.MI1) + 1.e-10 ;
 Eo=(abs rv.inco.MI2) + 1.e-10 ;

 mtinv.'XINIT'=Ko;
 To=kcht $mod scal sommet (Ko*(inve Eo));
 rv.inco.'TKTE'=To ;
 Fio=inve (To*To*Ko);

 Io=kcht $mod scal sommet (Eo*(inve Ko));
 rv.inco.'TETK'=Io ;

* Muto est repris dans inco 
 Si(exist rv.inco 'Muto');
 Muto=rv.inco.'Muto';
 Muto=Muto + 1.e-10 ;
 Sinon;
 Muto=Ro * cnuo * Ko * To ;
 Finsi;
 Nuto = Muto*iRo;

******************** Calcul de P ***********************
* mess 'AVANT CALCUL P';
MIU=rx.'ARG2' ;
UN=rv.inco.miu;

Si(ega iarg 4);
P = PRODT UN  $mod ;
Sinon ;
 Si(ega iarg 6);
 si(ega ('TYPE' rx.'ARG5') 'MOT     ');
 m4=rv.inco.(rx.'ARG5');
 sinon;
 m4=rx.'ARG5';
 finsi;
 si(ega ('TYPE' rx.'ARG6') 'MOT     ');
 m5=rv.inco.(rx.'ARG6');
 sinon;
 m5=rx.'ARG6';
 finsi;
 P= PRODT UN $mod m4 m5 ;
 Sinon;
 mess ' Nombre d argument incorrect dans KEPSILON ';
 erreur KEPSILON ;
 Finsi;
Finsi;

 lc = extr UN 'COMP';
 mdu2=  PSCA UN lc UN lc         ;
 mdu2= mdu2 + 1.e-10;

Si (non Kbrey);
* Pmax
* Réalisabilité sur P dés/activée [A]  (P < 10 epsilon/nut (Menter))
 Pmax =(ksim*Eo)*(inve Muto);
 a =Pmax    ; al=0.7 ; ala=al*a ;
*                               b=ala*al*((1.-al)**(-1.));
             b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
*b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
 ik=masq P 'INFERIEUR' ala ;
*P =(ik*P)+((1.-ik)*a*(P+b)*(inve(a+P+b)));
 P =(ik*P)+((1.-ik)*a*((P*P+b2)*(inve(a+(P*P+b2)))));
*P =(ik*P)+((1.-ik)*a*((P*P*P+b3)*(inve(a+(P*P*P+b3)))));
* Pmax
Finsi;
 
 rv.inco.'PRODT'= P ;
*mess ' P    = ' (maxi P   ) (mini P   );

 rP=(ABS(P+1.e-20))**0.5 ;

******************  FIN CALCUL P  **********************
********************************************************


 Cnu = kcht $mod scal sommet (cnuo);
 Ak = Cnu * To * To * P ;

 Si KCnu ;
* On recalcule Cnu *************************************
* alternative hyperbolique
 mess ' déséquilibre Ak ' (maxi Ak) (mini Ak) ;
 alfa = 0.103 ; beta = 0.15 ;
 Cnu = alfa*(inve (Ak + beta));
* alternative exponentielle Alfa = Log (cnu(0)/cnu(1))
* Ak =0 -> cnu=0.7   Ak = 1 cnu=0.09 -> alfa = 2.0512
* Cnu = cnuo*( (exp (2.0512*(1.-Ak))) + 1.e-10);
* On recalcule Cnu Fin *********************************
********************************************************
 Finsi ;

**************** Calcul constante c2 *******************

*f3=1.;
 ATETA=(c1-1.); BTETA=c2-1.;

Si RNG;
eta=rP*To;
eta3=eta*eta*eta;
cp2=(cnu*eta3*(1.-(etai*eta))) *
 (inve(1.+(beta*eta3))) ;
*cc2=c2 + cp2;
 c2=c2 + cp2;
oubli cp2; oubli eta; oubli eta3;
finsi;

 cnuo=0.09; c1=1.44; c2=1.92; sgk=1.; sge=1.3; alf=1.;
 Ret = Ro * Ko * To * iMu                         ;
 fmu = 1.                                         ;
 Rocnu=Ro*Cnu;
 Epsok=0.;
 Epsoe=0.;

**************** Modèle K-L ****************************

 Si Kkl;
 cnuo=0.09; c1=1.44; c2=1.92; sgk=1.; sge=1.3; alf=1.;
 fmu = 1. ;
 Cnu=Cnuo                                  ;
 Rocnu=Ro*Cnuo                             ;
 Ret = Muto * iMu * (1./Cnu);
 Finsi ;

**************** Modèle K-L Bas-Reynolds ***************

 Si KLbr;      
 cnuo=0.09; c1=1.44; c2=1.92; sgk=1.; sge=1.3; alf=1.;
 fmu = 1. ;
 Cnu=Cnuo ;
 Rocnu=Ro*Cnu                              ;
 dparoi=rv.'INCO'.'dparoi';
 yet=dparoi*(Ko**0.5)*(iNu);
 lmu=cl*dparoi*(1. - (exp ( (-1./Aetmu)*yet)));
 leps=cl*dparoi*(1. - (exp ( (-1./Aeteps)*yet)));
 leps=leps+1.e-10;
 Eo=(Ko**1.5)*(inve leps);
 ATETA=lmu*(inve leps)/2.; BTETA=0.5;
*f3=lmu*(inve leps);
 Finsi ;

**************** Modèle Bas-Reynolds de Chien **********

 Si Kchn;
 cnuo=0.09; c1=1.35; c2=1.8; sgk=1.; sge=1.3;
 yplus=rv.'INCO'.'yplus';
 dparoi=rv.'INCO'.'dparoi';
 fmu = 1. - (exp (-0.0115*yplus));
 Cnu=Cnuo*fmu;
 Rocnu= Ro*cnu;
 Ret = Ro*Ko*Ko*iMu*(inve Eo);
 f2  = 1. - (0.22 * (exp((-1./36.)*Ret*Ret)));
*E=(-2.)*Mu*Eo*(dparoi**(-2.))*(exp ((-0.5)*yplus));
*epso=2.*Mu*Ko*(dparoi**(-2.));
*alf=1. + ((inve Eo)*(E*To + epso));

 Epsoe=2.*Mu*(dparoi**(-2.))*(exp ((-0.5)*yplus));
 epsok=2.*Mu*(dparoi**(-2.));
 alf=1.;

 ATETA=(c1-1.); 
 BTETA = (c2*f2) - alf;

 Finsi ;

**************** Modèle Bas-Reynolds de Launder-Sharma *

 Si Brls;
 lux=mots 'UX' 'UY';
 si(EGA (VALE DIME) 3); lux=lux et (mots 'UZ');Finsi;
 cnuo=0.09; c1=1.44; c2=1.92; sgk=1.; sge=1.3;
 Ret = Ro*Ko*Ko*iMu*(inve Eo);
 fmu = exp(-3.5*(inve((1.+ (Ret*(1./50.)))**2.)));
 Cnu = Cnuo*fmu;
 Rocnu=Ro*Cnu  ;
 rKo=Ko**0.5;
 grk=kops 'GRADS' rKo $mod;
 epsok=2.*Mu*(psca grk lux grk lux)*(inve Ko);
  Grp=Kops 'GRADS' P $mod;
  D2u2=1./4.*(psca Grp lux Grp lux)*(inve P);

 E=2.*mu*Muto*D2u2;
 f2  = 1. - (0.3 * (exp((-1.)*Ret*Ret)));
 alf=1. + ((inve Eo)*E*To);

 ATETA=(c1-1.); 
 BTETA = (c2*f2) - alf;

 Finsi ;

**************** Modèle Bas-Reynolds de Jones-Launder **

 Si Brjl;
 lux=mots 'UX' 'UY';
 si(EGA (VALE DIME) 3); lux=lux et (mots 'UZ');Finsi;
 cnuo=0.09; c1=1.55; c2=2.; sgk=1.; sge=1.3;
 Ret = Ro*Ko*Ko*iMu*(inve Eo);
 fmu = exp(-2.5*(inve(1.+ (Ret*(1./50.)))));
 Cnu = Cnuo*fmu;
 Rocnu=Ro*Cnu  ;
 rKo=Ko**0.5;
 grk=kops 'GRADS' rKo $mod;
 epso=2.*Mu*(psca grk lux grk lux);
  Grp=Kops 'GRADS' P $mod;
  D2u2=1./4.*(psca Grp lux Grp lux)*(inve P);
*Grp=Kops 'GRADS' (exco un 'UX') $mod;
*Grp=Kops 'GRADS' (exco grp 'UY') $mod;
*Grp=exco grp 'UY';
*mess 'Grp ' (mini grp) (maxi grp);
*D2u2=Grp*Grp;

 E=2.*mu*Muto*D2u2;
 f2  = 1. - (0.3 * (exp((-1.)*Ret*Ret)));
 alf=1. + ((inve Eo)*(E*To + epso));

 ATETA=(c1-1.); 
 BTETA = (c2*f2) - alf;

 Finsi ;

**************** Modèle Bas-Reynolds de Lam-Bremhorst **

 Si Brlb;
 cnuo=0.09; c1=1.44; c2=1.92; sgk=1.; sge=1.3;
 Ry = (Ko**0.5)*dparoi*iMu*Ro;
 Ret = Muto * iMu * (1./cnuo);
 fmu = ((1.-(exp(-0.0165*Ry)))**2.)*(1.+(20.5*(inve Ret)));
 f1=1.+((0.05*(inve fmu))**3.);
 Cnu = Cnuo*fmu;
 Rocnu=Ro*Cnu  ;
 f2  = 1. - (exp((-1.)*Ret*Ret));
 c1 = c1 * f1 ;
 c2 = c2 * f2 ;
 E=0.;
 epso=0.;
 alf=1. + ((inve Eo)*(E*To + epso));

 ATETA=(c1-1.); 
 BTETA = (c2*f2) - alf;

 Finsi ;

***** Solution Teta(t) et Equation sur Nuti ************

Si TETA;
*................ EDP TETA .............................................

*TSCA
*mess ' TSCA sur TETA' ;
 rxt1.'NOMOPER'=mot 'TSCA_T';
 rxt1.'LISTINCO'=MOTS 'TKTE' ;
 rxt1.'IARG'=3;
 rxt1.'ARG1'=Num;
 rxt1.'ARG2'=UN       ;
 rxt1.'ARG3'=0. ;
 st  mat= TSCA rxt1 ;

 climt= nomc 'TKTE'
 ((exco MI1 (rv.'climk'))*(inve (exco MI2 (rv.'clime'))));

 T=To ;
 alt = 1.e-1;
 alt = 0.   ;
 REPETER BCLTETA NBTETA;

    Stp = BTETA + (To/dt) + (cnu*ATETA*P*T*T);
    Stn = (2.*cnu*ATETA*P*T) + (1./dt) ;

*MDIA
*mess ' MDIA sur TETA' ;
 rxt1.'NOMOPER'=mot 'MDIA_T';
 rxt1.'LISTINCO'=MOTS 'TKTE' ;
 rxt1.'IARG'=1;
 rxt1.'ARG1'=kcht $mod scal sommet (Stn+alt) ;
 st1 mat1 = MDIA rxt1 ;
 St1 = kcht $mod scal sommet comp 'TKTE'
                     (nomc 'TKTE'  (Dg*(Stp+(alt*T))));
 St1 = St et st1 ;
 mat1 = mat et mat1 ;

 mtinv.'MATASS'=mat1;
 mtinv.'MAPREC'=mat1;
 mtinv.'XINIT'=nomc 'TKTE'  T;
 mtinv.'IMPINV'=0;
 mtinv.'TYPINV'=TYPINV;

 Si PERI ;
  MPI1='TKTE';
  mat1=mat1 et (kops 'CHANINCO' mati
       ('MOTS' 'LX' 'UX')
       ('MOTS' (chai 'LX' MPI1) MPI1)
       ('MOTS' 'FLX' 'FX')
       ('MOTS' (chai 'LX' MPI1) MPI1) );
  st1 =st1 et (exco sti 'FLX' (chai 'LX' MPI1));
 Finsi;

*mess 'Equation sur TKTE';
 T1='KRES' mat1 'TYPI' mtinv 'CLIM' climt 'SMBR' st1 'IMPR' 0 ;
 T1= nomc 'SCAL' T1;

Si Tmin;
* Réalisabilité sur TETA dés/activée [A]
 iTetamin = (1./Num)*mdu2;
*iTetamin = (1./(Num*0.01))*mdu2;
 a =iTetamin ; al=0.3; ala=al*a ;
                                b=al*ala*((1.-al)**(-1.));
*            b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
*b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
 iT1 = inve T1;
 ik=masq iT1 'INFERIEUR' ala ;
 iT1=(ik*iT1)+((1.-ik)*a*(iT1+b)*(inve(a+iT1+b)));
*iT1=(ik*iT1)+((1.-ik)*a*((iT1*iT1+b2)*(inve(a+(iT1*iT1+b2)))));
*iT1=(ik*iT1)+((1.-ik)*a*((iT1*iT1*iT1+b3)*(inve(a+(iT1*iT1*iT1+b3)))));
 T1 = inve iT1 ;
Finsi ;

Si (non Kbrey);
* Tmax
* Réalisabilité sur T1 dés/activée [A]  (nut P < 10 epsilon (Menter))
*Tmax =(((ABS P)**0.5)*cnuo) inve ;
 Tmax =(inve(cnuo/ksim*(ABS P)))**0.5;
 a =Tmax    ; al=0.7 ; ala=al*a ;
*                               b=ala*al*((1.-al)**(-1.));
             b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
*b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
 ik=masq T1 'INFERIEUR' ala ;
*iE1=(ik*iE1)+((1.-ik)*a*(iE1+b)*(inve(a+iE1+b)));
 T1=(ik*T1)+((1.-ik)*a*((T1*T1+b2)*(inve(a+(T1*T1+b2)))));
*iE1=(ik*iE1)+((1.-ik)*a*((iE1*iE1*iE1+b3)*(inve(a+(iE1*iE1*iE1+b3)))));
* Tmax
Finsi;

 T1= (ABS T1) + 1.e-10 ;
 ERRT= (maxi (ABS (T1 - T)))/(maxi T1);
*mess ' iter =' &BCLTETA ' ERRT='  errt ;
 T = T1 ;

 Si (ERRT < 1.e-12) ; QUITTER BCLTETA ; Finsi ;
 Si ('EGA' &BCLTETA 12);
 mess ' Non convergence sur TETA en 10 itérations: ERRT=' errt ;
 Finsi ;
 FIN BCLTETA ;
*mess ' Iteration ' &BCLTETA ' ERRT=' errt ;
*............ Fin EDP TETA .............................................
Finsi;


*................ EDP Nut ..............................................
Si EDPNUT ;
Nutj  =     Nuto;

  Snp= cnu*T*P*(2.-c1)        ;
  Snn= (inve T)*((2.*alf)-c2) ;

*MDIA
*mess ' MDIA sur NUTI' ;
 rxt2.'NOMOPER'=mot 'MDIA_N';
 rxt2.'LISTINCO'=MOTS 'NUTI' ;
 rxt2.'IARG'=1;
 rxt2.'ARG1'=kcht $mod scal sommet ((1./dt) + Snn);
 st2 mat2 = MDIA rxt2 ;

 St2 = kcht $mod scal sommet comp 'NUTI'
            (nomc 'NUTI'  (Dg*Nuto*((1./dt) + Snp)));
*           (nomc 'NUTI'  (Dg*((Nuto*(1./dt)) + (Nutj*Snp))));

*TSCA
*mess ' TSCA sur NUTI' ;
 rxt2.'NOMOPER'=mot 'TSCA_N';
 rxt2.'LISTINCO'=MOTS 'NUTI' ;
 rxt2.'IARG'=4;
 rxt2.'ARG1'=1.;
 rxt2.'ARG2'=UN       ;
 rxt2.'ARG3'=(Mus2*iRo);
 rxt2.'ARG4'=0. ;
 st  mat = TSCA rxt2 ;
 St2= St et st2 ;
 mat2= mat et mat2 ;

 climn = nomc 'NUTI'
(cnu*((abs(nomc climk 'SCAL'))**2.)*(inve(abs (nomc clime 'SCAL'))));
 Nutv=Nutj;

 mtinv.'MATASS'=mat2;
 mtinv.'MAPREC'=mat2;
 mtinv.'XINIT'=nomc 'NUTI'  Nutv;
 mtinv.'IMPINV'=0;
 mtinv.'TYPINV'=TYPINV;

 Si PERI ;
  MPI1='NUTI';
  mat2=mat2 et (kops 'CHANINCO' mati
       ('MOTS' 'LX' 'UX')
       ('MOTS' (chai 'LX' MPI1) MPI1)
       ('MOTS' 'FLX' 'FX')
       ('MOTS' (chai 'LX' MPI1) MPI1) );
  st2 =st2 et (exco sti 'FLX' (chai 'LX' MPI1));
 Finsi;

*mess 'Equation sur NUTI';
 Nuti='KRES' mat2 'TYPI' mtinv 'CLIM' climn 'SMBR' st2 'IMPR' 0 ;
 Nuti=nomc 'SCAL' Nuti;
 Nuti = (masq Nuti 'SUPERIEUR' 0.)*Nuti + 1.e-10;

 Ki = kcht $mod scal sommet (Nuti * (inve(cnu * T)));

Si(NON Kbrey);
* Réalisabilité sur Ki dés/activée [A]
 a = 0.5*mdu2; al=0.9 ; ala=al*a ;
*                               b=ala*al*((1.-al)**(-1.));
*            b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
 b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
 ik=masq Ki 'INFERIEUR' ala ;
*Ki=(ik*Ki)+((1.-ik)*a*(Ki+b)*(inve(a+Ki+b)));
*Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki+b2)*(inve(a+(Ki*Ki+b2)))));
 Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki*Ki+b3)*(inve(a+(Ki*Ki*Ki+b3)))));
Finsi;

 Ki = kcht $mod scal sommet (Nuti * (inve(cnu * T)));
 Ei = kcht $mod scal sommet (Nuti * (inve(cnu * T * T)));

Finsi ;
*............ Fin EDP Nut ..............................................

*................ EDP Fi ...............................................
Si EDPFI  ;

 mp=-3. ; mq = 2. ;

 Sfn= ((inve To)*((mq*c2) + mp)) - (cnu*To*P*(mp + (mq*c1)));
 Sfp= 0.;

*Sfn= (cnu*To*P*(mp + (mq*c1)));
*Sfp= ((inve To)*((mq*c2) + mp)) ;

*MDIA
*mess ' MDIA sur FI' ;
 rxt2.'NOMOPER'=mot 'MDIA_F';
 rxt2.'LISTINCO'=MOTS 'FI' ;
 rxt2.'IARG'=1;
 rxt2.'ARG1'=kcht $mod scal sommet Sfn ;
 st2 mat2 = MDIA rxt2 ;
 St2 = kcht $mod scal sommet comp 'FI'
                     (nomc 'FI'  (Dg*Sfp));

*TSCA
*mess ' TSCA sur FI' ;
 rxt2.'NOMOPER'=mot 'TSCA_F';
 rxt2.'LISTINCO'=MOTS 'FI'   ;
 rxt2.'IARG'=3;
 rxt2.'ARG1'=Num;
 rxt2.'ARG2'=UN    ;
 rxt2.'ARG3'=0. ;
 st  mat= TSCA rxt2 ;
 St2 = St et st2 ;
 mat2 = mat et mat2 ;

 climf= nomc 'FI'
 ((exco MI1 (rv.'climk'))**mp)*((exco MI2 (rv.'clime'))**mq);

 mtinv.'MATASS'=mat2;
 mtinv.'MAPREC'=mat2;
 mtinv.'XINIT'=nomc 'FI'  Fio;
 mtinv.'IMPINV'=0;
 mtinv.'TYPINV'=TYPINV;
*mess 'Equation sur FI';
 Fi='KRES' mat2 'TYPI' mtinv 'CLIM' climf 'SMBR' st2 'IMPR' 0 ;
 Fi= nomc 'SCAL' Fi ;
 Fi= (ABS Fi) + 1.e-10 ;

 Ki  =kcht $mod scal sommet ((Fi**(1./(mp+mq)))*(T**(mq/(mp+mq))));

Si(NON Kbrey);
* Réalisabilité sur Ki dés/activée [A]
 a = 0.5*mdu2; al=0.9 ; ala=al*a ;
*                               b=ala*al*((1.-al)**(-1.));
*            b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
 b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
 ik=masq Ki 'INFERIEUR' ala ;
*Ki=(ik*Ki)+((1.-ik)*a*(Ki+b)*(inve(a+Ki+b)));
*Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki+b2)*(inve(a+(Ki*Ki+b2)))));
 Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki*Ki+b3)*(inve(a+(Ki*Ki*Ki+b3)))));
Finsi;

 Ei  =kcht $mod scal sommet (Ki*(inve  T ));
 Nuti=kcht $mod scal sommet (cnu*Ki*Ki*(inve  Ei));

Finsi ;
*............ Fin EDP Fi ...............................................

*................ modèle K-L ...........................................
 Si Kkl;
 Echl= rv.inco.'Echl';
 T2=T**(-2.);
 Ki  =kcht $mod scal sommet ((Echl**2.)*(T2));
** Réalisabilité sur Ki dés/activée [A]
* a = 0.5*mdu2; al=0.9 ; ala=al*a ;
**                               b=ala*al*((1.-al)**(-1.));
**            b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
* b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
* ik=masq Ki 'INFERIEUR' ala ;
**Ki=(ik*Ki)+((1.-ik)*a*(Ki+b)*(inve(a+Ki+b)));
**Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki+b2)*(inve(a+(Ki*Ki+b2)))));
* Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki*Ki+b3)*(inve(a+(Ki*Ki*Ki+b3)))));
 Ei = cd*(Ki**1.5)*(Echl**(-1.));

 Nuti=kcht $mod scal sommet (cnu*(Ki**0.5)*Echl);
 Finsi;
*............ Fin modèle K-L ...........................................

*................ modèle KLbr ..........................................
 Si KLbr;
 T2=T**(-2.);
 Ki  =kcht $mod scal sommet ((leps**2.)*(T2));

Si(NON Kbrey);
* Réalisabilité sur Ki dés/activée [A]
 a = 0.5*mdu2; al=0.9 ; ala=al*a ;
*                               b=ala*al*((1.-al)**(-1.));
*            b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
 b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
 ik=masq Ki 'INFERIEUR' ala ;
*Ki=(ik*Ki)+((1.-ik)*a*(Ki+b)*(inve(a+Ki+b)));
*Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki+b2)*(inve(a+(Ki*Ki+b2)))));
 Ki=(ik*Ki)+((1.-ik)*a*((Ki*Ki*Ki+b3)*(inve(a+(Ki*Ki*Ki+b3)))));
Finsi;

 Ei = (Ki**1.5)*(inve leps);
 Nuti=kcht $mod scal sommet (cnu*(Ki**0.5)*lmu);

 Finsi;
*............ Fin modèle KLbr ..........................................

*********************** Fin Equation sur Nuti et Teta(t) ************


********************** On conserve les résultats intermédiaires *****
 rv.inco.'TKTI'=kcht $mod scal sommet T ;
 rv.inco.'NUTI' = Nuti ;

****** Etape de diffusion equation sur K ***************************
*DFDT
* mess ' Etape de convection/diffusion + DFDT sur K' ;
 rxk.'KOPT'.'KFORM'=0 ;
 rxk.'KOPT'.'IDCEN'=1 ;
 rxk.'LISTINCO'=MOTS MI1 ;
 rxk.'IARG'=3;
 rxk.'ARG1'=Ro;
 rxk.'ARG2'=kcht $mod scal sommet Ki;
 rxk.'ARG3'=rx.'ARG4';
 sk1 mak1= DFDT rxk ;

 rxk.'KOPT'.'KFORM'=KFORM;
*TSCAL
 rxk.'IARG' = 1   ;
 rxk.'ARG1' = kcht $mod scal sommet (Mus2 + (Muto * (1./sgk)));
 sko mak = LAPN rxk ;

*MDIA
*mess ' MDIA pour les termes Bas-Reynolds' ;
 rxk.'NOMOPER'=mot 'MDIA_K';
 rxk.'IARG'=1;
 rxk.'ARG1'=kcht $mod scal sommet Epsok ;
 skm matk = MDIA rxk ;
*St2 = kcht $mod scal sommet comp MI1
*                    (nomc MI1  (Dg*Sfp));

 sk  = sk1 et sko ;
 mak = mak1 et mak et matk;

 climk= nomc MI1 (rv.'climk') ;

 Si PERI ;
  MPI1=MI1;
  mak=mak et (kops 'CHANINCO' mati
       ('MOTS' 'LX' 'UX')
       ('MOTS' (chai 'LX' MPI1) MPI1)
       ('MOTS' 'FLX' 'FX')
       ('MOTS' (chai 'LX' MPI1) MPI1) );
  sk =sk et (exco sti 'FLX' (chai 'LX' MPI1));
 Finsi;

 mtinv.'MATASS'=mak ;
 mtinv.'MAPREC'=mak ;
 mtinv.'XINIT'=nomc MI1 Ki ;
 mtinv.'IMPINV'=0;
 mtinv.'TYPINV'=TYPINV;

* mess 'Equation de diffusion sur K ';
 K1='KRES' mak 'TYPI' mtinv 'CLIM' climk 'SMBR' sk 'IMPR' 0 ;
 K1=(abs K1) + 1.e-10 ;
 K1=kcht $mod scal sommet (nomc 'SCAL' K1);
*mess 'maxi mini K1 apres Diff ' (maxi K1) (mini K1) ;

****** Etape de diffusion equation sur Epsilon *********************
Si ((NON Kkl) et (NON KLbr));
*DFDT
*mess ' Etape de diffusion DFDT sur EPSILON' ;
 rxe.'KOPT'.'KFORM'=0 ;
 rxe.'KOPT'.'IDCEN'=1 ;
 rxe.'LISTINCO'=MOTS MI2 ;
 rxe.'IARG'=3;
 rxe.'ARG1'=Ro;
 rxe.'ARG2'=kcht $mod scal sommet Ei;
 rxe.'ARG3'=rx.'ARG4';
 se1 mae1= DFDT rxe ;

 rxe.'KOPT'.'KFORM'=KFORM;
*TSCA
 rxe.'IARG'=1;
 rxe.'ARG1' = kcht $mod scal sommet (Mus2 + (Muto * (1./sge)));
 seo mae = LAPN rxe ;

*MDIA
*mess ' MDIA pour les termes Bas-Reynolds' ;
 rxe.'NOMOPER'=mot 'MDIA_E';
 rxe.'IARG'=1;
 rxe.'ARG1'=kcht $mod scal sommet Epsoe ;
 ske maed = MDIA rxe ;
*St2 = kcht $mod scal sommet comp MI2
*                    (nomc MI2  (Dg*Sfp));

 se  = se1 et seo ;
 mae = mae1 et mae et maed;

 clime= rv.'clime' ;

 Si PERI ;
  MPI1=MI2;
  mae=mae et (kops 'CHANINCO' mati
       ('MOTS' 'LX' 'UX')
       ('MOTS' (chai 'LX' MPI1) MPI1)
       ('MOTS' 'FLX' 'FX')
       ('MOTS' (chai 'LX' MPI1) MPI1) );
  se =se et (exco sti 'FLX' (chai 'LX' MPI1));
 Finsi;

 mtinv.'MATASS'=mae ;
 mtinv.'MAPREC'=mae ;
 mtinv.'XINIT'=nomc MI2 Eo ;
 mtinv.'IMPINV'=0;
 mtinv.'TYPINV'=TYPINV;

* mess 'Equation de diffusion sur Epsilon ';
 E1='KRES' mae 'TYPI' mtinv 'CLIM' clime 'SMBR' se 'IMPR' 0 ;
 E1=(Abs E1) + 1.e-10 ;
 E1=kcht $mod scal sommet (nomc 'SCAL' E1);
* mess 'maxi mini E1 apres Diff ' (maxi E1) (mini E1) ;
Finsi;

****** Cas K-L *****************************************************
Si Kkl ;
* mess 'Modèle K-L ';
 E1=cd*(K1**1.5)*(rv.inco.'Echl'**(-1.));
 E1=kcht $mod scal sommet (nomc 'SCAL' E1);
* mess 'maxi mini E1 apres Diff ' (maxi E1) (mini E1) ;
Finsi;

****** Avancement en temps ******************************

Si(NON Kbrey);
 Mut = Rocnu * K1 * K1 * (INVE E1) ;

* Réalisabilité sur K1 dés/activée [A]
 a = 0.5*mdu2; al=0.9 ; ala=al*a ;
*                               b=ala*al*((1.-al)**(-1.));
*            b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
 b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
 ik=masq K1 'INFERIEUR' ala ;
*K1=(ik*K1)+((1.-ik)*a*(K1+b)*(inve(a+K1+b)));
*K1=(ik*K1)+((1.-ik)*a*((K1*K1+b2)*(inve(a+(K1*K1+b2)))));
 K1=(ik*K1)+((1.-ik)*a*((K1*K1*K1+b3)*(inve(a+(K1*K1*K1+b3)))));

* Réalisabilité sur E1 dés/activée [A]  (nut P < 10 epsilon (Menter))
* iEmin =10*((K1*((ABS P)**0.5)) inve) ;
 iEmin =(((cnu*K1*K1*(ABS P)) inve)*ksim)**0.5 ;
 a =iEmin    ; al=0.7 ; ala=al*a ;
*                               b=ala*al*((1.-al)**(-1.));
             b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
*b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
 iE1 = inve E1 ;
 ik=masq iE1 'INFERIEUR' ala ;
*iE1=(ik*iE1)+((1.-ik)*a*(iE1+b)*(inve(a+iE1+b)));
 iE1=(ik*iE1)+((1.-ik)*a*((iE1*iE1+b2)*(inve(a+(iE1*iE1+b2)))));
*iE1=(ik*iE1)+((1.-ik)*a*((iE1*iE1*iE1+b3)*(inve(a+(iE1*iE1*iE1+b3)))));
 E1  = inve iE1 ;

finsi;


Si Kbw ;
*mess ' Réalisabilité Bw ';
* Réalisabilité sur Nut dés/activée [A]
Mutmax = Ro * Bw * K1 * (inve (P**0.5));
a =Mutmax    ; al=0.9 ; ala=al*a ;
*                               b=ala*al*((1.-al)**(-1.));
*            b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
 b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
ik=masq Mut 'INFERIEUR' ala ;
*Mut=(ik*Mut)+((1.-ik)*a*(Mut+b)*(inve(a+Mut+b)));
*Mut=(ik*Mut)+((1.-ik)*a*((Mut*Mut+b2)*(inve(a+(Mut*Nut+b2)))));
 Mut=(ik*Mut)+((1.-ik)*a*((Mut*Mut*Mut+b3)*(inve(a+(Mut*Mut*Mut+b3)))));
*Teta1 = K1 * (inve E1);
mess ' MUT1?';
E1 = Rocnu * K1 * K1 * (inve Mut) ;
mess ' MUT1?  OK';
Finsi  ;

  Si Filtre;
  Ech = (K1**1.5)*(inve E1);
  Echmax=rv.inco.'Echl';
 mess ' ** Modèle K-Epsilon filtré **';
 a =Echmax    ; al=0.9; ala=al*a ;
*                               b=ala*al*((1.-al)**(-1.));
*            b2=ala*(1.+ (ala*al) - ala)*((1.-al)**(-1.));
 b3 = ala*(1.+ (ala*ala*al) - (ala*ala))*((1.-al)**(-1.));
 ik=masq Ech 'INFERIEUR' ala ;
*Ech=(ik*Ech)+((1.-ik)*a*(Ech+b)*(inve(a+Ech+b)));
*Ech=(ik*Ech)+((1.-ik)*a*((Ech*Ech+b2)*(inve(a+(Ech*Ech+b2)))));
 Ech=(ik*Ech)+((1.-ik)*a*((Ech*Ech*Ech+b3)*(inve(a+(Ech*Ech*Ech+b3)))));
 E1 = (K1**1.5)*(inve Ech)       ;
 Mut = Rocnu * K1 * K1 * (inve E1) ;
  Finsi ;

Si (NON KLbr);
  Mut = Rocnu * K1 * K1 * (INVE E1) ;
Finsi;

****** Cas KLbr ****************************************************
Si KLbr ;
 E1=(K1**1.5)*(leps**(-1.));
 E1=kcht $mod scal sommet E1;
 Mut = Rocnu * (K1**0.5) * lmu;
Finsi;
****** Cas KLbr Fin ************************************************


  dk= maxi (abs ( K1 - Ko))/(maxi Ko) ;
  Ko=K1 ;
  de= maxi (abs ( E1 - Eo))/(maxi Eo) ;
  Eo=E1 ;
  Ksi = Mut * P * (inve E1) * iRo ;
  rv.inco.'Ksi'=Ksi ;
  rv.inco.MI1=kcht $mod scal sommet K1;
  rv.inco.MI2=kcht $mod scal sommet E1;


  rv.inco.'MUF' = kcht $mod scal sommet (Mut + Mu);
  rv.inco.'Muto'= kcht $mod scal sommet Mut;
  rv.inco.MI2=kcht $mod scal sommet E1;


   FIDT=RV.'FIDT';
   INDF = IPT - (IPT/FIDT * FIDT);
Si (EGA INDF 0);
Mess ' Ecart relatif entre deux pas de temps : sur K ' dk
 ' sur Epsilon ' de ;
Finsi ;
****** Historiques **************************************
Si (exist rv 'HIST') ;
   H=rv.'HIST';
   Si(non(exist H 'KFIH'));
   H.'KFIH'=20;
   Finsi ;
   KFIH=H.'KFIH';
   INDH = IPT - (IPT/KFIH * KFIH);
  Si (EGA INDH 0);

    Si (exist (rv.'HIST') MI1);

  M1=H.(CHAI '$' MI1);
  EV1=H.MI1;
  ltps=extr EV1 'ABSC' 1 ;
  nlt=dime ltps;
  Si(ega nlt 0);
  tps= 0. ;

    Sinon;

  tps=(extr ltps nlt) + dt ;
  finsi ;
  ltps=ltps et (prog tps) ;
  nb=nbel M1 ;
  repeter bloc1 nb;
  pi=poin M1 &bloc1;
  lki =(extr EV1 'ORDO' &bloc1) et (prog (extr Ko 'SCAL' pi));
  Si (EGA &bloc1 1);
  evi = (evol 'MANU' 'TEMPS' ltps MI1 lki) ;
  Sinon;
  evi = evi et (evol 'MANU' 'TEMPS' ltps MI1 lki) ;
  Finsi ;
  Fin bloc1 ;
  H.MI1=evi;

    Finsi ;

    Si (exist (rv.'HIST') MI2);

  M2=H.(CHAI '$' MI2);
  EV2=H.MI2;
  nb=nbel M2 ;
  repeter bloc1 nb;
  pi=poin M2 &bloc1;
  lki =(extr EV2 'ORDO' &bloc1) et (prog (extr Eo 'SCAL' pi));
  Si (EGA &bloc1 1);
  evi = (evol 'MANU' 'TEMPS' ltps MI2 lki) ;
  Sinon;
  evi = evi et (evol 'MANU' 'TEMPS' ltps MI2 lki) ;
  Finsi ;
  Fin bloc1 ;
  H.MI2=evi;
   Finsi ;
  Finsi ;
Finsi ;
****** Petit menage avant retour ************************
   rx.'KOPT'.'IDIV'=IDIV ;

 as2 ama1  = 'KOPS' 'MATRIK' ;

*Mess ' FIN PROC KEPSILON' ;
RESPRO as2 ama1 ;

FINPROC ;