Télécharger unpas.procedur

Retour à la liste

* UNPAS     PROCEDUR  JK148537  25/12/17    21:15:04     12424          
'DEBPROC' UNPAS PRECED*'TABLE';
*----------------------------------------------------------------------*
*                          PROCEDURE UNPAS                             *
*                                                                      *
* Calcul d'un increment de solution en grand deplacement plastique par *
* la methode des residus.                                              *
*                                                                      *
* Les differentes configurations qui interviennent sont :              *
* WTAB.'FOR0'    : configuration de debut de calcul                    *
* WTAB.'GE0_DEB' :                                                     *
*      -> en GRANDS_DEPLACEMENTS ou FEFP_FORMULATION (1) :             *
*         configuration de debut de pas                                *
*      -> dans les autres cas (2) :                                    *
*         configuration de debut de calcul (= WTAB.'FOR0)              *
*  Dans le cas (1), la configuration WTAB.'GE0_DEB' est actualisee     * 
*  dans PASAPAS avec GEOM2.                                            *
*  Dans le cas (2), il n'y a pas d'actualisation de la configuration   *
*  au cours du calcul.                                                 * 
*                                                                      *
*  Pour les calculs DYNAMIQUE, un schema de Newmark (implicite)        *
*  avec gamma = 1/2 et beta = 1/4 est utilise.                         *
*                                                                      *
*----------------------------------------------------------------------*
LOG_CNTRL = FAUX;
*CB215821 : Recuperation de XPETIT (07/12/2016)
XPETIT   = 'VALE' 'PETI'    ;
XGRAND   = 'VALE' 'GRAND'   ;
*
WTAB  = PRECED.'WTABLE'      ;
conti = PRECED.'CONTINUATION';
estim = PRECED.'ESTIMATION'  ;
*
* Liste de composantes utiles
MXMYMZ = 'MOTS' 'MX'  'MY'  'MZ'  'MT'  'FP'  'FPQ'  'FTP'
                'IMX' 'IMY' 'IMZ' 'IMT' 'IFP' 'IFPQ' 'IFTP' ;
MXMFLX = 'MOTS' 'MX'  'MY'  'MZ'  'MT'  'FLX' 'FP'  'FPQ'  'FTP'
                'IMX' 'IMY' 'IMZ' 'IMT' 'FLX' 'IFP' 'IFPQ' 'IFTP' ;
MNPRIM = 'MOTS' 'UX'  'UY'  'UZ'  'UR'  'UT'   'RX'   'RY'  'RZ'  'RT'
                      'P'   'PQ'  'TP'   'ALFA' 'BETA'
                'IUX' 'IUY' 'IUZ' 'IUR'  'IUT'  'IRX'  'IRY' 'IRZ' 'IRT'
                'IP' 'IPQ'  'ITP' 'IALF' 'IBET' ;
MLPRIM = MNPRIM 'ET' ('MOTS' 'LX');
MNDUAL = 'MOTS' 'FX'  'FY'  'FZ'  'FR'   'FT'   'MX'   'MY'  'MZ'  'MT'
                      'FP'  'FPQ' 'FTP'  'FALF' 'FBET'
                'IFX' 'IFY' 'IFZ' 'IFR'  'IFT'  'IMX'  'IMY' 'IMZ' 'IMT'
                'IFP' 'IFPQ' 'IFTP' 'IFAL' 'IFBE' ;
MLDUAL = MNDUAL 'ET' ('MOTS' 'FLX');
MLDEPL = 'MOTS' 'UX' 'UY' 'UZ' 'UR' 'UT'  'IUX' 'IUY' 'IUZ' 'IUR' 'IUT'
                'ALFA' 'BETA' 'IALF' 'IBET' ;
MLROTA = 'MOTS' 'RX' 'RY' 'RZ' 'RT' 'P' 'PQ' 'TP'
                'IRX' 'IRY' 'IRZ' 'IRT' 'IP' 'IPQ' 'ITP';
MLDEFOR = 'MOTS'  'EPXX' 'EPYY' 'EPZZ' 'EPSS' 'EPTT' 'EPRR'
                  'GAXY' 'GAXZ' 'GAYZ' 'GAST' 'GASN' 'GATN'
                  'GARZ' 'GART' 'GAZT' 'GXY '
*                 'CX  ' 'CY  ' 'CZ  '
                  'EPSE' 'EPS' ;
MLPRINC = 'MOTS' 'SI11' 'SI22' 'SI33';
* Pour MOD_LIAISON :
MLPRIM_LIA = MLPRIM ;
MVPRIM_LIA = 'MOTS' 'VTX' 'VTY' 'VTZ' 'VTR' 'VTT' 'VWX' 'VWY' 'VWZ' 'VWT'
                    'VVP' 'VVPQ' 'VVTP' 'VALF' 'VBET'
                    'IVTX' 'IVTY' 'IVTZ' 'IVTR' 'IVTT' 'IVWX' 'IVWY' 'IVWZ' 'IVWT'
                    'IVVP' 'IVPQ' 'IVTP' 'IVAL' 'IVBE' 'VLX';
MLDUAL_LIA = MLDUAL ;
 
*=DEB===== Formulation HHO = Definition des DDLs ===============================
* Definition des DDLs primaux : inconnues de deplacements des cellules et des Faces
MNPRIM_HHO = 'MOTS'
*                'UXC0' 'UXC1' 'UXC2' 'UXC3' 'UXC4' 'UXC5' 'UXC6' 'UXC7' 'UXC8' 'UXC9'
*                'UYC0' 'UYC1' 'UYC2' 'UYC3' 'UYC4' 'UYC5' 'UYC6' 'UYC7' 'UYC8' 'UYC9'
*                'UZC0' 'UZC1' 'UZC2' 'UZC3' 'UZC4' 'UZC5' 'UZC6' 'UZC7' 'UZC8' 'UZC9'
                'UXF0' 'UXF1' 'UXF2' 'UXF3' 'UXF4' 'UXF5' 'UXF6' 'UXF7' 'UXF8' 'UXF9'
                'UYF0' 'UYF1' 'UYF2' 'UYF3' 'UYF4' 'UYF5' 'UYF6' 'UYF7' 'UYF8' 'UYF9'
                'UZF0' 'UZF1' 'UZF2' 'UZF3' 'UZF4' 'UZF5' 'UZF6' 'UZF7' 'UZF8' 'UZF9' ;
MNDUAL_HHO = 'MOTS'
*                'FXC0' 'FXC1' 'FXC2' 'FXC3' 'FXC4' 'FXC5' 'FXC6' 'FXC7' 'FXC8' 'FXC9'
*                'FYC0' 'FYC1' 'FYC2' 'FYC3' 'FYC4' 'FYC5' 'FYC6' 'FYC7' 'FYC8' 'FYC9'
*                'FZC0' 'FZC1' 'FZC2' 'FZC3' 'FZC4' 'FZC5' 'FZC6' 'FZC7' 'FZC8' 'FZC9'
                'FXF0' 'FXF1' 'FXF2' 'FXF3' 'FXF4' 'FXF5' 'FXF6' 'FXF7' 'FXF8' 'FXF9'
                'FYF0' 'FYF1' 'FYF2' 'FYF3' 'FYF4' 'FYF5' 'FYF6' 'FYF7' 'FYF8' 'FYF9'
                'FZF0' 'FZF1' 'FZF2' 'FZF3' 'FZF4' 'FZF5' 'FZF6' 'FZF7' 'FZF8' 'FZF9' ;
* Pour l'istant on ne considere que les inconnues (DEPL et FORC) de face pour
* le calcul de l'equilibre et l'acceleration de convergence
* D'ou les lignes commentees ci-dessus dans la definition de MNPRIM_HHO et MNDUAL_HHO
MNPRIM = MNPRIM 'ET' MNPRIM_HHO ;
MLPRIM = MLPRIM 'ET' MNPRIM_HHO ;
MNDUAL = MNDUAL 'ET' MNDUAL_HHO ;
MLDUAL = MLDUAL 'ET' MNDUAL_HHO ;
MLDEPL = MLDEPL 'ET' MNPRIM_HHO ;
*A Ameliorer : ne faire que si presence HHO et n'ajouter que les DDLs utilises par HHO !
*=FIN===== Formulation HHO =====================================================
*
* Pour stocker des informations necessaires aux calculs d'usure
WTAB.'POST_COFR' = 'TABL' ;
*----------------------------------------------------------------------*
*             Options de pasapas - parametres du calcul                *
*----------------------------------------------------------------------*
* 1- rigidite
IKSIA  = WTAB.'K_SIGMA';
IKTAN  = WTAB.'K_TANGENT' ;
ZKTASYM = 'TEXTE' '    ' ;
'SI' WTAB.'K_TANGENT_SYME' ;
  ZKTASYM = 'MOT' 'SYME' ;
'FINS';
IMPLP   = WTAB.'LIAISON_PERSISTANTE';
AUTAUG  = WTAB.'AUTOAUGM';
IRAUGLU = WTAB.'RAIDAUGM';
IRAUG   = FAUX;
'SI' IRAUGLU ;
  RIG_AUG = WTAB.'RIGIDITE_AUGMENTEE';
  'SI' ('NON' AUTAUG) ; 
    IRAUG = VRAI ;
  'FINS';
'FINS';
IRCON = WTAB.'RAIDCONST';
'SI' IRCON;
  RIG_CONS = WTAB.'RIGIDITE_CONSTANTE';
  MAI_CONS ='EXTR' RIG_CONS 'MAIL'    ;
'FINS';
*  declenchement du recalcul de la matrice
DLTAIT = WTAB.'DELTAITER' ;
ITRCLC =  1  '*' DLTAIT ;
'SI' (ITRCLC <  20) ; ITRCLC =  20; 'FINSI';
*
* 2- Type de formulation
IFEFP  = WTAB.'FEFP_FORMULATION' ;
LAG_TOT= WTAB.'LAG_TOT';
ISSTE  = WTAB.'SUBSTEPPING';
IFEFPUL= WTAB.'UPDATE_LAGRANGIAN';
LNLOC  = WTAB.'NLOC';
'SI' (LNLOC 'ET' ('EGA' WTAB.'NON_LOCAL' 'HELM'));
  TAHELM = WTAB.'HELMHOLTZ' ;
  NHELM =  TAHELM . 'N_VARI_NL' ;
'FINSI' ;
*
* 3- Type de materiau
ICERAM = WTAB.'CERAMIQUE' ;
IENDOM = WTAB.'ENDOMMAGEMENT';
IPLAVI = WTAB.'IPLAVI';
IVIEXT = WTAB.'VISCO_EXTERNE';
IVIDOM = WTAB.'VISCODOMMAGE';
IVISCO = WTAB.'VISCOPLASTIQUE';
*
* 4- Critere de convergence/non-convergence
ZMAXIT = WTAB.'MAXITERATION' ;
NITMA  = WTAB.'NITERINTER_MAX';
ZPREK  = WTAB.'PRECISINTER' ;
ZPREC  = WTAB.'PRECISION' ;
ZPRECM = WTAB.'PRECFLEX' ;
ZPRECHPP = ZPREC ;
zpreceps = WTAB.'PRECDEFO';
IFTOL  = 'NEG' WTAB.'FTOL' 'INCONNU' ;
'SI' IFTOL ;
  ZFTOL = 'ABS' WTAB.'FTOL' ;
'FINS';
IMTOL  = 'NEG' WTAB.'MTOL' 'INCONNU';
'SI' IMTOL ;
  ZMTOL = 'ABS' WTAB.'MTOL' ;
'FINS';
'SI' ('OU' IVISCO IVIDOM IVIEXT); ZPREK = 5.E-7 ; 'FINS';
'SI' IENDOM;                      ZPREK = ZPREC ; 'FINS';
NSOINCRN = WTAB.'SOUS_INCREMENT' ;
*
* 5- Type de calcul
IPILOT  = WTAB.'AUTOMATIQUE';
ISOL    = WTAB.'CONSOLIDATION';
IDYN    = WTAB.'DYNAMIQUE';
IGRD    = WTAB.'GRANDS_DEPLACEMENTS';
HYPDEF  = 'MOT' WTAB.'HYPOTHESE_DEFORMATIONS';
IJAUMA  = 'EGA' HYPDEF 'JAUMANN' ;
'OPTIO' 'EPSI' HYPDEF;
HPP_EPS = 'EGA' WTAB.'PREDICTEUR' 'HPP';
*
* 6- Chargement particulier
LOGDEF = WTAB.'CHAR_DEFI';
LOGPRE = WTAB.'CHAR_PRES' ;
LOGPIL = WTAB.'CHAR_PILO' ;
ITHER  = WTAB.'ITHER' ;
*
* 7- Acceleration de convergence
ZNACCE = VRAI;
*
* 8- Instant de debut et de fin de pas, pas de temps
TEMPS0 = conti.'TEMPS';
TI     = estim.'TEMPS';
PASDT  = WTAB.'DT' ;
TDEBUT = TEMPS0 ;
DTTAN  = 0. ;
*
* 9- Blocages mecaniques (BLOM)
'SI' WTAB.'CHAR_BLOM' ;
  BLOM1 = 'TIRE' PRECED.'BLOCAGES_MECANIQUES' 'BLOM' TI ;
  WTAB.'BLOCAGES_MECANIQUES' = BLOM1 ;
'FINSI' ;
ZCLIM0 = WTAB.'BLOCAGES_MECANIQUES';
ZCLIM  = ZCLIM0 ;
*----------------------------------------------------------------------*
*                Initialisation (Configuration GEOREF0)                *
*----------------------------------------------------------------------*
'SI' (IKTAN 'ET' ('NON' IPLAVI)) ;
  'MESS' 'IPLAVI faux : pas de matrice tangente ->'
         ' on utilise la rigidite elastique' ;
  IKTAN = FAUX ;
'FINS' ;
IPERT = WTAB.'K_TANGENT_PERT' 'ET' ('NON' LNLOC) 'ET' IPLAVI ;
* Matrice tangente : pas d'acceleration en cas de modele FEFP ou SSTE
'SI' IKTAN ;
  'SI' (IFEFP 'OU' ISSTE) ; ZNACCE = FAUX; 'FINS';
'FINS' ;
*
* Configurations de reference et de debut de pas
GEOREF0 = WTAB.'FOR0' ;
GEOM1   = WTAB.'GE0_DEB' ;
'FORM' GEOREF0 ;
*----------------------------------------------------------------------*
*           Modele complet                                             *
*                                                                      *
* Notations utilisees :                                                *
*   ZMODLI : modeles mecanique + poreux (non //)                       *
*   -> aucune parallelisation       alors ZMODL = ZMODLI               *
*   -> parallelisation comportement alors ZMODL = ZMODLI               *
*                                         MODRELOC est //              *
*   -> parallelisation automatique  alors ZMODL est //                 *
*----------------------------------------------------------------------*
* PAS_MODL : mise a jour des indices de PRECED.WTABLE relatifs aux 
*            modeles si PRECED.WTABLE.MODELE a ete modifie.
PAS_MODL PRECED ;
ZMODLI = WTAB.'MO_TOT' ;
ZMODL  = ZMODLI ;
'SI' WTAB.'FOR_MECA' ;
  PARALLEL = FAUX ; PARTLOCA = FAUX ;
  'SI' ('EGA' WTAB.'PROCESSEURS' 'COMPORTEMENT') ;
    PARTLOCA = VRAI ;
    MODRELOC = 'PART' 'ARLE' ZMODLI WTAB.'NBPART' ;
  'FINS' ;
  'SI' ('EGA' WTAB.'PROCESSEURS' 'AUTOMATIQUE') ;
    PARALLEL = VRAI ;
    ZMODL = 'PART' 'ARLE' ZMODLI WTAB.'NBPART' ;
    'OPTI' 'PARA' VRAI ;
  'FINS' ;
*
* Si le modele a change, adapter egalement les MCHAMLS
  'SI' ('NEG' ZMODLI WTAB.'MO_TOT_PREC') ;
*
    CONTI.'CONTRAINTES' = ('REDU' CONTI.'CONTRAINTES'  ZMODLI) '+' ('ZERO' ZMODLI 'CONTRAINTES') ;
    CONTI.'DEFORMATIONS'= ('REDU' CONTI.'DEFORMATIONS' ZMODLI) '+' ('ZERO' ZMODLI 'DEFORMAT'   ) ;
    'SI' IPLAVI ;
      CONTI.'DEFORMATIONS_INELASTIQUES'= ('REDU' CONTI.'DEFORMATIONS_INELASTIQUES' ZMODLI) '+' ('ZERO' ZMODLI 'DEFINELA') ;
      CONTI.'VARIABLES_INTERNES'       = ('REDU' CONTI.'VARIABLES_INTERNES'        ZMODLI) '+' ('ZERO' ZMODLI 'VARINTER') ;
    'FINS';
*   
    ETAT1       = PAS_ETAT PRECED TEMPS0 ;
    WTAB.'MAT1' = PAS_MATE PRECED ETAT1  ;
*  
    'SI' ITHER ;
      WTAB.'ETHER1'    = ('REDU' WTAB.'ETHER1'   ZMODLI) '+' ('ZERO' ZMODLI 'DEFORMAT') ;
      'SI' WTAB.'POR1';
        WTAB.'MSRTHER1'= ('REDU' WTAB.'MSRTHER1' ZMODLI) '+' ('ZERO' ZMODLI 'DEFORMAT') ;
      'FINS';
    'FINS';
*  
    'SI' LOGDEF ; 
      WTAB.'DEFOR1' = ('REDU' WTAB.'DEFOR1' ZMODLI) '+' ('ZERO' ZMODLI 'DEFORMAT') ;
    'FINS';
*
  'FINS';
'FINSI';
*----------------------------------------------------------------------*
*                 Champs materiau et caracteristiques                  *
*                                                                      *
* Notations :                                                          *
*   ZMAT1  : champs materiau et caracteristiques a TEMPS0              *
*   ZMAT20 : champs materiau et caracteristiques a TI dans GEOREF0     *
*   ZMAT21 : champs materiau et caracteristiques a TI dans GEOM1       *
*   ZMAT22 : champs materiau et caracteristiques a TI dans GEOM2       *
*----------------------------------------------------------------------*
ZMAT1 = WTAB.'MAT1' 'REDU' ZMODL;
*
ETAT2  = PAS_ETAT PRECED TI;
ZMAT20 = PAS_MATE PRECED ETAT2;
*
'SI' WTAB.'FOR_MECA' ;
  ZMAT21 = 'REDU' ZMAT20 ZMODL ;
  'SI' IGRD;
    'FORM' GEOM1 ;
    ZMAT21 = 'CONF' ZMODL ZMAT21 ;
    'FORM' GEOREF0 ;
  'FINSI' ;
  'SI' WTAB.'POR1';
    MA_POR = 'REDU' ZMAT20 WTAB.'MOD_POR' ;
  'FINS';
'FINSI' ;
*----------------------------------------------------------------------*
*           Chargement a imposer a l'instant TI      (F^ext_n+1)       *
*----------------------------------------------------------------------*
ZFEXT2 = 'VIDE' 'CHPOINT' / 'DISCRET';
*
'SI' ('EXIS' WTAB.'CHARGEMENT' 'MECA');
  ZFEXT2 = 'TIRE' WTAB.'CHARGEMENT' 'MECA' TI;
  'SI' ('NEG' ('TYPE' ZFEXT2) 'CHPOINT ');
    'MESS' '*** Erreur dans la definition du chargement ( MECA ) ***';
    'ERRE' '*** Le type du champ n est pas CHPOINT ***';
  'FINS';
'FINS';
*
'SI' ('EXIS' WTAB.'CHARGEMENT' 'FORC');
  F2_FOR = 'TIRE' WTAB.'CHARGEMENT' 'FORC' TI;
  'SI' ('NEG' ('TYPE' F2_FOR) 'CHPOINT ');
    'MESS' '*** Erreur dans la definition du chargement ( FORC ) ***';
    'ERRE' '*** Le type du champ n est pas CHPOINT ***';
  'FINS';
  ZFEXT2 = ZFEXT2 '+' F2_FOR ;
'FINS';
*
'SI' (LOGPRE 'ET' ('NON' IGRD)) ;
  ZPEXT = 'TIRE' WTAB.'CHARGEMENT' 'PRES' TI ;
  'SI' ('EXIS' WTAB 'MAT_PRE') ;
    ZFPEXT = 'BSIG' WTAB.'MOD_PRE' ZPEXT WTAB.'MAT_PRE' ;
  'SINON' ;
    ZFPEXT = 'BSIG' WTAB.'MOD_PRE' ZPEXT ;
  'FINS' ;
  ZFEXT2 = ZFEXT2 '+' ZFPEXT ;
'FINS' ;
*
'SI' ('EXIS' WTAB.'CHARGEMENT' 'DIMP');
  F2_DEP = 'TIRE' WTAB.'CHARGEMENT' 'DIMP' TI;
  ZFEXT2 = ZFEXT2 + F2_DEP;
'FINS';
*
'SI' ('EXIS' WTAB.'CHARGEMENT' 'DINC');
  F2_INC = 'TIRE' WTAB.'CHARGEMENT' 'DINC' TI;
  F2_mail = extr f2_inc maillage;
  zclim_mail = zclim extrai 'MAILLAGE';
  F2_Mcmp = zclim_mail elem appu larg f2_mail noverif ;
  zclim_inc = zclim redu f2_mcmp;
  f2_base = (zclim_inc * conti.'DEPLACEMENTS') exco 'FLX' 'NOID' 'FLX';
  ZFEXT2 = ZFEXT2 + F2_INC + F2_base ;
'FINS';
*
'SI' LOGDEF ;
  WTAB.'DEFOR2' = 'TIRE' WTAB.'CHARGEMENT' 'DEFI' TI;
'FINS';
*
'SI' WTAB.'FOR_MECA' ;
  'SI' ITHER ;
*
    'SI' WTAB.'CHAR_THE';
      WTAB.'TET2' = 'TIRE' WTAB.'CHARGEMENT' 'T   ' TI;
    'FINSI';
    'SI' WTAB.'FOR_THER';
      WTAB.'TET2' = estim.'TEMPERATURES' ;
    'FINSI';
*
*   Champs de temperature en debut (ZTEMP1) et fin de pas (ZTEMP2)
*   Si les maillages mecanique/thermique sont differents,
*   on projette le champ de temperature sur le modele mecanique
    'SI' WTAB.'PROJECTION';
      ZMOMAI = 'EXTR' ZMODL 'MAILLAGE' ;
      'SI' ('EGA' ('TYPE' WTAB.'TET1') 'CHPOINT') ;
        ZTEMP1 = 'PROI' ZMOMAI ('CHAN' 'CHAM' WTAB.'TET1' WTAB.'MOD_THE' 'NOEUD') ;
        ZTEMP2 = 'PROI' ZMOMAI ('CHAN' 'CHAM' WTAB.'TET2' WTAB.'MOD_THE' 'NOEUD') ;
      'SINON' ;
        ZTEMP1 = 'PROI' ZMOMAI WTAB.'TET1' ;
        ZTEMP2 = 'PROI' ZMOMAI WTAB.'TET2' ;
      'FINSI' ;
*   Sinon, on vient seulement recuperer les champs de temperature
    'SINON' ;
      ZTEMP1 = WTAB.'TET1';
      ZTEMP2 = WTAB.'TET2';
    'FINSI' ;
*
    'SI' ('EGA' ('TYPE' ZTEMP2) 'CHPOINT');
      MCHTETA2 = 'CHAN' 'CHAM' ZTEMP2 ZMODLI 'STRESSES' 'TEMPERATURES' ;
    'SINON';
      MCHTETA2 = ZTEMP2 ;
    'FINS';
*
*   REDU au cas ou le modele soit //
    ETT0 = 'REDU' WTAB.'ETHER1' ZMODL;
    'SI' WTAB.'POR1';
      MSRTT0 = 'REDU' WTAB.'MSRTHER1' ZMODL;
    'FINS';
*
  'FINSI';
'FINSI';
*----------------------------------------------------------------------*
*         Dynamique : ajout au second membre de                        *
*                     F0 + 4/DT*M*V0 - B0*SIG0                         *
*----------------------------------------------------------------------*
* Matrice de masse
'SI' (IDYN 'OU' WTAB.'FREQUENTIEL') ;
  'SI' (('NON' ('EXIS' WTAB 'MASSE')) 'OU' (IGRD 'ET' IDYN)) ;
    'SI' IGRD; 'FORM' GEOM1 ; 'FINSI';
    WTAB.'MASSE' = 'MASS' ZMODLI ('CHAN' 'MASSE' ZMAT21 ZMODLI) ;
    'SI' IGRD; 'FORM' GEOREF0 ; 'FINSI';
    'SI' (IDYN 'ET' WTAB.'MASSCONST') ;
      WTAB.'MASSE' = WTAB.'MASSE' 'ET' WTAB.'MASSE_CONSTANTE';
    'FINS';
  'FINS';
'FINS';
*
'SI' IDYN ;
  UNSURH = 1.D0 '/' PASDT ;
  'SI' ('EGA' WTAB.'FREA1' 'INCONNU');
    F1 = 'VIDE' 'CHPOINT' / 'DISCRET';
*
*   Forces exterieures a TEMPS0
    'SI' ('EXIS' WTAB.'CHARGEMENT' 'MECA');
      F1M = 'TIRE' WTAB.'CHARGEMENT' 'MECA' TEMPS0 ;
      F1  = F1 '+' F1M ;
    'FINS';
*
    'SI' ('EXIS' WTAB.'CHARGEMENT' 'FORC');
      F1F = 'TIRE' WTAB.'CHARGEMENT' 'FORC' TEMPS0 ;
      F1  = F1 '+' F1F ;
    'FINS';
*
    'SI' WTAB.'PROCEDURE_CHARMECA';
      TFF1 = CHARMECA PRECED TEMPS0 ;
      'SI' ('EXIS' TFF1 'ADDI_SECOND');
        FF1 = TFF1.'ADDI_SECOND';
      'FINS';
      F1 = F1 '+' FF1;
    'FINS';
*
    'SI' LOGPRE ;
      ZPEXT0 = 'TIRE' WTAB.'CHARGEMENT' 'PRES' TEMPS0 ;
      'SI' IGRD; 'FORM' GEOM1 ; 'FINSI';
      'SI' ('EXIS' WTAB 'MAT_PRE') ;
        FF1 = 'BSIG' WTAB.'MOD_PRE' ZPEXT0 WTAB.'MAT_PRE' ;
      'SINON' ;
        FF1 = 'BSIG' WTAB.'MOD_PRE' ZPEXT0 ;
      'FINS' ;
      'SI' IGRD; 'FORM' GEOREF0 ; 'FINSI';
      F1 = F1 '+' FF1;
    'FINS';
*
*   Forces interieures a TEMPS0
*HHO : Modifications appel a BSIGMA (ajout du champ de deplacements necessaire en HHO)
    'SI' IGRD; 'FORM' GEOM1 ; 'FINSI';
    LAF0 = 'BSIG' ZMODLI conti.'CONTRAINTES' ZMAT21 conti.'DEPLACEMENTS' ;
    'SI' IGRD; 'FORM' GEOREF0 ; 'FINSI';
    'SI'  IRCON  ;
      LAF0 = LAF0 'ET' ('REDU' MAI_CONS (RIG_CONS '*' conti.'DEPLACEMENTS'));
    'FINS';
*
*   FREA1 : (masse*acceleration initiale)+(amortissement*vitesse initiale)
    WTAB.'FREA1' = F1 '-' LAF0 ;
  'FINS';
*
  'SI' IMPLP ;
*  forces d'acceleration au debut du pas
    'SI' ('NEG' WTAB.'AMORTISSEMENT' 'INCONNU');
      FF4 = WTAB.'AMORTISSEMENT' '*' conti.'VITESSES';
      WTAB.'FMAN' = WTAB.'FREA1' - FF4 ;
    'SINON' ;
      WTAB.'FMAN' = WTAB.'FREA1' ;
    'FINS';
  'FINS';
*
  FF = WTAB.'MASSE' '*' conti.'VITESSES';
  FF4 = 4. '/' PASDT '*' FF; 'DETR' FF;
*
* partie du second membre qui ne depend que des informations du pas prec
  WTAB.'FREA1' = FF4 '+' WTAB.'FREA1';
'FINS';
*
* Advection
'SI' WTAB.'NVSTNL'                                        ;
  LOGADV = 'EXISTE' WTAB.'CHARGEMENT' 'ADVE'              ;
  'SI' LOGADV                                             ;
    UADV = 'TIRE' WTAB.'CHARGEMENT' 'ADVE' TI             ;
  'FINSI'                                                 ; 
'FINSI'                                                   ;
*----------------------------------------------------------------------*
*     Consolidation : ajout au second membre de                        *
*                     -B0*SIG0 + DT*(1-TETA)*FI0 + DT*H*P              *
*----------------------------------------------------------------------*
'SI' ISOL ;
  FF  = 'BSIG' WTAB.'MOD_POR' conti.'CONTRAINTES' ;
  FF4 = 'EXCO' WTAB.'MOT_POR' FF
               WTAB.'MOT_POR' 'NOID' 'NATURE' 'DISCRET' ;
  ZFEXT2 = ZFEXT2 '+' FF4  ; 'DETR' FF ;
*
*           ---- traitement des flux si besoin ----
  'SI' ('EXIS' WTAB.'CHARGEMENT' 'FLUX' ) ;
    FLUXT0= 'TIRE' WTAB.'CHARGEMENT' 'FLUX' TEMPS0;
    FLUXTI= 'TIRE' WTAB.'CHARGEMENT' 'FLUX' TI  ;
    FACFLU =  -1. * PASDT ;
    FLUXT =   ( FACFLU * (1 - WTAB.'TETA') * FLUXT0 )
            + ( FACFLU * WTAB.'TETA' *FLUXTI ) ;
    ZFEXT2 = ZFEXT2 '+' FLUXT ;
    'DETR' FLUXT ; 'DETR' FLUXT0; 'DETR' FLUXTI;
  'FINS' ;
'FINS';
*----------------------------------------------------------------------*
*                  Pilotage indirect : preparation                     *
*----------------------------------------------------------------------*
'SI' LOGPIL ;
  ZFPILIN = WTAB.'FORCES_PILOTEES' 'ET' WTAB.'DEPLACEMENTS_PILOTES' ;
  DTAU = ('IPOL' WTAB.'PARAMETRE_DE_PILOTAGE' TI) -
         ('IPOL' WTAB.'PARAMETRE_DE_PILOTAGE' TEMPS0) ;
  ETA0 = 'EXTR' PRECED.'COEFFICIENT_DE_PILOTAGE' 
                             ('DIME' PRECED.'COEFFICIENT_DE_PILOTAGE') ;
  D_ETA=0.;                                                   
'FINS' ;
*----------------------------------------------------------------------*
*         Partie constante du second membre sur le pas de temps        *
*----------------------------------------------------------------------*
ZFCONST1 = ZFEXT2 ;
'SI' IDYN ;
  ZDYFEXT = 'ENLE' ZFEXT2 'FLX';
  ZFEXT2  = ZFEXT2 '+' WTAB.'FREA1';
'FINSI' ;
ZFLXB = 'EXCO' ZFEXT2 'FLX' 'NOID' 'FLX' ;
*----------------------------------------------------------------------*
*   Recuperation de valeurs de WTAB pour initialisations eventuelles   *
*----------------------------------------------------------------------*
ZRIBLO_M = 'MOT' 'INCONNU' ;
'SI' WTAB.'FOR_MECA' ;
*
  ZXDENO  = WTAB.'XDENO'  ;
  ZXDENOM = WTAB.'XDENOM' ;
* XDENOo  = WTAB.'XDENO';
* XDENOMo = WTAB.'XDENOM';
  'SI' LOGDEF ;
    ZDEFOR1 = WTAB.'DEFOR1' ;
    ZDEFOR2 = WTAB.'DEFOR2' ;
  'FINS';
  'SI' ITHER ;
    ZTET1 = WTAB.'TET1' ;
    ZTET2 = WTAB.'TET2' ;
  'FINS';
  ZLASTKTAN = WTAB.'LASTKTAN' ;
*
  'SI' ('NEG' WTAB.'RIBLO_M' 'INCONNU') ;
    ZRIBLO_M = WTAB.'RIBLO_M' ;
    ZLISEA_M = WTAB.'LISEA_M' ;
  'FINS' ;
*
  ZFNONL = WTAB.'FNONL' ;
  'SI' ('NEG' WTAB.'INCREMENT' 'INCONNU');
    INCRPREC = WTAB.'INCREMENT';
  'FINS';
*
* teste t'on les moments  ?
  TSTMOM = 'EGA' ('VALE' 'MODE') 'PLANGENE' ;
  TSTMOM = TSTMOM 'OU' ('EXIS' ZMAT21 'EPAI') 'OU' ('EXIS' ZMAT21 'INRY') 
                  'OU' ('EXIS' ZMAT21 'MODS') 'OU'
           (('EXIS' ZMODLI 'ELEM' 'JOI1') 'ET' ('EXIS' ZMAT21 'V1X ')
                                          'ET' ('EXIS' ZMAT21 'V1Y '));
* teste t'on les POREUX   ?
  'SI' WTAB.'POR1' ; TSTMOM=VRAI ; 'FINS';
*
  mocrit = 'CHAI' 'Critere' ;
  moflex = 'CHAI' 'Crit.flex' ;
  IKLFFF = FAUX ;
  'SI'TSTMOM;    
    'SI' IFTOL; 'SI' IMTOL;
       IKLFFF=VRAI;
    'FINS'; 'FINS';
  'SINON' ;
    'SI' IFTOL; IKLFFF=VRAI; 'FINS';
  'FINS';
  'SI' IKLFFF ;
    mocrit = 'CHAI' 'Fresidu' ;
    moflex = 'CHAI' 'MResidu' ;
  'FINSI' ;
*
* Jacobien du modele pour ponderer les champs de act3
  'SI' IGRD; 'FORM' GEOM1 ; 'FINSI';
  CH1  = 'MANU' 'CHML' zmodli 'SCAL' 1 'GRAVITE';
  zjac = 'INTG' 'ELEM' CH1  zmodli;
  zjac = zjac zmodli chan 'CHPO' 'SOMM' + xpetit;
  zjac = zjac **-0.5;
  'SI' IGRD; 'FORM' GEOREF0 ; 'FINSI';
*
* Initialisation CHArgement SANS T :
  CHASANST= WTAB.'CHARGEMENT' 'ENLE' 'T' ;
*
* isoucomp : Indicateur SOUci dans COMPortement, FAUX par defaut 
  isoucomp = faux ;
  PASUNIL  = FAUX;
  residc   = 'VIDE' 'CHPOINT' / 'DISCRET';
  RFNS     = 'VIDE' 'RIGIDITE' / 'RIGIDITE' ;
  RTRSF    = 'VIDE' 'RIGIDITE' / 'RIGIDITE' ;
  ZNITE    = 0 ;
  ISOUSPPP = WTAB.'ISOUSPAS' ;
  WTAB.'ISOUSPAS' = 0;
'FINSI' ;
*
RED_URG = 0;
resmul  = 1;
augmult = 0.30000000;
augm = 1.     ;
augk    = 1.;
XCONV   = 0.;
DEPSTDM = 0.;
 dpsmax = 0;
KNOCONV = 0 ;
WTAB.'CONV' = VRAI;
*----------------------------------------------------------------------*
*              Valeurs des champs en debut de pas TEMPS0               *
*----------------------------------------------------------------------*
'SI' WTAB.'FOR_MECA'  ;
  ZDEP0   = CONTI.'DEPLACEMENTS' ;
  ZDEF0   = 'REDU' CONTI.'DEFORMATIONS' ZMODL ;
  ZSIG0   = 'REDU' CONTI.'CONTRAINTES'  ZMODL ;
  com_sig = 'EXTR' ZMODLI 'CONT';
  'SI' IPLAVI ;
    ZDEI0 = 'REDU' CONTI.'DEFORMATIONS_INELASTIQUES' ZMODL ;
    ZVAR0 = 'REDU' CONTI.'VARIABLES_INTERNES'        ZMODL ;
    'SI' ('NEG' WTAB.'MOVA' 'RIEN') ;
      ACC0 = 'EXCO' ZVAR0 WTAB.'MOVA' ;
    'FINS' ;
    com_var = 'EXTR' ZMODLI 'VARI' ;
    com_dei = 'EXTR' ZMODLI 'DEIN' ;
    'SI' ISOL  ;
      com_maa = 'EXTR' ZMODLI 'MATE' ;
    'FINS' ;
  'FINS' ;
  ZGRDU0 = 'MOT' 'INCONNU';
  ZGRDUF = 'MOT' 'INCONNU';
  'SI' WTAB.'MEC1' ;
    'SI' (IGRD 'ET' WTAB.'UTILISATEUR');
      ZGRDU0 ='GRAD' ZMODL ZMAT20 ZDEP0 ;
    'FINS' ;
  'FINS' ;
'FINSI' ;
ZMAT2I = ZMAT20 ;
 
'SI' IFEFP ; 'FORM' GEOM1; FINSI;
*######################################################################*
*----------------------------------------------------------------------*
*       Boucle de non convergence BONOCONV (Configuration GEOM1)       *
*----------------------------------------------------------------------*
'REPETER' BONOCONV WTAB.'MAXSOUSPAS' ;
*
* Initialisation objectif de fin de pas. Forces et deplacements
  ZFCONSTA = ZFEXT2 ;
*
  excfconv = 0;
  'SI' ('EGA' RED_URG 0);
    'SI' ((XCONV < ZPREC) 'ET' (DEPSTDM < ZPRECeps) 'ET' (resmul > 0.99) 'ET' (AUGMULT < 100.));
      augmult = 0.30000000;
      augm = 1.     ;
      augk    = 1.;
      iraug = faux;
*pv resmul = 1;
 
      'SI' autaug ; 
        RE_CNTRL = 'VIDE' 'RIGIDITE' / 'RIGIDITE';
        DE_CNTRL  = FAUX;
        'SI' (dpsmax  < zprec);
          iraug = faux;
        'FINSI';
      'FINSI';
      znacce=VRAI;
      'SI' (IFEFP 'OU' ISSTE) ; ZNACCE = FAUX; 'FINS';
    'SINON';
      augmult = augmult / 1.8;
      si (augmult < 1d-3); augmult = 1d-3; finsi;
    'FINSI';
  'SINON';
    augmult = augmult * 1.5;
    'SI' (autaug 'ET' (augmult > 1d-2)); IRAUG = VRAI; 
*       'MESS' 'IRAUG a vrai';
    'FINSI';
  'FINSI';
*
  'SI' WTAB.'FOR_MECA' ;
    KNOCONV = KNOCONV+1 ;
*
    ZD0SLX = 'ENLE' ZDEP0 'LX';
    ZDEPT  = ZDEP0 '-' ZD0SLX ;
    DEPST  = 'ZERO' ZMODL 'DEFORMATIONS' ;
    ZSIGF  = ZSIG0 ;
*
*   Transporter le champ de contrainte de GEOREF0 -> GEOM1
    ZSIG1 = ZSIG0 ;
    'SI' IGRD ;
      'FORM' GEOM1;
      ZSIG1 = 'CONF' ZMODL ZSIG0 ;
    'FINSI' ;
*  
*   --------------------------------------------------------------------
*   Reevaluer les materiaux en convergence forcee
    'SI' (knoconv '>' 1) ;
      HPP_EPS = FAUX ;
      ZMAT21 = ZMAT20 ;
      'SI' IGRD;
        ZMAT21 = 'CONF' ZMODL ZMAT21 ;
      'FINSI';
      ZMAT1 = ZMAT21 ;
    'FINS' ;
*  
*   --------------------------------------------------------------------
*   Rigidite a la fin du pas 
    'SI' ('OU' WTAB.'RECARI' IRAUG ('NON' ('EXIS' WTAB 'RRRR')));
*
      'SI' ('OU' IENDOM IVIDOM ICERAM);
        HOOKENDO = 'HOOK' ZMODL ZMAT21   ZVAR0  ;
        RH       = 'RIGI' ZMODL HOOKENDO ZMAT21 ;
        'DETR' HOOKENDO;
      'SINON';
        'SI' (IGRD 'ET' ('EGA' LAG_TOT 1)) ;
          'FORM' GEOREF0;
          HOOKRH  = 'HOOK' ZMODL ZMAT21 ;
          'FORM' GEOM1;
          HOOKRH2 = 'CONF' HOOKRH ZMODL ;
          RH = 'RIGI' HOOKRH2 ZMODL ZMAT21;
          'DETR' HOOKRH; 'DETR' HOOKRH2;
        'SINON';
          RH  = 'RIGI' ZMODL ZMAT21 ;
        'FINSI';
      'FINS';
*
*     RH peut contenir des CL
      ZCL = 'EXTR' RH 'RIGI' 'MULT' ;
      'SI' ('NEG' ('DIME' ZCL) 0) ;
        ZCLIM = ZCLIM0 'ET' ZCL ;
      'FINSI' ;
      RRRR = RH 'ET' ZCLIM0 ;
*
*     Option 'RIGIDITE_CONSTANTE'
      'SI' IRCON;
        RRRR = RRRR 'ET' RIG_CONS;
      'FINS';
*
*     Option 'RIGIDITE_AUGMENTEE'
      'SI' IRAUGLU ;
        'SI' ('NON' (AUTAUG 'ET' IRAUG)) ;
          RRRR = RRRR 'ET' RIG_AUG  ;
        'FINSI';
      'FINSI';
*  
*     Option 'AUGMENTATION_AUTOMATIQUE'
      'SI' (AUTAUG 'ET' ('NON' IRAUGLU));
        RIG_AUG = 'MASSE' ZMODL ZMAT21 ;
***     mess 'actualisation rig_aug';
      'FINSI';
      'SI' (AUTAUG 'ET' IRAUG);
        RRRR = RRRR 'ET' (RIG_AUG * (augm * augmult)) 'ET' (RH   * (augk * augmult)) ;
**      'MESS' 'multiplicateur d augmentation masse' (augm * augmult) ' raideur' (augk * augmult) ;
      'FINS';
*  
*     Advection mecanique
      'SI' WTAB.'ADVECTION_MECANIQUE' ;
         RRRR = RRRR 'ET' ('ADVE' WTAB.'MOMEC_ADV' ZMAT21) ;
      'FINSI' ;
*
*     Stockage de la rigidite pour eviter de la recalculer
      'SI' ('EGA' ('DIME' ZCL) 0) ;
        WTAB.'RRRR'=RRRR;
      'FINSI';
*
    'FINS';
    'DETR' ZRAID;
    ZRAID = WTAB.'RRRR' ;
 
  'FINSI' ;
*FOR_MECA 
*
  'SI' WTAB.'NVSTNL'                                        ;
     'SI' ('NON' ('EXISTE' WTAB 'RNSL'))                    ;
   'SI' (PASDT '&lt;EG' 0.) ; 'MESS' 'PAS DE TEMPS !! ' PASDT  ;
                                        'ERRE' 21 ; 'FINSI' ;
        WTAB.'MAT_NSL' = WTAB.'MAT_NSL' 'ET' 
           ('MATE' ZMODL 'FREQ' (PASDT **(-1))) ;
        RINSL = 'RIGI' ZMODL WTAB.'MAT_NSL'        ;
        ndiv = 'RIGI' 'DIVE' ZMODL WTAB.'MAT_NSL'  ;
        WTAB.'NSMA' = 'MASSE' ZMODL WTAB.'MAT_NSL' ;
        WTAB.'RNSL' = RINSL 'ET' ZCLIM0 'ET' ndiv           ;
     'FINSI'                                                ;
  'FINSI'                                                   ;
*
*------------ consolidation ou dynamique faut-il recalculer l'operateur?
  'SI' (ISOL 'OU' IDYN);
    'SI' ( '>' (WTAB.'DTPREC' '*' 0.9999) PASDT) ;
      WTAB.'RECAOP' = VRAI;
    'FINS';
    'SI' ( '<' (WTAB.'DTPREC' '*' 1.0001) PASDT) ;
      WTAB.'RECAOP' = VRAI;
    'FINS';
    'SI' WTAB.'MATVAR';
      WTAB.'RECAOP' = VRAI;
    'FINS';
    'SI' ('NON' WTAB.'RECAOP') ;
      'SI' ('NEG' WTAB.'OPERATEUR' 'INCONNU');
        ZRAID = WTAB.'OPERATEUR';
      'FINS';
    'FINS';
  'FINS';
*
*------------ operateur amortissement en frequentiel
  'SI' WTAB.'FREQUENTIEL' ;
    RRR2 = 'AMOR' ZMODL ZMAT21 ;
    RR2  = 'CHAN' 'INCO' RRR2
            ('MOTS' 'ALFA' 'BETA') ('MOTS' 'IALF' 'IBET')
            ('MOTS' 'FALF' 'FBET') ('MOTS' 'FALF' 'FBET') 'QUEL' ;
    RR3  = 'CHAN' 'INCO' RRR2
              ('MOTS' 'ALFA' 'BETA') ('MOTS' 'ALFA' 'BETA')
              ('MOTS' 'FALF' 'FBET') ('MOTS' 'IFAL' 'IFBE') 'QUEL' ;
    RRR2 = RR2 'ET' RR3 ;
 
    RR1  = ZRAID ;
    OMEGI= 2.* PI * TI ;
    RRR1 = OMEGI * OMEGI * (-1.) * WTAB.'MASSE' ;
    RR1  = ZRAID 'ET' RRR1 ;
    RR4  = 'CHAN' 'INCO' (RR1 '*' (-1.D0))
            ('MOTS' 'ALFA' 'BETA') ('MOTS' 'IALF' 'IBET')
            ('MOTS' 'FALF' 'FBET') ('MOTS' 'IFAL' 'IFBE') 'QUEL';
    ZRAID= RR1   'ET' RR4  ;
    RR5  = OMEGI '*'  RRR2 ;
    ZRAID= ZRAID 'ET' RR5  ;
  'FINS' ;
 
*--------------- et la perméabilité  ----------------------------------
  'SI' ISOL ;
    'SI' (IGRD 'OU' WTAB.'MATVAR') ;
      WTAB.'PERMEABILITE' = 'PERM' WTAB.'MOD_POR' ZMAT21;
      WTAB.'RECAOP'       =  VRAI ;
    'FINS';
  'FINS';
 
*------------- Cas de la consolidation ou de la dynamique -------------
*------------- il faut recalculer l'operateur d'iteration -------------
  'SI' WTAB.'RECAOP' ;
*
    'SI' IDYN;
      ZRMAS = 4.D0 '/' (PASDT '**' 2) '*' WTAB.'MASSE' ;
      ZRAID = ZRMAS 'ET' ZRAID;
      'SI' ('NEG' WTAB.'AMORTISSEMENT' 'INCONNU');
        ZRAID = WTAB.'AMORTISSEMENT' '*' (2.D0 '/' PASDT ) 'ET' ZRAID;
      'FINS';
    'FINS' ;
*
    'SI' ISOL ;
      ZRAID =-1.* PASDT * WTAB.'TETA'* WTAB.'PERMEABILITE' 'ET' ZRAID ;
    'FINS' ;
*
    WTAB.'OPERATEUR'= ZRAID ;
  'FINS';
*
  'SI' WTAB.'FOR_MECA' ;
*
*   --------------------------------------------------------------------
*   Conditions de contact-frottement
    'SI' WTAB.'CONTACT';
      MODCON = WTAB.'MODCONTA';
      'SI' ('EXIS' WTAB 'MATCONTA') ;
        CJEU CRR = 'RFCO' MODCON WTAB.'CONV' WTAB.'MATCONTA' ;
      'SINON' ;
        CJEU CRR = 'RFCO' MODCON WTAB.'CONV'                 ;
      'FINS';
*
*     Separer ce qui concerne l'adherence (FADH) et les jeux (FLX)
      'SI' ('NEG' CJEU 0) ;
        FADHE = 'EXCO' CJEU 'FADH' 'NOID' 'FADH';
        CDAP  = 'EXCO' CJEU 'FLX'  'NOID' 'FLX' ;
      'FINSI' ;
*
      'SI' WTAB.'MODAL' ;
        CRR  = 'PJBA' CRR   ZMODL ZMAT21 ;
        MCRR = 'EXTR' CRR  'MAIL' 'MULT' ;
        MCDAP= 'EXTR' CDAP 'MAIL'        ;
        'REPETER' BCDA ('NBNO' MCDAP) ;
          PBCDA ='POIN' MCDAP &BCDA ;
          PCRR  ='POIN' MCRR  &BCDA ;
          CHCR  ='MANU' 'CHPO' PCRR 1 'FLX' ('EXTR' CDAP 'FLX' PBCDA)
                 'NATURE' 'DISCRETE' ;
          'SI' ('EGA' 1 &BCDA) ;
            CCDA = CHCR ;
          'SINON' ;
            CCDA = CHCR 'ET' CCDA ;
          'FINS' ;
        'FIN' BCDA ;
        CDAP = CCDA ;
      'FINS' ;
*
      ZFCONT = 'VIDE' 'CHPOINT' / 'DISCRET' ; 
      'SI' ('NEG' CRR 0) ;
        ZFCONT = CRR '*' ZD0SLX 'EXCO' 'FLX' 'NOID' 'FLX';
        ZCLIM  = CRR 'ET' ZCLIM ;
        ZRAID  = CRR 'ET' ZRAID ;
      'FINS';
*
      'SI' ('NEG' CJEU 0) ;
        ZFCONT = ZFCONT '+' CDAP ;
      'FINS';
*
      'SI' (('NEG' CRR 0) 'OU' ('NEG' CJEU 0)) ;
        ZFCONSTA = ZFCONSTA '+' ZFCONT ;
      'FINS';
*
    'FINS';
*
*   --------------------------------------------------------------------
*   Matrice tangente et fefp
    'SI' IKTAN ;
      'SI' IFEFP ;
        IKT_SAUV = VRAI ;
        'SI' ('NEG' WTAB.'LASTKTAN' 'INCONNU') ;
          'MESS' 'FEFP: Start with LASTKTAN' ;
          ZRIKTA = ZLASTKTAN ;
          ZRAID = ZCLIM 'ET' ZRIKTA ;
        'SINON' ;
          'MESS' 'FEFP: Previous KTAN not available' ;
          ZRAID = ZRAID 'ET' ('KSIGMA' ZMODL ZSIG0 ZMAT21) ;
        'FINS' ;
      'SINON' ;
        'SI' ('NEG' WTAB.'LASTKTAN' 'INCONNU') ;
          IKT_SAUV = VRAI ;
          'SI' IPERT ;
            'MESS' 'Matrice tangente par perturbation - '
            'Demarrage avec KTAN = LASTKTAN' ;
          'SINON' ;
            'MESS' 'Matrice tangente "coherente" - '
            'Demarrage avec KTAN = LASTKTAN' ;
          'FINS' ;
          ZRIKTA = ZLASTKTAN ;
          ZRAID = ZCLIM 'ET' ZRIKTA ;
        'SINON' ;
          'SI' IPERT ;
            IKT_SAUV = 'NEG' WTAB.'K_TANGENT_ITER0' 'MAT_ELASTIQUE';
            'MESS' 'Matrice tangente par perturbation - '
            'Demarrage avec KTAN = rigidite elastique' ;
          'SINON' ;
            IKT_SAUV = ('NEG' WTAB.'K_TANGENT_ITER0' 'MAT_ELASTIQUE')
                  'ET' ('NEG' WTAB.'K_TANGENT_ITER0' 'MAT_TANGENTE') ;
            'SI' ('EGA' WTAB.'K_TANGENT_ITER0' 'MAT_ELASTIQUE') ;
              'MESS' 'Matrice tangente "coherente" - '
                              'Demarrage avec KTAN = rigidite elastique' ;
            'SINON' ;
              'MESS' 'Matrice tangente "coherente" - '
                                 'Demarrage avec KTAN (DTTAN = 0.)' ;
              ZRIKTA = 'KTAN' ZMODL ZSIG0 ZVAR0 ZMAT21
                              'PREC' ZPREK 'DT  ' DTTAN ZKTASYM ;
              ZRAID = ZCLIM 'ET' ZRIKTA ;
            'FINS' ;
          'FINS' ;
        'FINS' ;
      'FINS' ;
    'FINS' ;
    'SI' WTAB.'MAN' ; ZRAIDINI = ZRAID ; 'FINSI' ;
*
*   --------------------------------------------------------------------
*   Raideur geometrique (option K_SIGMA)
    'SI' (IGRD 'ET' ('NON' HPP_EPS)) ;
      'SI' (IKSIA 'ET' ('NON' IFEFP)) ;
        KSIGTC = 'KSIGMA' ZMODL ZSIG1 ZMAT21;
        ZRAID  = ZRAID 'ET' KSIGTC;
      'FINS' ;
    'FINS' ;
*----------------------------------------------------------------------*
*                  Pilotage automatique                                *
*----------------------------------------------------------------------*
    COEPI = 1.D0 ;
*
    'SI' IPILOT ;
      'SI' WTAB.'AUTODEUX' ;
        COEPI = WTAB.'AUTOCOEF' ;
        COEPI = 'MINI' (COEPI / (1. '-' COEPI)) 1.D0 ;
      'FINS' ;
      ZAUTOREDU = 1.D0 ;
      RED2 = 0 ;
    'FINS';
*
  'FINSI' ;
*FOR_MECA 
*----------------------------------------------------------------------*
*             Partie variable du chargement a imposer a TI             *
*          (geometrie, materiaux, chargement en deformations)          *
*----------------------------------------------------------------------*
  DEPSI = 0. ; EPSIF = 0. ; 
  DSIGM = 0. ; SIGFI = 0. ; 
  DFDEF = 0. ; FSIGF = 0. ; 
  ZFSUIV= 0. ;
*
* - Cas de materiaux a proprietes variables
  'SI' WTAB.'MATVAR' ;
    'SI' ('OU' IENDOM IVIDOM ICERAM);
      XXX3 = 'ELAS' ZMODL ZSIG1 ZMAT20 ZVAR0;
      XXX4 = 'ELAS' ZMODL ZSIG1 ZMAT1  ZVAR0;
    'SINON';
      XXX3 = 'ELAS' ZMODL ZSIG1 ZMAT20      ;
      XXX4 = 'ELAS' ZMODL ZSIG1 ZMAT1       ;
    'FINS';
    DDEF0 = XXX4 - XXX3;
    'DETR' XXX3; 'DETR' XXX4;
    DEPSI = -1.* DDEF0;
  'FINS';
*
* - Chargement thermique
  'SI' ITHER ;
    DTETD = ZTEMP2 '-' ZTEMP1 ;
    ETT   =('EPTH' ZMODL ZMAT20 MCHTETA2 ) - ('REDU' WTAB.'ETREF' ZMODL) ;
    DTT   = ETT '-' ETT0 ;
    DEPSI = DEPSI '+' DTT ;
    EPSIF = EPSIF '+' ETT ;
  'FINS';
*
* - Chargement en deformations imposees
  'SI' LOGDEF;
    DDEFO = ZDEFOR2 '-' ZDEFOR1 ;
    DEPSI = DEPSI '+' DDEFO ;
    EPSIF = EPSIF '+' ZDEFOR2 ;
  'FINS';
*
* - Contraintes associees a ces deformations
  'SI' (DEPSI 'NEG' 0);
    'SI' ('OU' IENDOM IVIDOM ICERAM);
      DSIGM = 'ELAS' ZMODL DEPSI ZMAT20 ZVAR0;
      'SI' (EPSIF 'NEG' 0) ;
        SIGFI = 'ELAS' ZMODL EPSIF ZMAT20 ZVAR0;
      'FINSI' ;
    'SINON';
      DSIGM = 'ELAS' ZMODL DEPSI ZMAT20 ;
      'SI' (EPSIF 'NEG' 0) ;
        SIGFI = 'ELAS' ZMODL EPSIF ZMAT20 ;
      'FINSI' ;
    'FINS' ;
  'FINS';
 
* - Cas d'un milieu poreux avec chargement thermique
* Cas isotrope seulement et on ne s'occupe pas du alpha-reference !!
  'SI' (ITHER 'ET' WTAB.'POR1') ;
    MSRTT  = PAS_EPTH PRECED ZMODL ZMAT20 MCHTETA2 ;
    MSRTT1 = 'REDU' MSRTT ZMODL ;
    DMSRT0 = MSRTT '-' MSRTT0 ;
*
    DSIGM  = DSIGM '+' DMSRT0 ;
    SIGFI  = SIGFI '+' MSRTT  ;
  'FINS' ;
*
* - Forces equivalentes a ces contraintes
  'SI' (DSIGM 'NEG' 0);
    DSIGT = DSIGM ;
    'SI' IGRD ; 
      DSIGT = 'CONF' ZMODL DSIGM ; 
    'FINSI' ;
    DFDEF = 'BSIG' ZMODL DSIGT ZMAT21 ;
    'SI' (SIGFI 'NEG' 0) ;
      SIGFT = SIGFI ;
      'SI' IGRD ; 
        SIGFT = 'CONF' ZMODL SIGFI ; 
      'FINSI' ;
      FSIGF = 'BSIG' ZMODL SIGFT ZMAT21 ;
    'FINSI' ;
  'FINS';
*
* - Chargement de pression suiveuse en grands deplacements
  'SI' (LOGPRE 'ET' IGRD) ;
    ZPEXTF = 'TIRE' WTAB.'CHARGEMENT' 'PRES' TI ;
    'SI' ('EXIS' WTAB 'MAT_PRE') ;
      ZFPEXTF = 'BSIG' WTAB.'MOD_PRE' ZPEXTF WTAB.'MAT_PRE' ;
    'SINON' ;
      ZFPEXTF = 'BSIG' WTAB.'MOD_PRE' ZPEXTF ;
    'FINS' ;
    ZFSUIV = ZFSUIV '+' ZFPEXTF ;
  'FINS';
*----------------------------------------------------------------------*
*     Procedure utilisateur : y a-t-il des forces non conservatives    *
*----------------------------------------------------------------------*
  'SI' WTAB.'PROCEDURE_CHARMECA';
*   on ajoute l indice ADDI_MATRICE pour signaler a charmeca qu on
*   souhaite aussi l operateur linearisé des Forces NL de charmeca
    PRECED.'ADDI_MATRICE' = vrai;
    TFP22  = CHARMECA PRECED TI ;
    PRECED.'ADDI_MATRICE' = faux;
*
*   FP22 = F^suiv_n+1
    ADDISEC2 = 'EXIS' TFP22 'ADDI_SECOND' ;
    'SI' ADDISEC2 ;
      FP22   = TFP22.'ADDI_SECOND' ;
      FP022  = 'COPIER' FP22       ;
      ZFSUIV = ZFSUIV '+' FP22     ;
    'FINS';
    'SI' ('EXIS' TFP22 'ADDI_MATRICE');
      ZRAID = ZRAID 'ET' TFP22.'ADDI_MATRICE';
    'FINS';
  'FINS';
*----------------------------------------------------------------------*
*                      Second membre    RESIDU                         *
*                                                                      *
* Calcul du premier residu : desequilibre entre les forces externes et *
* le calcul B*SIGMA. Le sigma qui sert est celui qui existerait si     *
*   le champ de deplacement ne changeait pas (ZDEP0).                  *
*                                                                      *
*          | forces exterieures sans reactions - forces interieures |  *
*          |   F^ext_n+1 + DF^suiv + DF^ther + DF^defi - F^int_n    |  *
*          |                         |-------------------------|    |  *
*          |                                  -1.*ZFPLO             |  *
*          |          (composantes de forces FX FY FZ ...)          |  *  
* RESIDU = |                                                        |  *
*          |      increment des relations imposees Du^imp           |  *
*          |          (composantes de depl. FLX)                    |  *
*                                                                      *
*  A F^ext peuvent s'ajouter des termes supplementaires p.ex. en       *
*   dynamique ou en poreux.                                            *
*----------------------------------------------------------------------*
* A-t-on des C.L. unilaterales ?
  MCLIMU = 'EXTR' ZCLIM 'MAIL' 'UNIL';
  IMPO12 = ('NBEL' MCLIMU) '>' 0 ;
*
* ZFCONSTA = [ F^ext_n+1 ; u^imp_n+1 ]
  ZFEXT = 'ENLE' ZFCONSTA 'FLX';
  ZFLX1 = 'EXCO' ZFCONSTA 'FLX' 'NOID' 'FLX' 'NATURE' 'DISCRET';
*
  'SI' WTAB.'FOR_MECA' ;
*
*   Forces externes deja equilibrees au debut du pas par B*sigma
*   ZF1 = F^int_n = B*sigma_n + K^cst*u_n
*HHO : Modifications appel a BSIGMA (ajout du champ de deplacements necessaire en HHO)
    ZF1 = 'BSIG' ZMODL ZSIG1 ZMAT1 ZDEP0 ;
    'SI' IRCON;
      ZF1 = ZF1 '+' ('REDU' MAI_CONS (RIG_CONS '*' ZDEP0));
    'FINS';
*
    'SI' IDYN ; 
      FFDYN = 'COPIER' ZF1; 
    'FINS';
*
    'SI' ISOL ;
      GRAP0= 'GRAD' WTAB.'MOD_POR' ZDEP0 WTAB.'MAT1' 'CONS' ;
      XXX1 = 'GRAD' WTAB.'MOD_POR' ZDEP0 MA_POR      'CONS' ;
      XXXS =((1.- WTAB. 'TETA' )*GRAP0)+ (WTAB. 'TETA' * XXX1);
      XXX2 = PASDT '*' ('GNFL' WTAB.'MOD_POR' XXXS) ;
      XXX3 = ZF1  ;
      ZF1  = XXX3 - XXX2;
      'DETR' XXX3; 'DETR' XXX2 ;
    'FINS';
*
*   [ FLXINI ; FREAP ] = [ u^imp_n ; -F^reac_n ]
    FZDEP0 = ZDEP0 '*' ZCLIM;
    FLXINI = 'EXCO' FZDEP0 'FLX' 'NOID' 'FLX' 'NATURE' 'DISCRET';
    FLXREC = FLXINI ;
    FREAP  = 'ENLE' FZDEP0 'FLX';
*
*   XXX1 = [ F^ext_n+1 ; u^imp_n+1 - u^imp_n ]
    XXX1 = ZFCONSTA '-' FLXINI ;
*
*   FEXT0 : chargement externe (sans reactions) au pas precedent
*   FEXT0 = F^int_n - F^reac_n = F^ext_n
    FEXT0 = ZF1 '+' FREAP;
*
*   RESIDU : forces exterieures sans reactions (avec des termes 
*   supplementaires le cas echeant p.ex. en dynamique ou en poreux)
*   - forces interieures et increment des relations imposees
*
*   RESIDU = [DF^tot ; Du^imp]
*   RESIDU = [F^ext_n+1 + F^suiv_n+1 + DF^ther + DF^defi - F^int_n ; Du^imp]
    RESIDU = XXX1 '+' DFDEF '-' ZF1;
    RESIDU = RESIDU '+' ZFSUIV ;
*
*   ZFPLO  = -1.* [- F^int_n + DF^ther + DF^defi ]
    ZFPLO  = ZF1 '-' DFDEF ;
*
*   ZDFORC = [F^ext_n+1 + DF^ther + DF^defi - (F^int_n - F^reac_n) ; Du^imp]
    ZDFORC = XXX1 '-' FEXT0 ;
*
*   Inc. de forces et de deplacements en distinguant CL unil et autres
    ZLX     = 'EXCO' ZDEP0  'LX' 'NOID' 'LX' 'NATURE' 'DIFFUS';
    DFEXT0F = 'ENLE' ZDFORC 'FLX';
    DFEXT0L = 'EXCO' ZDFORC 'FLX' 'NOID' 'FLX';
    'SI' IMPO12 ;
      DUUNIL = 'REDU' DFEXT0L MCLIMU ;
      XXX2   = 'DIFF' ('EXTR' ZCLIM 'MAIL' 'MULT') MCLIMU ;
      DUIMPO = 'REDU' DFEXT0L XXX2 ;
    'FINS';
*
    'SI' ('EGA' KNOCONV 1) ;
      ZDFINI = 1. '*' ZDFORC ;
      'SI' (DFDEF 'NEG' 0.) ;
        ZDFINI = ZDFINI '+' DFDEF;
      'FINS';
    'FINS';
*
    ZSDMBR = RESIDU '*' 1.D0;
*
  'FINSI' ;
* FOR_MECA
*
  'SI' WTAB.'NVSTNL'                          ;
    ZVIFL = conti.'VITESSES_FLUIDE'           ;
    movifl = 'EXTR' ZMODL 'DEPL'     ;
    UN = 'EXCO' ZVIFL movifl                  ;
    ch_un = 'MANU' 'CHML' ZMODL 'ADVE' un 'TYPE' 'CARACTERISTIQUES';
    nugrad = 'ADVE' ZMODL ch_un      ;
    'SI' ('EXIS' PRECED 'OPTION_BDF2') ;
      ZRNSL = WTAB.'RNSL' 'ET' nugrad 'ET' (WTAB.'NSMA' '*' 1.5)  ;
      unm = 'EXCO' conti.'VITESSES_FLUIDE_0' movifl ;
      fmass = WTAB.'NSMA' '*' ((2. '*' un) '-' (0.5 '*' unm))  ;
    'SINON' ;
      ZRNSL = WTAB.'RNSL' 'ET' nugrad 'ET' WTAB.'NSMA' ;
      fmass = WTAB.'NSMA' * un                ;
    'FINSI' ;
    ftnsl  = ZFCONSTA 'ET' fmass              ;
 
    'SI' LOGADV                               ;
      ch_ad = 'MANU' 'CHML' ZMODL 'ADVE' UADV 'TYPE' 'CARACTERISTIQUES';
      nuadv = 'ADVE' ZMODL ch_ad     ;
      ZRNSL = ZRNSL 'ET' nuadv                ;
      fadv = nugrad * UADV                    ;
      mess 'advec / inertie' (maxi (abs fadv)) (maxi (abs ftnsl)) (maxi (abs UADV)) ;
      ftnsl = ftnsl '-' fadv                  ;
    'FINSI'                                   ;
 
    resnsl = ftnsl '-' (ZRNSL * ZVIFL)        ;
 
* resolution 
    ZVIF1  = 'KRES' ZRNSL resnsl              ; 
    ZVIFL = ZVIFL '+' ZVIF1                   ;
  'FINSI'                                     ;
*----------------------------------------------------------------------*
*           1ere Resolution  - ZDEP1 = [ Du^0 ; LX_n+1 ]               *
*                                                                      *
*           -> estimation a TI a partir de l'instant TDEBUT            *
*                                                                      *
*             K*Du^0 + At*LX_n+1 = F^ext_n+1 - Bsigma_n                *
*             A*Du^0             = d1 - A*u_n                          *
*             |----------------|  |-------------------|                *
*               ZRAID * ZDEP1    =       RESIDU                        *
*----------------------------------------------------------------------*
  'SI' WTAB.'FOR_MECA' ;
    FEXCI = 'VIDE' 'CHPOINT' / 'DISCRET' ;
    'SI' WTAB.'ADHERENCE' ;
      'SI' ('EXIS' WTAB 'MAIL_BLOM') ;
        FADRE = 'REDU' FADHE WTAB.'MAIL_BLOM' ;
      'SINON' ;
        FADRE = FADHE ;
      'FINSI' ;
      FEXCI = FEXCI 'ET' FADRE ;
    'FINSI' ;
*  
*   Modification du residu si automatique
    'SI' LOGPIL ;      
      RESIDU = RESIDU 'ET' (ETA0 '*' ZFPILIN);        
    'FINSI' ;
*  
*   Force limite de frottement
    'SI' (WTAB.'CAFROTTE' 'ET' IMPO12);
      'SI' WTAB.'FROCABL' ;
        zsigfT = 'REDU' zsigf ZMODLI;
        FFROT = 'EXCC' ZCLIM ZDEPT MODCON WTAB.'MATCONTA' ZSIGFT;
        RTRSF = 'VIDE' 'RIGIDITE' / 'RIGIDITE' ;
        RFNS  = 'VIDE' 'RIGIDITE' / 'RIGIDITE' ;
        FEXCI = FEXCI '+' FFROT ;
      'FINSI';
      'SI' WTAB.'FROCOUL' ;
        excfconv = excfconv + 1 ;  
        RTRSF*'RIGIDITE' RFNS*'RIGIDITE'  FFROT = 'EXCF' ZCLIM ZDEPT MODCON WTAB.'MATCONTA'  excfconv;
        FEXCI = FEXCI '+' FFROT ;
        MFROT = 'EXTR' FFROT 'MAIL' ;
      'FINSI';
    'FINS';
*  
    RESIDU = 'CHAN' 'ATTRIBUT' RESIDU 'NATURE' 'DISCRET' ;
    'SI' ('NEG' ZRIBLO_M 'INCONNU') ;
      'SI' ((WTAB.'CAFROTTE' 'OU' WTAB.'ADHERENCE') 'ET' IMPO12);
        ZDEP1 BID BID BID BID     = 'RESO' ZRAID 'SOUC' RESIDU 'INIB'
                                       ZRIBLO_M ZLISEA_M FEXCI RFNS RTRSF ;
      'SINON';
        ZDEP1 BID BID BID     = 'RESO' ZRAID 'SOUC' RESIDU 'INIB'
                                              ZRIBLO_M ZLISEA_M ;
      'FINS';
      ZRIBLO_M = 'MOT' 'INCONNU' ;
    'SINON';
      'SI' ((WTAB.'CAFROTTE' 'OU' WTAB.'ADHERENCE') 'ET' IMPO12);
        ZDEP1 = 'RESO' ZRAID 'SOUC' RESIDU FEXCI RFNS RTRSF          ;
      'SINON';
        ZDEP1 = 'RESO' ZRAID 'SOUC' RESIDU;
      'FINS';
    'FINS';
*  
    'SI' LOGPIL ; 
      ZDEPII  = 'RESO' ZRAID ZFPILIN ; 
      ZDEPILO = 0. * ZDEPII ;     
*     Calcul de D_eta par appel a la procedure PILOINDI
      D_ETA = PILOINDI PRECED ZDEP0 ZDEPILO ZDEP1 ZDEPII DTAU; 
*     Mise à jour  
      ZDEP1 = ZDEP1 '+' (D_ETA '*' ZDEPII) ;       
      ETA = ETA0 ;  
    'FINSI' ;
*
*   Pour initialiser le statut des CL unil lors du prochain appel a RESO
    ZRAID_T = 'EXTR' ZRAID 'CONT';
    'SI' ('EXIS' ZRAID_T 'NITER');
*     Conserver les conditions unilaterales pour VITEUNIL
      'SI' IDYN;
        WTAB.'ZRAIDV' = ZRAID;
      'FINS';
      ZRIBLO_M = ZRAID_T. 7 ;
      ZLISEA_M = ZRAID_T. 6 ;
    'FINS';
*
*   Initialisation zdepl (sert pour xnum)
    ZDEPL = 'EXCO' ZDEP1 'LX' 'NOID' 'LX' ;     
*----------------------------------------------------------------------*
*              Calcul d'une norme pour la convergence                  *
*----------------------------------------------------------------------*
    XXX1 = ZFEXT '+' ZFSUIV ;
    'SI' LOGPIL ;      
      XXX1 = XXX1 '+' ((ETA '+' D_ETA) '*' ZFPILIN);  
    'FINSI' ;
*
    'SI' ('EGA' knoconv 1) ;
      ZDEP1P50 = ZDEP1 ;
      XDENO='XTY' ZDEP1P50  ( XXX1 -( RESIDU  'EXCO'
          'FLX' 'NOID' 'FLX' 'NATURE' 'DISCRET')) MLPRIM MLDUAL;
      'SI' ('VERI' xdeno) ; 'SINON'; xdeno = 1; 'FINSI';
      MZDEP1M = 'MAXI' ZDEP1P50 'ABS' 'AVEC' MLDEPL;
      MZFM = 'MAXI' (FSIGF + DFDEF + ZF1) 'ABS' 'AVEC' MLDUAL;
      XDENO1 = 'ABS' XDENO  + (MZFM * MZDEP1M);
      MZDEP1M = MZDEP1M + XPETIT;
      XDENO=XDENO1/MZDEP1M;
      XDENO = XDENO + MZFM;
*     XDENO = 'MAXI' XDENO ('MAXI' XXX1 'ABS' 'AVEC' MLDUAL);
      XDENO =            'MAXI' XDENO ('MAXI' XXX1 'ABS' 'AVEC' MNDUAL);
       'SI' IGRD; XDENO = 'MAXI' XDENO ('MAXI' (RH * zdep0) 'ABS' 'AVEC' MNDUAL); 'FINSI';
      'SI' (XDENO < XPETIT); XDENO = 1.; 'FINSI';
      XDENOM=XDENO;
      'SI' TSTMOM ;
        ZDEP1P50 = ZDEP1 + XPETIT;          ;
        XDENOM = XDENO1/('MAXI' ZDEP1P50 'ABS' 'AVEC' MLROTA);
        XDENOM = XDENOM + XPETIT       ;
      'FINS' ;
    'FINS' ;
*
**  mess 'xdeno propose' ' ' xdeno ' xdeno precedent' ' ' xdenoo;
**     XDENO = MAXI XDENOo XDENO;
**     XDENOo = XDENO;
**     XDENOM = MAXI XDENOMo XDENOM;
**     XDENOMo = XDENOM;
*
    'SI' WTAB.'CONV';
      ZINCREMENT = ZDFINI '+' ZFSUIV ; 
    'FINS';
*
    RESIDNOR = 'COPIER' RESIDU ;
*
    'SI' IPILOT;
*     Objectif non atteint : conserver le XDENO
      'SI' WTAB.'AUTODEUX' ;
        XDENO  = ZXDENO ;
        XDENOM = ZXDENOM ;
      'FINS';
*
*     Forces/deplacements en fin de pas
      DFEXT = COEPI '*' DFEXT0F ;
      ZFEXT = DFEXT '+' FEXT0 ;
      'SI' IMPO12;
        DUIMP = (COEPI '*' DUIMPO) '+' DUUNIL ;
      'SINON' ;
        DUIMP = COEPI '*' DFEXT0L ;
      'FINS';
      ZFLX1 = DUIMP '+' FLXINI ;
*
      XFORC  = ZFEXT '+' (COEPI '*' ZFSUIV) '-' ZF1 ;
      RESIDU = XFORC '+' DUIMP ;
    'FINS';
*----------------------------------------------------------------------*
*            Corriger le residu a partir du pas precedent              *
*----------------------------------------------------------------------*
*   petite correction du residu pour esperer gagner du temps ***********
    INIT = FAUX ;
    'SI'(('NEG' ZFNONL 'INCONNU') 'ET' WTAB.'INITIALISATION');
      'SI' IPILOT;
        'SI' (WTAB.'AUTODEUX' 'ET' (COEPI 'NEG' 1.D0)) ;
          'MESS' 'Initialisation a partir du pas precedent ';
          INIT   = VRAI ;
          RESIDU = RESIDU '+' ZFNONL ;
        'FINS' ;
      'SINON';
*       Faire la correction si le pas precedent a converge sans non convergence
        'SI' (WTAB.'CONV' 'ET' (ISOUSPPP 'EGA' 0));
*         Faire la correction si le pas precedent etait non lineaire.
          'SI' (('MAXI' 'ABS' ZFNONL) > (ZPREC * XDENO)) ;
* on enleve le residu du pas precedent pour recuperer l'increment nominal du 
* second membre a imposer et l'increment du second membre du pas precedent
            ZINCREMENT = ZINCREMENT - WTAB.'RESIDU' ;
            zdeps = WTAB.'ZDEP1' + zdep1;
            FFNO  = 'XTY' ZFNONL     zdeps MNDUAL MNPRIM;   
            f12   = 'XTY' ZINCREMENT zdeps MNDUAL MNPRIM;                
            f22   = 'XTY' INCRPREC   zdeps MNDUAL MNPRIM;                
            AMPL = f12/(f22 + XPETIT); AMPLT = 0;
            'SI' (WTAB.'DTPREC' '>' XPETIT);
              AMPLT = PASDT '/' WTAB.'DTPREC' ;
            'FINS';
*  Le chargement n'est il pas de fluage ou de thermique ?';
            XDCOMP = ('XTY' ZDEP1 ZINCREMENT MNPRIM MNDUAL) ;  
            'SI' (('ABS' XDCOMP) < (ZPREC * XDENO * mzdep1m)
                      'OU' (('ABS' FFNO) > (('ABS' F22) * 2.e2   )));
              AMPL=AMPLT;
* la decharge est-elle significative
              'SI'((F12/(f22 + XPETIT)) < -0.05);ampl=0.;'FINS';
**            'MESS' 'F12 F22_' f12 f22;
              'MESS' 'Pas d increment de charge, initialisation calculee avec le temps' ;
            'SINON';
              AMPL = F12 / (F22 + XPETIT);
            'FINS';
* changement de modele on n'initialise pas
            'SI' ('NEG' ZMODLI WTAB.'MO_TOT_PREC');  AMPL = 0; 'FINSI';
              AMPL = MINI (prog AMPL AMPLT);
              'SI' ((AMPL > 0) 'ET' (AMPL < 2e1)) ;
      'MESS' 'Initialisation a partir de la solution precedente Coeff'AMPL;
                XXX1 = AMPL * ZFNONL ;
                XXX2 =RESIDU+ XXX1;
                'DETR' XXX1;'DETR' RESIDU;
                RESIDU = XXX2;
                INIT = VRAI;
            'FINS';
          'FINS';
        'FINS';
      'FINS';
    'FINS';
    WTAB.'ZDEP1'=zdep1;
*
*   Initialisation pour l'acceleration de convergence
    ACFP1  = 'VIDE' 'CHPOINT';
    ACFP2  = ACFP1 ;
    ACFP3  = ACFP1 ;
    ACFEP1 = ACFP1 ;
    ACFEP2 = ACFP1 ;
    FCORF  = 'COPIER' FREAP ;
    CORREC = 0;
*
    zdept = (zdept 'EXCO' 'LX' 'NOID' 'LX') + (ZDEPT * 0.);
*----------------------------------------------------------------------*
*       Initialisation des messages pour le pas de temps courant       *
*----------------------------------------------------------------------*
    'SI' IPILOT ;
      'MESSAGE' 'Iter'*13 'Nplas'*26 mocrit*39 'Deps.max'*52 'Eps.max'*65 moflex*78 'Alpha'*91;
    'SINON' ;
      'MESSAGE' 'Iter'*13 'Nplas'*26 mocrit*39 'Deps.max'*52 'Eps.max'*65 moflex*78 ;
    'FINSI' ;
  'FINSI' ;
*FOR_MECA
*
  NONCONV = FAUX; ZICONV = VRAI; PASTEST = FAUX; PASREINI = VRAI;
  ITACC = 0 ; RECA_N = 0 ; IT_RECA = 0 ;
  RED_URG = 0 ; DPSMAXP = 1; DEPSTDM = 0. ; MMCMAX = 0 ;
  IT = 0 ; MMC = 0 ; XCONV = 0. ; DPSMAX = xpetit ; EPSM = 0. ;
  IPREM = VRAI ; IRATE = FAUX; URG = FAUX ;
  IKT = FAUX ; zprecnc=1e-5 ; TABCONV = 'TABL';
  RE_CNTRL = 'VIDE' 'RIGIDITE' / 'RIGIDITE'; DE_CNTRL = FAUX;
*######################################################################*
*----------------------------------------------------------------------*
*           Boucle de convergence ETIQ (Configuration GEOM1)           *
*----------------------------------------------------------------------*
* preparation pour la gestion du pas avant etiq
  depstp zdeptp zsigfp fcorfp correcp geom2p zmat22p = depst zdept zsigf fcorf 0 geom1 zmat21 ;
*
  ITURG = 0; TFIN = TI ;
  'REPETER' ETIQ ;
*
*   IT est le compteur de ETIQ
    IT = IT + 1 ;
*
    'SI' WTAB.'FOR_MECA' ;
*
      depst zdept zsigf fcorf correc geom2 zmat22 = depstp zdeptp zsigfp fcorfp correcp geom2p zmat22p ;
*
*     zdep1d doit etre coherent avec residu
      'SI' IPREM; 
        zdep1d = ZDEP0 'ENLE' 'LX';
      'SINON' ;
        'SI' ('NON' IRATE) ;
          zdep1d = ZDEP1 'ENLE' 'LX';
        'FINSI';
      'FINSI';
*
*     Cas ou l'iteration precedente ne s'est pas faite entierement
*     -> partir du bon etat (mis a jour si le sous increment est acceptable) 
      'SI' IRATE;
        HPP_EPS = FAUX ;
        URG     = VRAI ;
        ITACC   = 4 ;
        'SI' AUTAUG ; IRAUG = VRAI ; 'FINSI' ;
      'FINSI' ;
      IRATE = VRAI;
*              
      'SI' autaug;
        resmul = 2. * resmul;
        'SI' ((dpsmax <   1d-3) 'ET' (resmul > 1D-20) 'ET' ('NON' PASTEST)); 
          resmul = resmul / (dpsmax + xpetit) * 1d-3; 
        'FINSI';
        'SI' (resmul > 1.0); resmul = 1.0; 'SINON'; itacc = 4; 'FINSI';
        'SI' (resmul < 1.d-40); 
              RE_CNTRL = 'VIDE' 'RIGIDITE' / 'RIGIDITE';
              DE_CNTRL = FAUX;
        'FINSI';
      'FINSI';            
*
      zdep1 = (zdept 'EXCO' 'LX' 'NOID' 'LX') + (ZDEPT * 0.);
*
*     Pas de temps courant eventuellement sous-decoupe via COEPI
      ZDT = PASDT '*' COEPI ;
*
*     ITACC doit etre =< 0 pour qu'on accelere
      ITACC = ITACC - 1;
*
      PASTEST = FAUX;
      RECA_K  = FAUX;
*
*     Pour conserver les criteres de convergence
      tabconv. it = 1;
*----------------------------------------------------------------------*
*                Recalcul de la rigidite si necessaire                 *
*----------------------------------------------------------------------*
      'SI' IGRD ;
        'FORM' GEOM1;
*
        ITURG = ITURG + 1;
        'SI' ((IT > 1) 'ET' (URG 'OU' ('MULT' (IT - IT_RECA) ITRCLC))) ;
*
          URG     = FAUX ;
          ITURG = 0;
          IT_RECA = IT;
          RECA_K  = VRAI;
          RECA_N  = RECA_N + 1;
          HPP_EPS = FAUX;
          'SI' (RECA_N > 20) ;
            nonconv = vrai;
          'FINS';
          PASREINI = FAUX;
*
*         ----------------------------------------------------------------
*         Recalcul de la rigidite a la fin du pas 
          'FORM' GEOM2 ;
          txt_k ='CHAI' ' Recalcul de K (= K^elast';
          'SOUC' 0;
          'SI' ('EGA' LAG_TOT 1) ;
            'FORM' GEOREF0;
            HOOKRH  = 'HOOK' ZMODL ZMAT2I ;
            'FORM' GEOM2;
            HOOKRH2 = 'CONF' HOOKRH ZMODL ;
            RITC = 'RIGI' HOOKRH2 ZMODL ZMAT22 'NOER';
            'DETR' HOOKRH; 'DETR' HOOKRH2;
          'SINON';
            RITC = 'RIGI' ZMODL ZMAT22 'NOER' ;
          'FINSI';
          'FORM' GEOM1;
          'SI' (SOUCI) ;
            'MESS' 'rigi rate. on reessaye avec la configuration geom1';
            RITC = 'RIGI' ZMODL ZMAT21 ;
          'FINSI';
*
*         RITC peut contenir des CL
          ZRI = 'EXTR' RITC 'RIGI' 'NOMU' ;
          ZCL = 'EXTR' RITC 'RIGI' 'MULT' ;
          'SI' ('NEG' ('DIME' ZCL) 0) ;
            ZCLIM = ZCLIM0 'ET' ZCL ;
          'SINON' ;
            ZCLIM = ZCLIM0 ;
          'FINSI' ;
*
*         ----------------------------------------------------------------
*         Option 'AUGMENTATION_AUTOMATIQUE' : actualiser la rig augm
          'SI' (AUTAUG 'ET' ('NON' IRAUGLU));
            'FORM' GEOM2 ;
            RIG_AUG = 'MASSE' ZMODL ZMAT22 ;
            RH      = RITC;
            'FORM' GEOM1;
***         mess 'actualisation rig_aug';
          'FINSI';
*
*         ----------------------------------------------------------------
*         Ajout de relations issues de la procedure DEFO_IMP
          ZFCONT   = 0 ;
          IACTURES = FAUX ;
          ZDEPTR   = ZDEP0 '+' ZDEPT ;
          CDEPSLX  = 'ENLE' ZDEPTR 'LX' ;
 
**        DE_CNTRL = FAUX;
          'SI' (AUTAUG 'ET' DE_CNTRL);
            IACTURES = VRAI ;
            ZFCONT = RE_CNTRL '*' CDEPSLX 'EXCO' 'FLX' 'NOID' 'FLX'  ;
            ZCLIM  = RE_CNTRL 'ET' ZCLIM ;
            txt_k  = 'CHAI' txt_k ' + K^control';
          'FINSI';
*
*         ----------------------------------------------------------------
*         Recalcul des conditions de contact-frottement ds la conf finale
          'SI' WTAB.'CONTACT';
            txt_k  = 'CHAI' txt_k ' + K^contact';
            'FORM' GEOM2 ;
*
*     
            MODCON  = wtab.'MODCONTA';
            'SI' ('EXIS' WTAB 'MATCONTA');
              CJEU CRR = 'RFCO' MODCON WTAB.'CONV' WTAB.'MATCONTA';
            'SINON' ;
              CJEU CRR = 'RFCO' MODCON WTAB.'CONV'                ;
            'FINS' ; 
            'SI' WTAB.'CAFROTTE';
              RTRSF*'RIGIDITE' RFNS*'RIGIDITE'  FFROT = 'EXCF' CRR ZDEPT MODCON WTAB.'MATCONTA'  excfconv;
              MFROT = 'EXTR' FFROT 'MAIL' ;
            'FINSI';
            ZFLX = 'EXCO' (CRR '*' ZDEPT) 'FLX' 'NOID' 'FLX' ;
*
*           Separer ce qui concerne l'adherence (FADH) et les jeux (FLX)
            'SI' ('NEG' CJEU 0) ;
              FADHE = 'EXCO' CJEU 'FADH' 'NOID' 'FADH';
              CDAP  = 'EXCO' CJEU 'FLX'  'NOID' 'FLX' ;
*
*             Ajout du glissement deja parcouru sur les conditions de frot.
              'SI' WTAB.'CAFROTTE';
                'SI' ('EGA' ('TYPE' MFROT) 'MAILLAGE');
                  CGLI = 'REDU' ZFLX MFROT ;
                  CDAP = CDAP ET (-1*CGLI) ;
                'FINSI'; 
              'FINSI'; 
            'FINSI'; 
*
            'SI' WTAB.'MODAL' ;
              CRR  = 'PJBA' CRR  ZMODL ZMAT22 ;
              MCRR = 'EXTR' CRR 'MAIL' 'MULT';
              MCDAP = 'EXTR' CDAP 'MAIL' ;
 
              'REPETER' BCDA ('NBNO' MCDAP) ;
                 PBCDA = MCDAP 'POINT' &BCDA ;
                 PCRR  ='POINT' MCRR &BCDA   ;
                 CHCR  ='MANU' 'CHPO' PCRR 1 'FLX' ('EXTR' CDAP 'FLX' PBCDA)
                        'NATURE' 'DISCRETE' ;
                'SI' ('EGA' 1 &BCDA) ;
                   CCDA = CHCR ;
                'SINON' ;
                   CCDA = CCDA 'ET' CHCR ;
                'FINS' ;
              'FIN' BCDA ;
 
              CDAP = CCDA ;
            'FINS' ;
*
            'SI' ('NEG' CRR 0) ;
              IACTURES = VRAI ;
              CHPZ     = CRR '*' CDEPSLX 'EXCO' 'FLX' 'NOID' 'FLX'  ;
              'SI' ('NEG' ZFCONT 0) ;
              ZFCONT = ZFCONT '+' CHPZ ;
              'SINON' ;
              ZFCONT = CHPZ ;
              'FINSI' ;
              ZCLIM  = CRR 'ET' ZCLIM ;
            'FINS';
*
            'SI' ('NEG' CJEU 0) ;
              IACTURES = VRAI ;
              ZFCONT = ZFCONT '+' CDAP ;
            'FINS';
            'FORM' GEOM1 ;
          'FINS';
*
*         Mise a jour des jeux dans le residu (partie force inchangee)
          ZFCONSTA = ZFEXT2 ;
          ZFLX1    = ZFLXB ;
          'SI' IACTURES ;
            ZFCONSTA = ZFCONSTA '+' ZFCONT ;
            RFORCE   = 'ENLE' RESIDU 'FLX' ;
            FLXPRE   = 'EXCO' (ZCLIM '*' ZDEPTR) 'FLX' 'NOID' 'FLX' 'NATURE' 'DISCRET' ;
            RJEUX    = ZFCONT '+' ZFLXB '-' FLXPRE ;
            RESIDU   = RFORCE 'ET' RJEUX ;
            ZFLX1    = ZFLX1  'ET' ZFCONT ;
          'FINS';
*
          'SI' WTAB.'CONTACT' ;
            XXX1   = 'EXTR' CRR 'MAIL' 'UNIL';
            XXX2   = 'DIFF' ('EXTR' ZCLIM 'MAIL' 'MULT') XXX1 ;
            FLXREC = 'REDU' FLXINI XXX2 ;
            FLXREC = FLXREC 'ET' ('REDU' FLXPRE XXX1) ;
            DUUNIL = CDAP ;
          'FINS';
*
          'DETR' ZRAID;
* A-t-on des C.L. unilaterales ?
          MCLIMU = 'EXTR' ZCLIM 'MAIL' 'UNIL';
          IMPO12 = ('NBEL' MCLIMU) '>' 0 ;
          'SI' IMPO12 ;
            DUUNIL = 'REDU' DFEXT0L MCLIMU ;
            'SI' ( WTAB.'CONTACT' et ('NEG' CJEU 0)) ;
              DUUNIL = CDAP 'ET' DUUNIL ;
            'FINSI' ;
            XXX2   = 'DIFF' ('EXTR' ZCLIM 'MAIL' 'MULT') MCLIMU ;
            DUIMPO = 'REDU' DFEXT0L XXX2 ;
          'FINS';
          ZRAID = ZCLIM 'ET' ZRI;
*         Mise a jour de FREAP (faut-il faire cette actualisation?)
          FREAP = (ZDEP0 'EXCO' 'LX' 'LX')  * ZCLIM;
*
*         ----------------------------------------------------------------
*         Ajout de la partie dynamique
          'SI' IDYN;
*           bp : en toute rigueur, il faudrait aussi recalculer la MASSE ...
*                et ajouter l'amortissement le cas échéant ...
            ZRAID = ZRAID 'ET' ZRMAS ;
          'FINS';
*
*         ----------------------------------------------------------------
*         Recalcul de la raideur geometrique (option K_SIGMA)
          'SI' IKSIA ;
            'FORM' GEOM2;
            ZSIGKS = 'CONF' ZSIGF ZMODL ;
            KSIGTC = 'KSIGMA' ZMODL ZSIGKS ZMAT22;
            ZRAID  = ZRAID 'ET' ksigtc ;
            txt_k  = 'CHAI' txt_k ' + K^sigma';
            'FORM' GEOM1;
          'FINS';
*
*         ----------------------------------------------------------------
*         Prise en compte d'eventuelles 'RIGIDITE_CONSTANTE'
          'SI' IRCON;
            ZRAID = ZRAID 'ET' RIG_CONS;
            txt_k = 'CHAI' txt_k ' + K^constant';
          'FINS';
*
*         ----------------------------------------------------------------
*         Traitement en de la rigidite augmentee
          'SI' IRAUG;
            'SI' AUTAUG ;
              xkx = (xtmx zdep1d rh)  + xpetit;
              xmx = (xtmx zdep1d rig_aug) + xpetit;
              'SI' (('VERI' xkx) 'ET' ('VERI' xmx) 'NON');
                zu1l = ZDEPT 'ENLE' 'LX';
                xkx = (xtmx zu1l rh)  + xpetit;
                xmx = (xtmx zu1l rig_aug) + xpetit;
              'FINSI';
              augm = xkx / xmx; augk = 1.d-1 ; augm = augm * 1d1;
**              augk    = augk    * (1. + (wtab.'ISOUSPAS' * 2.) );
              augm = abs augm;
              augk    = abs augk;
              augmult = augmult * 1.02;
              'SI' (augmult > 1D3 ); augmult=1D3 ; 'FINSI';
 
              ZRAID = ZRAID 'ET' (RIG_AUG * (augm * augmult)) 'ET' (RH * (augk * augmult)) ;
 
              'MESS' 'multiplicateur d augmentation masse' (augm * augmult) ' raideur' (augk * augmult) ;
            'SINON';
              ZRAID = ZRAID 'ET' RIG_AUG ;
            'FINSI';
            txt_k = 'CHAI' txt_k ' + K^augmente';
***         ZNACCE = FAUX;
          'SINON';
            augmult = augmult * 0.55 ;
            si (augmult < 1d-3); augmult = 1d-3; finsi;
          'FINSI';
          'FORM' GEOM1;
*  
*         ----------------------------------------------------------------
*         Procedure UTILISATEUR
          'SI' WTAB.'PROCEDURE_CHARMECA';
            PRECED.'ADDI_MATRICE' = vrai;
            TFP22= CHARMECA PRECED TFIN ;
            PRECED.'ADDI_MATRICE' = faux;
            'SI' (EXIS TFP22 'ADDI_MATRICE');
              zraid = zraid 'ET' TFP22.'ADDI_MATRICE';
              txt_k = 'CHAI' txt_k ' + K^cent';
            'FINS';
          'FINS';
*
          'MESS' ('CHAI' txt_k ' ) dans config deformee ') ;
*
*         on impose le recalcul de K a la prochaine iteration si it=2 (pq ?)
          'SI' (IT 'EGA' 2) ;
            ITACC = 4 ;
          'FINS' ;
        'FINS';
      'FINS';
*----------------------------------------------------------------------*
*            Evaluation de la matrice tangente si demandee             *
*----------------------------------------------------------------------*
      'SI' IKTAN ;
        'SI' ('NON' IFEFP) ;
          'SI' (('MULT' IT WTAB.'NITER_KTANGENT') 'ET' (IT > 1)) ;
*
            IKT = VRAI ;
            'SI' (IGRD 'ET' ('NON' HPP_EPS));
              'FORM' GEOM2 ;
            'FINS' ;
 
            'SI' ('NON' ISSTE) ;
              'SI' IPERT ;
                'SI' PARTLOCA ;
                  zktap = 'ASSIS' 'TOUS' 'KTAN' 'PERT' MODRELOC Z1COMP Z2COMP
                           'C1' WTAB.'K_TANG_PERT_C1' 
                           'C2' WTAB.'K_TANG_PERT_C2'
                           ZKTASYM ;
                  ZRIKTA = 'ET' zktap ;
                'SINON';
                  ZRIKTA = 'KTAN' 'PERT' ZMODL Z1COMP Z2COMP
                           'C1' WTAB.'K_TANG_PERT_C1'
                           'C2' WTAB.'K_TANG_PERT_C2'
                           ZKTASYM ;
                'FINS';
              'SINON' ;
                'SI' (IVISCO 'OU' IVIDOM) ;
                  DTTAN = ZDT ;
                'FINS' ;
                ZRIKTA = 'KTAN' ZMODL ZSIGF ZVARF ZMAT22
                                'PREC' ZPREK 'DT  ' DTTAN ZKTASYM ;
              'FINS' ;
            'FINS' ;
 
            'SI' IKSIA ;
              ZKSIG  = 'KSIGMA' ZMODL ZSIGF ZMAT22 ;
              ZRIKTA = ZRIKTA 'ET' ZKSIG ;
            'FINS' ;
 
            ZRAID = ZRIKTA 'ET' ZCLIM ;
 
            'SI' IGRD ; 'FORM' GEOM1 ; 'FINS' ;
 
          'FINS' ;
        'FINS' ;
      'FINS' ;
    'FINSI' ;
*FOR_MECA
*
    'SI' WTAB.'NVSTNL'                          ;
      UN = 'EXCO' ZVIFL movifl                  ;
      ch_un = 'MANU' 'CHML' ZMODL 'ADVE' un 'TYPE' 'CARACTERISTIQUES';
      nugrad = 'ADVE' ZMODL ch_un      ;
      'SI' ('EXIS' PRECED 'OPTION_BDF2') ;
        ZRNSL = WTAB.'RNSL' 'ET' nugrad 'ET' (WTAB.'NSMA' '*' 1.5)  ;
      'SINON'                                   ;
        ZRNSL = WTAB.'RNSL' 'ET' nugrad 'ET' WTAB.'NSMA' ;
      'FINSI' ;
      ftnsl  = ZFCONSTA 'ET' fmass              ;
 
      'SI' LOGADV                               ;
        ch_ad = 'MANU' 'CHML' ZMODL 'ADVE' UADV 'TYPE' 'CARACTERISTIQUES';
        NUADV = 'ADVE' ZMODL ch_ad     ;
        ZRNSL = ZRNSL 'ET' NUADV                ;
        ftnsl = ftnsl  '-' fadv                 ;
      'FINSI'                                   ;
 
      resnsl = ftnsl '-' (ZRNSL * ZVIFL)        ;
 
    'FINSI'                                       ;
*----------------------------------------------------------------------*
*                Acceleration de convergence effective                 *
*----------------------------------------------------------------------*
    'SI' WTAB.'FOR_MECA' ;
      CORRECA = CORREC ;
      CORREC = 0;
      'SI' (resmul < 0.99);
        'MESS' 'Reduction du chargement. Coefficient: ' resmul;
        residu = ((residu -fcorf) * resmul) +  fcorf;
      'FINSI';
      ACFEP0 = RESIDU - FCORF;
*      'SI' ('EGA' &etiq 1); ACFEP0 = 'VIDE' 'CHPOINT'; 'FINSI';
      ACFEP0 = ACFEP0 'ENLE' 'FLX';
      ACFEPT = ACFEP0;
      ACFP0  = ACFEP0 * ZJAC ;
      ACFEP0 = ACFEP0 - CORRECA;
      'SI' ((ITACC '&lt;EG' 0) 'ET' (IT '>' 3) 'ET' ZNACCE);
        ITACC = 2;
 
 
        CORREC = 'ACT3' ACFEP2 ACFEP1 ACFEP0
                        ACFP3  ACFP2  ACFP1  ACFP0 ;
*
*       verif que l'acceleration ne renvoie pas en arriere
        'SI' (WTAB.'STABILITE' 'OU' IPILOT) ;
          acc_ref = xty acfept zdep1d MNDUAL MNPRIM;
          acc_ref = acc_ref + xpetit;
 
          acc_dir = xty (acfept - correc) zdep1d MNDUAL MNPRIM;
          acc_rap = acc_dir/acc_ref;
          acc_lim = 1E5  ; 'SI' (acc_ref < 0.); acc_lim = 1.; 'FINSI';
          'SI' (acc_rap  '>' acc_lim) ;
*       en cas d'acceleration trop grande, on la limite
            'MESS' 'Limitation acceleration ' ' 'acc_rap;  CORREC = correc * (acc_lim/acc_rap);
          'FINSI';
          'SI' (acc_rap '&lt;EG' 0.) ;
*       en cas d'acceleration retrograde, on n'accelere pas
            'MESS' 'Annulation acceleration: retrograde' ' ' acc_rap;
            correc =  0.;
          'FINSI';
        'FINSI';
*
        RESIDU = RESIDU '-' CORREC;
      'FINS';
*
      ACFP3  = ACFP2  ;
      ACFP2  = ACFP1  ;
      ACFP1  = ACFP0  ;
      ACFEP2 = ACFEP1 ;
      ACFEP1 = ACFEP0 ;
*
*----------------------------------------------------------------------*
*                 Systeme a resoudre a l'iteration (i)                 *
*                                                                      *
*           | K  At | |  du   |   | Fsur + Fvol - Bsigma^(i) |         *
*           |       | |       | = |                          |         *
*           | A  0  | | LX(i) |   | d1 - A*u^(i)             |         *
*----------------------------------------------------------------------*
      'SOUCI' 0;
      FEXCI = 'VIDE' 'CHPOINT' / 'DISCRET' ;
      'SI' WTAB.'ADHERENCE' ;
        'SI' ('EXIS' WTAB 'MAIL_BLOM') ;
          FADRE = 'REDU' FADHE WTAB.'MAIL_BLOM' ;
        'SINON' ;
          FADRE = FADHE ;
        'FINSI' ;
        FEXCI = FEXCI 'ET' FADRE ;
      'FINSI' ;
*
*     Force limite de frottement
      'SI' (WTAB.'CAFROTTE' 'ET' IMPO12);
        'SI' WTAB.'FROCABL' ;
          zsigfT = 'REDU' zsigf ZMODLI;
          FFROT = 'EXCC' ZCLIM ZDEPT MODCON WTAB.'MATCONTA' ZSIGFT;
          RTRSF = 'VIDE' 'RIGIDITE' / 'RIGIDITE';
          RFNS  = 'VIDE' 'RIGIDITE' / 'RIGIDITE';
          FEXCI = FEXCI '+' FFROT ;
        'FINSI';
        'SI' WTAB.'FROCOUL' ;
          excfconv = excfconv + 1 ;  
          RTRSF*'RIGIDITE' RFNS*'RIGIDITE'  FFROT = 'EXCF' ZCLIM ZDEPT MODCON WTAB.'MATCONTA' excfconv;
          FEXCI = FEXCI '+' FFROT ;
          MFROT = 'EXTR' FFROT 'MAIL' ;
        'FINSI';
      'FINS';
*
      RESIDU = 'CHAN' 'ATTRIBUT' RESIDU 'NATURE' 'DISCRET' ;
      'SI' ( ('NEG' ZRIBLO_M 'INCONNU') 'ET' (IPREM 'OU' RECA_K) ) ;
        'SI' ((WTAB.'CAFROTTE' 'OU' WTAB.'ADHERENCE') 'ET' IMPO12);
          ZDEP1 BID BID BID BID     = 'RESO' ZRAID 'SOUC' RESIDU 'INIB'
                                   ZRIBLO_M ZLISEA_M
                                   FEXCI RFNS RTRSF;
        'SINON';
          ZDEP1 BID BID BID     = 'RESO' ZRAID 'SOUC' RESIDU 'INIB'
                                   ZRIBLO_M ZLISEA_M ;
        'FINS';
      'SINON';
        'SI' ((WTAB.'CAFROTTE' 'OU' WTAB.'ADHERENCE') 'ET' IMPO12);
          ZDEP1 = 'RESO' ZRAID 'SOUC' RESIDU FEXCI       RFNS RTRSF;
        'SINON';
          ZDEP1 = 'RESO' ZRAID 'SOUC' RESIDU 'NOID';
          'SI' (WTAB.'MAN' 'ET' IPREM);
            ORDRE      = WTAB.'ORDRE' ;
            ZDEP2 IOUT = CORMAN ZRAIDINI ZMODL ZMAT22 ORDRE
                          ZDEP0 ZSIG0      RESIDNOR  WTAB ;
            'SI' (IOUT 'EGA' 1) ;
              ZDEP1=ZDEP2;
            'FINS';
          'FINSI' ;
        'FINS';
      'FINS';
*
*     verif stabilite matrice elastique
      'SI' (wtab.'STABILITE' 'ET' ('NON' HPP_EPS) 'ET' AUTAUG);
        RESIDD = RESIDU - (ZRAID * (ZDEP1 'EXCO' 'LX' 'LX'));
        acc_ref= xty (RESIDD       ) zdep1 MNDUAL MNPRIM;
** diag refactorise la matrise. On evite donc provisoirement de l'appeler
**      nbng = 'DIAG' ZRAID;
        nbng = 0;
        'SI' ((acc_ref <   0.) 'OU' (nbng > 0));
          'MESS' 'Raideur negative' ' ' acc_ref ;
**        resmul = resmul * 0.25;
**          augmult = augmult * 1.4;
          'SI' LOG_CNTRL;
*  controler la deformation max
            DEPST = 'EPSI' 'LINE' ZMODL (ZDEPT + ZDEP1) ZMAT22;
            DPSMAX = 'MAXI' 'ABS' DEPST 'AVEC' MLDEFOR;
*            PRIMAX = DEPST 'PRIN' ZMODL ZVAR0;
*            DPSMAX = 'MAXI' 'ABS' PRIMAX 'AVEC' MLPRINC;
*            'DETR' PRIMAX;
            RE_CNTRL2*'RIGIDITE'   = DEFO_IMP DEPST  (zdept + zdep1)  ZRAID ;
            DE_CNTRL = VRAI;
            ZFCONT = RE_CNTRL2 '*' (        ZDEP1) 'EXCO'  'FLX' 'NOID' 'FLX'  ;
**          cmul = MINI  1. ( 1d-4 / DPSMAX);
**          ZFEXT2   = ZFEXT2   + (ZFCONT * cmul);
            ZFLXB    = 'EXCO' ZFEXT2 'FLX' 'NOID' 'FLX' ;
            ZFCONSTA = ZFEXT2 ;
            'SI' (RED_URG < 3);  RE_CNTRL = RE_CNTRL 'ET' RE_CNTRL2; 'FINSI';
          'FINSI';
**        'ITERER' ETIQ;
** essai en inversant zdep1 pour avoir la positivite du travail
          zdep1 = ((zdep1 'ENLE' 'LX') * -1.) 'ET' (zdep1 'EXCO' 'LX' 'LX');
        'FINSI';
      'FINSI';
*
      'SI' ('SOUCI') ; 
        'MESS' 'Souci dans la resolution :'('VALE' 'SOUC');
        'SI' AUTAUG ;
          resmul  = resmul * 0.25;
          augmult = augmult * 1.4;
        'FINSI';
        pastest = VRAI;
        HPP_EPS = FAUX;
        URG     = VRAI;
        ITACC   = 4;
***     'ITERER' ETIQ;
      'FINSI';
*
      'SI' LOGPIL ; 
*       Calcul de D_eta par appel a la procedure PILOINDI
        D_ETA = PILOINDI PRECED ZDEP0 ZDEPT ZDEP1 ZDEPII DTAU ;       
*       Mise à jour  
        ZDEP1 = ZDEP1 '+' (D_ETA '*' ZDEPII) ;
        ETA = ETA '+' D_ETA ;
      'FINSI' ;
*
      'SI' (resmul < 1d-2) ; souci 0; 'FINSI';
      pasunil = 'SOUCI' ;
*
      ZRAID_T = 'EXTR' ZRAID 'CONT';
      'SI' ('EXIS' ZRAID_T 'NITER');
*       Conserver les conditions unilaterales pour VITEUNIL
        'SI' IDYN ; WTAB.'ZRAIDV' = ZRAID; 'FINS';
        ZRIBLO_M = ZRAID_T.7;
        ZLISEA_M = ZRAID_T.6;
        WTAB.'MAIL_BLOM' = 'EXTR' ZRIBLO_M 'MAIL' 'MULT' ;
      'FINS';
    'FINSI' ;
*FOR_MECA
 
    'SI' WTAB.'NVSTNL'             ;
      ZVIF1  = 'KRES' ZRNSL Resnsl ;  
      ZVIFL = ZVIFL '+' ZVIF1      ;
    'FINSI'                        ;
*----------------------------------------------------------------------*
*            Increment de deplacements et mult. de Lagrange            *
*                                                                      *
* ZDEP1 est un increment d'increment de deplacements, note du, et les  *
* multiplicateurs de Lagrange sont resolus "en total".                 *
* La solution de la boucle ETIQ a l'iteration i a fourni :             *
*                ZDEP1  = [ du  ; LX(i) ]                              *
* Le cumul est realise de maniere a avoir :                            *
*                ZDEPT  = [ Du^(i) ; LX(i) ]                           *
* avec Du^(i) = Du^(i-1) + du                                          *
* Le champ de deplacements total s'ecrit :                             *
*                ZDETOT = ZDEP0 '+' ZDEPT                              *
*----------------------------------------------------------------------*
    'SI' WTAB.'FOR_MECA' ;
      'SI' IPREM ;
        ZDEPT = 'COPIER' ZDEP1 ;
        ZDELA = 'COPIER' ZDEPT ;
      'SINON';
        XXX1  = ZDEPT 'ENLEVER' 'LX' ;
        ZDEPT = XXX1 '+' ZDEP1 ;
        'DETR' XXX1 ;
      'FINS' ;
*
      'SI' WTAB.'MODAL' ;
        'SI' ('EXIS' WTAB 'MODCONTA') ;
          'SI' IPREM ;
            mamoco1 = 'EXTR' (ZDEPT 'ENLEVER' 'LX') 'MAIL' ;
          'SINON' ;
*           mettre les point materiels dans zdept
            zdeptu1 =  'REDU' zdept   mamoco1      ;
            ch_dco  =  'RECO' zdeptu1 ZMODL ZMAT22 ;
            ZDEPT   = ('EXCO' zdept 'LX' 'LX') 'ET' zdeptu1 'ET' ch_dco ;
          'FINS';
        'FINS';
      'FINS' ;
*
*     Option automatique : determiner le coef de normalisation
      'SI' IPILOT;
        PASTEST = '<'    IT WTAB.'AUTORECA' ;
        ACCEL   = 'MULT' IT WTAB.'AUTORECA' ;
        'SI' ACCEL ;
*
*         MODIFICATION CB215821 : 18/06/2015
          'SI' (RED2 '<' (IT '/' 20)) ;
            RED2 = RED2 '+' 1 ;
            ZAUTOREDU = ZAUTOREDU '*' 3.D0 ;
          'FINS';
          'SI' (ZAUTOREDU '>' 1.D0 ) ;
            'MESS' 'On divise le critere de pilotage par 'ZAUTOREDU;
          'FINS' ;
*
          OO   = WTAB.'AUTOCRIT' '/' ZAUTOREDU ;
          U1MA = AUTOPILO ZDEPT (COEPI'*'ZDELA) ZMODLI ZMAT2I WTAB;
          AL1  = OO '/' U1MA ;
*
*         AL1 : coefficient de normalisation
          'SI'((AL1 '>' 1.D0) 'ET' (COEPI '>' 0.D0));
            'SI' (AL1 '>' (1.D0 '/' COEPI));
              AL1 = 1.D0 '/' COEPI;
            'FINS';
          'FINS';
*
*         Normalisation
          XXX1  = (1.D0 '-' AL1) '*' ZLX ;
          XXX3  = AL1 '*' ZDEPT ;
          ZDEPT = XXX3  '+' XXX1;
          'DETR' XXX3;
        'FINSI' ;
*
      'SINON' ;
*       On part dans le decor : redemarrer a 0
        'SI' ((XCONV > 1E8) 'ET' PASREINI) ;
          'MESS' 'Reinitialisation du schema';
          PASREINI=FAUX;
          'ITERER' ETIQ;
        'FINS';
      'FINS';
*
*     garder les reactions pour le test de convergence
      ZDEPLP = ZDEPL                         ;
      ZDEPL  = ZDEPT 'EXCO' 'LX' 'NOID' 'LX' ;
*----------------------------------------------------------------------
*     le nouveau champ est fixe on va tester l'equilibre(convergence)
*     et calculer la force motrice pour l'iteration suivante
*----------------------------------------------------------------------
      'DETRUIRE' FCORF;
*
*     Champ de deplacement total
      ZDETOT = ZD0SLX '+' ZDEPT ;
*
      'SI' IGRD;
*
*       Configuration fin de pas
        GEOM2  = 'FORM' GEOREF0 ZDETOT ;
        ZMAT22 = 'CONF' ZMODL ZMAT21 ;
        'FORM' GEOM1 ;
*
*       Gradient des deplacements (total) a la fin du pas de temps 
*       par rapport a la configuration de reference GEOREF0
        'SI' ('NEG' ZGRDU0 'INCONNU');
          'FORM' GEOREF0 ;
          ZGRDUF = 'GRAD' ZMODL ZMAT2I ZDETOT ;
          D_GR_U = ZGRDUF '-' ZGRDU0 ;
          'FORM' GEOM1 ;
        'FINS';
*
      'FINS';
*
*     XXX1 = [ FCORF ; FCORU ] = [ -Freac_it ; u^imp_it ]
      XXX1  = ZCLIM '*' ZDETOT;
      FCORF = 'ENLE' XXX1 'FLX';
      FCORU = 'EXCO' XXX1 'FLX' 'NOID' 'FLX';
      'DETR' XXX1 ;
*----------------------------------------------------------------------*
*              Increment de deformations totales DEPST                 *
*----------------------------------------------------------------------*
*     Update or total lagrangian ---------------------------------------
      'SI' IFEFP;
        GEOM2  = 'FORM' GEOREF0 ZDETOT ;
        'SI' IFEFPUL ;
*         mess ' update lagrangian ZRIKTA';
***       GEOM2 = 'FORM' ZDEPT ;
          ZRIKTA ZSIGF ZVARF ZDEIF = 'ECFEFP'
                ZMODL ZDEI0 ZVAR0 ZDEPT ZMAT22 ZPREK NITMA 1 ;
        'SINON';
*         mess ' total lagrangian ZRIKTA';
***       GEOM2 = 'FORM' ZDETOT ;
          ZRIKTA ZSIGF ZVARF ZDEIF = 'ECFEFP'
                ZMODL ZDEI0 ZVAR0 ZDETOT ZMAT22 ZPREK NITMA ;
        'FINS';
        FEQU2 = 'BSIG' ZMODL ZSIGF ZMAT22 ;
        ZRAID = ZRIKTA 'ET' ZCLIM ;
        'FORM' GEOM1 ;
*
        XXX1 = 'EXCO' ZVARF WTAB.'MOVA' ;
        EPSM = 'MAXI' 'ABS' XXX1 'AVEC' MLDEFOR ;
**        EPSPR= XXX1 'PRIN' ZMODL ZVAR0;
**        EPSM = 'MAXI' 'ABS' EPSPR 'AVEC' MLPRINC;
**        'DETR' EPSPR;
        ACC  = 'ABS' ( XXX1 - ACC0 ) ;
        MMC  = 'MASQUE' ACC 'SUPERIEUR' 1.D-10 'SOMME' ;
        'SI' (MMC '>' MMCMAX) ; MMCMAX = MMC ; 'FINS' ;
        DPSMAX = 'MAXI' 'ABS' ACC 'AVEC'  MLDEFOR ; 
**        DPSPRO = ACC 'PRIN' ZMODL ZVAR0;
**        DPSMAX = 'MAXI' 'ABS' DPSPRO 'AVEC' MLPRINC;
**        'DETR' DPSPRO;
        DEPST  = 'CHAN' 'TYPE' ZDEIF 'DEFORMATIONS' ;
        ZDEFF  = ZDEF0 '+' DEPST ;
*
*     cas standard -----------------------------------------------------
      'SINON';
*
        ZDEF0a = ZDEF0 ;
        ZSIG0a = ZSIG0 ; 
        ZVAR0a = ZVAR0 ;
        ZDEI0a = ZDEI0 ;
        DEPST  = 0. ;
        DSIGT  = 0. ;
        ZMATT  = ZMAT22 ;
*
        nsoincr = nsoincrn ;
        ZDEPTI  = ZDEPT '/' nsoincr ;
        ZDEPTI2 = ZDEPTi / 2. ;
        ZDEPTI2n = ZDEPTi2 * -1. ;
*
*       ---------------------------------------------------------------
*       Sous-incrementation du comportement
        'REPETER' sousinc nsoincr;
*
          codeb = 'FLOT' (&sousinc '-' 1) '/' nsoincr ;
          cofin = 'FLOT'  &sousinc        '/' nsoincr ;
*
          'SI' IGRD;
            ZDEPF = ZDEPT '*' cofin ;
            ZDEPM = ZDEPT '*' ((&sousinc - 0.5) '/' nsoincr) ;
            GEOMF = 'FORM' GEOM1 ZDEPF  ;
            GEOMM = 'FORM' GEOM1 ZDEPM  ;
            ZMATT = 'CONF' ZMODL ZMAT21 ;
            'FORM' GEOM1 ;
          'FINS';
*
*         --------------------------------------------------------------
*         Calcul de l'increment de deformations DEPST
*
*         Pour les hypotheses de deformations non lineaires,
*         on calcule l'increment de deformation sur le pas avec EPSI 'LINE'
*         mais en se placant sur une configuration intermediaire (mi pas)
          iDEFOK = FAUX ;
          'SI' (('NEG' HYPDEF 'LINEAIRE') 'ET' IGRD);
*
            'FORM' GEOMM ;
            'SOUC' 0;
            'SI' ('OU' IENDOM IVIDOM ICERAM);
              HOOKENDO = 'HOOK' ZMODL ZMAT22 ZVAR0 ;
              DEPSTa   = 'EPSI' 'LINE' 'NOER' ZMODL ZDEPTi HOOKENDO ZMAT22;
              'DETR' HOOKENDO;
            'SINON';
**              DEPSTa   = 'EPSI' 'LINE' 'NOER' ZMODL ZDEPTi ZMATT ;
 
              DEPSTa1  ='EPSI' 'QUAD' 'NOER' ZMODL ZDEPTi2 ZMATT;
              DEPSTA2 = 'EPSI' 'QUAD' 'NOER' ZMODL ZDEPTi2n ZMATT;
              DEPSTa   = DEPSTA1 - DEPSTA2;
              'DETR' depsta1; 'DETR' depsta2;
            'FINSI';
            iDEFOK = 'NON' ('SOUCI') ;
            DPSMAX=1001 ;
            'SI' iDEFOK ;
              DPSMAX = 'MAXI' 'ABS' DEPSTa 'AVEC' MLDEFOR;
**              DPSPRI = DEPSTa 'PRIN' ZMODL ZVAR0;
**              DPSMAX = 'MAXI' 'ABS' DPSPRI 'AVEC' MLPRINC;
**              'DETR' DPSPRI;
 
            'FINSI' ;
            'FORM' GEOM1; 
*
            'SI' (('NON' iDEFOK) 'OU' (DPSMAX > 1000))  ;
              MSGZ = 'CHAI' 'Deformation non lineaire non calculable' ;
              'SI' autaug;  
                mess ('CHAI' MSGZ ' : annulation sous-increment'); 
                augmult = augmult * 1.4;
                resmul = resmul * 0.25;
                'ITERER' ETIQ;
              'FINSI';
              MESS ('CHAI' MSGZ ' : deformation lineaire'); 
              pastest = vrai;
            'SINON';
*             Increment de deformations transporte sur GEOREF0
              'SI' IGRD ;
                'SI' IJAUMA ;
                  DEPSTa = 'CONF' DEPSTa ZMODL ;
                  'FORM' GEOREF0 ;
                  DEPSTa = 'CONF' DEPSTa ZMODL ;
                'SINON' ;
                  'FORM' GEOREF0;
                  DEPSTa = 'CONF' DEPSTa ZMODL ;
*                 'FORM' GEOM1 ;
                'FINSI' ;
              'FINSI';
              iDEFOK = VRAI ;
            'FINS';
          'FINS';
*
*         Calcul des deformations lineaires (demande ou pas reussi)
          'SI' ('NON' iDEFOK) ;
            'SI' ('OU' IENDOM IVIDOM  ICERAM);
              HOOKENDO = 'HOOK'        ZMODL ZMAT22 ZVAR0         ;
              DEPSTa   = 'EPSI' 'LINE' ZMODL ZDEPTi HOOKENDO ZMAT22;
              'DETR' HOOKENDO;
            'SINON';
              DEPSTa = 'EPSI' 'LINE' ZMODL ZDEPTi ZMAT22;
            'FINSI';
            'SI' (('NEG' HYPDEF 'LINEAIRE') 'ET' IGRD);
              'MESS' 'Attention utilisation des deformations lineaires';
            'FINS';
*
*           Increment de deformations transporte sur GEOREF0
            'SI' IGRD ;
              'FORM' GEOREF0;
              DEPSTa = 'CONF' DEPSTa ZMODL ;
            'FINSI';
          'FINS';
*
          ZDEFF  = ZDEF0a '+' DEPSTa ;
          DEPSTZ = DEPST  '+' DEPSTa ;
          'DETR' DEPST ;
          DEPST  = DEPSTZ ;
          DPSMAX = 'MAXI' 'ABS' DEPST 'AVEC'  MLDEFOR;
**          DPSPRI = DEPST 'PRIN' ZMODL ZVAR0;
**          DPSMAX = 'MAXI' 'ABS' DPSPRI 'AVEC' MLPRINC;
**          'DETR' DPSPRI;
*
*         --------------------------------------------------------------
*         Transport des quantites sur la configuration adequate
*
          'SI' IGRD ;
*           Lagrangien total : sur GEOREF0
            'SI' ('EGA' LAG_TOT 1);
              'FORM' GEOREF0;
            'FINSI';
*           Lagrangien fin pas : GEOREF0 -> GEOMF
            'SI' ('EGA' LAG_TOT 2);
              'FORM' GEOMF ;
            'FINS';
*           Lagrangien mi pas : GEOREF0 -> GEOMM
            'SI' ('EGA' LAG_TOT 3);
              'FORM' GEOMM;
            'FINS';
            ZSIG0Z ZDEF0Z DEPSTZ ZMATT = 'CONF' ZMODL ZSIG0a ZDEF0a DEPSTa ZMAT21 ;
            'SI' ('NEG' &sousinc 1) ;
              'DETR' ZSIG0a ; 'DETR' ZDEF0a ; 'DETR' DEPSTa ;
            'FINSI' ;
            ZSIG0a ZDEF0a DEPSTa = ZSIG0Z ZDEF0Z DEPSTZ ;
          'FINS';
*
          'SI' IDYN ;
            'SI' (('NEG' nsoincr 1) 'ET' ('EGA' &sousinc 1));
              ZSIGM0 = ZSIG0a ;
              ZSIGM  = ZSIGM0 '/' 2;
            'FINS';
          'FINS';
*
*         Verifier si les deformations ne sont pas trop importantes
          'SI' IGRD ;
            ref = 20;
*           'SI' ((DPSMAX > 1D-2) 'ET' (DPSMAXP < 1D-3)); REF = 100; 'FINSI';
            'SI' (DPSMAXP * 3 < DPSMAX); ref = 50; DPSMAXP = 1.1 * DPSMAXP ;'FINSI';
            'SI' WTAB.'RECALCUL_SOUSINC' ;
              nsoincrn = MAXI 1 (nsoincrn - 3) (((dpsmax/5d-3) + 1) ENTIER);
              nsoincrn = MINI 20 nsoincrn ;
**            mess 'sous decoupage pour epsi : ' nsoincrn;
            'FINSI' ;
            coefmul = 1.D0 / (DPSMAX + XPETIT);
            'SI' (((coefmul < ref) 'OU'  PASUNIL) 'ET' AUTAUG);
              'MESS' 'annulation du sous-increment: trop grand ' dpsmax;
                si LOG_CNTRL;
                RE_CNTRL2*'RIGIDITE'   = DEFO_IMP DEPSTa  (zdept + zdep1)  ZRAID ;
                DE_CNTRL = VRAI ;
                ZFCONT = RE_CNTRL2 '*' ( ZDEP1) 'EXCO'  'FLX' 'NOID' 'FLX'  ;
**              cmul = MINI  1. ( 1d-4 / DPSMAX);
**              ZFEXT2   = ZFEXT2   + (ZFCONT * cmul);
                ZFLXB    = 'EXCO' ZFEXT2 'FLX' 'NOID' 'FLX' ;
                ZFCONSTA = ZFEXT2 ;
            'SI' (RED_URG < 3);  RE_CNTRL = RE_CNTRL 'ET' RE_CNTRL2; 'FINSI';
                finsi;
*            il faut reduire d'urgence les deplacements
              'SI' (red_urg > 2);
                resmul = resmul * ('MINI' 0.25 (1d-3 / dpsmax));
              'FINSI';
              augmult = augmult * 1.5;
              zdep1d = zdep1 'ENLE' 'LX';
              'SI' (RED_URG > 2) ;
                'SI' ('NON' nonconv);  
                  'MESS' ' non convergence detectee 1' ;
                  ZMAXIT = IT+5;
                'FINS';
                nonconv = vrai;
              'FINS';
**            'SI' (resmul < 0.99); 
**              resmul = resmul * ('MINI' 0.25 (1d-3 / dpsmax));
**            'FINSI';
              'SI' ((dpsmaxp > 1.8d-2) 'ET' (augmult > 100));
                zprecnc= 2.1d-2;
                'SI' (zmaxit > it) ;
                  zmaxit = it - 1;
                'FINSI';
              'FINSI';
              'SI' ((zmaxit > it) 'ET' (RED_URG > 3)); ZMAXIT = IT-1; 'FINS';
              RED_URG = RED_URG + 1;
              'ITERER' etiq;
            'SINON';
**            RE_CNTRL = 'VIDE' 'RIGIDITE' / 'RIGIDITE';
**            DE_CNTRL = FAUX;
            'FINS';
          'FINS';
*
*         --------------------------------------------------------------
*         Loi de comportement sur la configuration definie par LAG_TOT 
*         (Operateur ELAS ou COMP)
*         ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*         - Comportement elastique lineaire
*         ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
          'SI' ('NON' IPLAVI) ;
*
*           Increment de contraintes
            DSIGTa = 'ELAS' ZMODL DEPSTa ZMAT22 ;
*
*           Soustraire les parties issues des chargements
            'SI' (ITHER 'OU' WTAB.'MATVAR' 'OU' LOGDEF) ;
              XXX1 = DSIGM '*' (COEPI '/' nsoincr) ;
              XXX2 = DSIGTa '-' XXX1;
              DSIGTa =  XXX2  ;
            'FINS';
*         ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*         - Comportement non lineaire : integration par l'operateur COMP
*         ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
          'SINON';
*
*  dans le cas des modeles endommageables de Lemaitre, on ecoule
*  en tenant compte, dans les iterations internes, de la variation du
*  materiau avec la temperature
*
*           Soustraire les parties issues des chargements
            'SI' (ITHER 'OU' WTAB.'MATVAR' 'OU' LOGDEF) ;
              XXX1   = DEPSI '*' (COEPI '/' nsoincr) ;
              DEPST  = DEPST  '-' XXX1 ;
              DEPSTa = DEPSTa '-' XXX1 ;
            'FINS' ;
*
            'SI' ('ET' ('ET' ITHER ('OU' IENDOM IVIDOM )) WTAB.'MATVAR') ;
* on recupere certain materiau avec les parametres fct de la temperature
* voir PAS_mate (il ne faut que la dependance thermique)
              ZMATT = 'REDU' WTAB.'MA_COMP' ZMODL;
            'FINS';
*
*           ...cas SSTE...
            'SI' ISSTE;
              ZRIKTA ZSIGFa ZVARF ZDEIF =
                'SSTE' ZMODL ZSIG0a ZVAR0 DEPSTa ZMATT
                       ZPREK WTAB.'NMAXSUBSTEPS' NITMA;
              zvar0 = ('EXCO' zvar0 ('ENLE' com_var ('DIME' com_var))) 
                        'ET' ('EXCO' zvarf ssii);
              zvar0 = 'CHANGER' 'TYPE' zvar0 'VARIABLES INTERNES';
*
*           ...cas non SSTE...
            'SINON';
*
              tsodeb = TDEBUT '+' (ZDT '*' codeb) ;
              tsofin = TDEBUT '+' (ZDT '*' cofin) ;
*
              che1 = 'ADET' 'NOUV' ZMODL CHASANST tsodeb 'TEMP' tsodeb;
              che2 = 'ADET' 'NOUV' ZMODL CHASANST tsofin 'TEMP' tsofin;
*
              'SI' ('EXIS' WTAB 'MOD_LIA') ;
                 vite1 =  'CHAN' 'COMP' conti.'VITESSES' MLPRIM_LIA MVPRIM_LIA ;
                 che1  = che1
                   'ET' ('CHAN' 'CHAM' ZMODL vite1                'STRESSES')
                   'ET' ('CHAN' 'CHAM' ZMODL conti.'DEPLACEMENTS' 'STRESSES');
                 che2  = che2
                   'ET' ('CHAN' 'CHAM' ZMODL  ZDETOT              'STRESSES')
                   'ET' ('CHAN' 'CHAM' ZMODL  ZFCONST1            'STRESSES');
              'FINS';
 
              'SI' ither ;
                TETDE = ZTET1 + (DTETD '*' codeb) ;
                TETDF = ZTET1 + (DTETD '*' cofin) ;
                che3 = 'CHAN' 'CHAM' ZMODL TETDE 'STRESSES' ;
                che4 = 'CHAN' 'CHAM' ZMODL TETDF 'STRESSES' ;
                che1 = che1 'ET' che3 ;
                che2 = che2 'ET' che4 ;
              'SINON' ;
                che1 = 'ADET' che1 ZMODL 'T' 20. ;
                che2 = 'ADET' che2 ZMODL 'T' 20. ;
              'FINS' ;
 
*On met les deformations en tete des champs pour COMP ('KTAN' 'PERT')
              che1 = ZDEF0a 'ET' che1 ;
              che6 = ZDEF0a '+' DEPSTa ;
              che2 = che6 'ET' che2 ;
* pour les materiaux on garde toujours la valeur a la fin du pas
* et au debut du pas
              'SI' ('EXISTE' PRECED 'MODELE_STATIONNAIRE') ;
                 che1 = che1 'ET' ('REDU' ZMAT1 PRECED.'MODELE_STATIONNAIRE') ;
               'SINON' ;
                 che1 = che1 'ET' ZMAT1 ;
              'FINSI' ;  
*
*  Modele NON LINEAIRE UTILISATEUR + GRANDES DEFORMATIONS
*  On ajoute les gradients de deplacements qui seront transformes en
*  gradients de transformation avant appel a UMAT (dans WKUMA1.ESO)
              'SI' (( 'NEG' ZGRDU0 'INCONNU') 'ET' igrd) ;
                gru1 = ZGRDU0 + (D_GR_U '*' codeb) ;
                che1 = che1 'ET' gru1 ;
                gru2 = ZGRDU0 + (D_GR_U '*' cofin) ;
                che2 = che2 'ET' gru2 ;
              'FINS';
*
              'SI' LNLOC ;
                ISTEP = 1 ;
                chm_z = 'ADET' 'NOUV' ZMODL CHASANST tsodeb 'STEP' ISTEP ;
                Che1a = che1 ;
                che1  = che1 'ET' chm_z ;
                chm_z = 'ADET' 'NOUV' ZMODL CHASANST tsofin 'STEP' ISTEP ;
                Che2a = che2 ;
                che2  = che2 'ET' chm_z ;
              'FINS';
*
              che11 = che1 'ET' ZSIG0a 'ET' ZVAR0a 'ET' ZDEI0a ;
              che22 = che2 'ET' ZMATT ;
              'SI' PARTLOCA;
                souci 0 ;
                cho22 = 'ASSIS' 'TOUS' 'COMP' MODRELOC che11 che22 ;
                isoucomp = souci ;
                cho2  = 'REDU' ZMODL ('ETG' cho22) ;
                cho22 = 0 ;
              'SINON';
                souci 0 ;
                cho2  = 'COMP' ZMODL che11 che22 ;
                isoucomp = souci ;
              'FINS';
*
*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*             DEBUT MODELE NON LOCAL
              'SI' LNLOC ;
                ZVARF = 'CHANGER' ('EXCO' cho2 com_var 'NOID') 'TYPE' 'VARIABLES INTERNES' ;
* NLOC ne traitant pas des champs //, on reduit tout au modele initial.
                'SI' PARALLEL ;
                  ZVARF = 'REDU' ZVARF ZMODLI ;
                'FINS' ;
 
********** boucle sous-iterations Helmholtz **************
 
*               donnees pour acceleration Helmholtz 
*               periode de recours a  l acceleration ACT3 pour Helmholtz     
                ZNACCEHL = 3;
*               NBR min de sous-iterations d Helmholtz  a faire avant d appeler ACT3 pour la 1ere fois
                ITDEPHL  = 4;
 
                'REPE' BIHL WTAB.'MAXSOUSITERATION' ;
 
                  'SI' ('NEG'  WTAB.'NON_LOCAL' 'HELM');
*                   cas non-local MOYE ou SB
                    'SI' ('EGA' WTAB.'NON_LOCAL' 'SB') ;
                      MOD_SB =  WTAB.'NLOC_MODL' ;
                      CONTP = 'PRIN' ('REDU' ZSIG0a MOD_SB) MOD_SB ;
                      ZVARF = ZVARF '+' CONTP '+' WTAB.'NLOC_SB_REGU' ;
                    'FINS' ;
                    ZVARN = 'NLOC'  ZVARF  WTAB.'CONN'  ;
                    ZVARN = 'CHANGER' ('EXCO' ZVARN com_var) 'TYPE' 'VARIABLES INTERNES' ;
                    LIERR1='PROG' 0 ; 
                    LOGHLIN = VRAI ;
                  'SINON' ;
 
*                   cas non-local HELM Sellier Millard
                    MOD_HELM =  WTAB.'NLOC_MODL' ;
                    ZVARN = 'COPI' ZVARF ;
                    ZVARF1 = 'REDU' ZVARN MOD_HELM ;
                    LICOHELM = 'EXTR' ZVARF1 'COMP' ;
*                   (re) lecture du materiau si init ou evolutif
                    ZMATHL = PAS_MATE PRECED cho2 ; 
*                   (re) calcul des matrices de rigidite de Helmholtz Sellier        
                    PAS_HELM PRECED ZMATHL;
*                   variable pour la liste des residus de convergence sur Helmholtz
                    LIERR1='PROG'; 
 
*                   ----- debut boucle sur les formulations de Helmholtz ---
*                   on reste en lineaire si toutes les vari sont lineaires (sinon on aura des istep a 3)
                    LOGHLIN = VRAI ;
                    'REPE' BH NHELM ;
                      LOGHLIN=LOGHLIN 'ET' TAHELM . &BH. 'LINEAIRE' ;
*                     nom de la variable a delocaliser
                      LEMO = TAHELM . &BH. 'NOM' ;
*                     recuperation de la variable sur le modele mecanique
                      ZVAUX = 'EXCO' ZVARF1 LEMO 'SCAL';
                      ZVAUX2 = 'CHAN' 'NOEUD' MOD_HELM ZVAUX;
*                     projection sur le modele de Helmholtz associe
                      ZVAUX2 = 'PROI' TAHELM . &BH . 'H_MODELE'  ZVAUX2 'MINI';
*                     TRAC zvaux2 TAHELM . &BH . 'H_MODELE' titre ('CHAI' 'Vari a diffuser:'LEMO );
*                     la variable devient la source
                      FSOUR = 'SOUR' (TAHELM . &BH . 'H_MODELE') ZVAUX2 ;
*                     trac fsour (extr TAHELM . &BH . 'H_MODELE' 'MAIL') titre 'Source';
*                     trac (TAHELM . &BH . 'H_OPER')  titre 'H_Oper';
*                     ***  pour istep=1 *** 1ere etape non locale ***             
                      'SI' ('EGA' &BIHL 1)  ;
*                       &BIHL est le compteur de sous iteration de Helmholtz
*                       il vaut forcement 1 si istep vaut 1, on fait la resolution classique
                        ZVNEW = 'RESO' (TAHELM . &BH . 'H_OPER') FSOUR ;
*                       sauvegarde pour possible acceleration ulterieure 
                        TAHELM . &BH . 'XN' = ZVNEW ;
                      'SINON';
*                       recuperation de la solution de la sous iteration de Helmholtz precedente
                        ZVNEW = TAHELM . &BH . 'XN' ;
                      'FINSI';
*                     *** traitement non lineaire de la formulation ***
                      'SI' ( 'NON'  TAHELM . &BH . 'LINEAIRE' ) ;
*                       calcul du residu sur l evolution de la source         
                        RESN = FSOUR - (TAHELM . &BH . 'H_OPER' '*' TAHELM . &BH . 'XN' );              
*                       stockage du residu pour la mesure d erreur                   
                        RESN0='NOMC' RESN 'SCAL' ;
*                       normalisation des residus par l inverse de la capacite 
*                       (calculee dans pas_helm.procedur)
                        ZERR0 = RESN0 '*'  (TAHELM . &BH . 'INV_CAPA' );
*                       erreur residu
                        ZERR2 = 'MAXI' ( 'ABS' (ZERR0) );
*                       'SI' ('EGA' &BIHL 1);
*                         TAHELM . &BH . 'RES1' = ZERR2;
*                       'FINSI';   
*                       carre de norme du residu 
*                       ZERR1 = 'XTX' ZERR0; 
*                       Residu 
*                       ZERR2 = (ZERR1)**0.5;
*                       normalisation de la source par la capacite
                        ZERR3 = FSOUR '*' (TAHELM . &BH . 'INV_CAPA' )('MOTS' 'Q') ('MOTS' 'SCAL') ('MOTS' 'SCAL');
*                       calcul erreur
                        ZERR6 = 'MAXI' ( 'ABS' (ZERR3) );
*                       carre de la norme de la source
*                       ZERR5 = 'XTX' ZERR3; 
*                       norme de la source
*                       ZERR6 = ZERR5 **0.5;
*                       valeur maximale du champs
*                       ZERR3='MAXI' (TAHELM . &BH . 'XN') 'ABS';
*                       mess 'variable Helmholtz', LEMO;
*                       mess 'residu', ZERR2;
*                       mess 'norme de la source ', ZERR5;                     
*                       mess 'norme Source', ZERR3;
*                       on teste si la valeur maximale est capable de normaliser
                        'SI' (ZERR6 '>' 0.);
                          ZERR4 = ZERR2 '/' ZERR6;
                        'SINON';
                          ZERR4 = ZERR2;
                        'FINS';
*                       preparation acceleration                              
                        'SI' ( 'EGA' &BIHL 1 ) ;         
                          TAHELM . &BH . 'H_RESO' = TAHELM . &BH . 'H_OPER' ;            
                          FWOR = FSOUR * 0. ;
                          TAHELM . &BH . 'ACFP1'  = FWOR ;
                          TAHELM . &BH . 'ACFP2'  = FWOR ;
                          TAHELM . &BH . 'RESNP'  = FWOR ;
                          TAHELM . &BH . 'ACFP3'  = FWOR ;
                          TAHELM . &BH . 'ACFEP1' = FWOR ;
                          TAHELM . &BH . 'ACFEP2' = FWOR ;
                          TAHELM . &BH . 'CORREC' = 0.;
                        'SINON' ;
*                         acceleration de convergence act3
                          CORRECPHL = TAHELM . &BH . 'CORREC';
                          CORRECHL = 0;
                          ACFP0HL = RESN - TAHELM . &BH . 'RESNP' ;
                          ACFEP0HL = ACFP0HL - CORRECPHL ;
                          'SI' ('MULT' &BIHL ZNACCEHL)   ;
                            'SI' (&BIHL > ITDEPHL);
                              CORRECHL   = 'ACT3' 
                              TAHELM . &BH . 'ACFEP2' 
                              TAHELM . &BH . 'ACFEP1'  ACFEP0HL 
                              TAHELM . &BH . 'ACFP3' 
                              TAHELM . &BH . 'ACFP2' 
                              TAHELM . &BH . 'ACFP1'   ACFP0HL ;
                              RESN = RESN - CORRECHL ;
                            'FINS';
                          'FINS';
                          TAHELM . &BH . 'RESNP'  = RESN ; 
                          TAHELM . &BH . 'ACFP3'  = TAHELM . &BH . 'ACFP2'  ; 
                          TAHELM . &BH . 'ACFP2'  = TAHELM . &BH . 'ACFP1'  ; 
                          TAHELM . &BH . 'ACFP1'  = ACFP0HL ;
                          TAHELM . &BH . 'ACFEP2' = TAHELM . &BH . 'ACFEP1' ; 
                          TAHELM . &BH . 'ACFEP1' = ACFEP0HL ;    
*                         calcul correction delta-sigma                               
                          DZVN = 'RESO' TAHELM . &BH . 'H_RESO' RESN;                 
                          ZVNEW = ZVNEW + DZVN ;
                          TAHELM . &BH . 'XN'   = ZVNEW ;
                        'FINSI' ;                          
                      'SINON';
*                       cette variable est traitee de facon lineaire                       
                        ZERR4=0.;
                      'FINSI';
                      ZVNEW1 = 'NOMC' ZVNEW 'SCAL' ;              
                      ZVNEW = 'CHAN' 'CHAM' ZVNEW1 MOD_HELM 'STRESSES' 'VARIABLES INTERNES';
*                     trac zvnew mod_helm titre 'ZVNEW';
                      DZVN = ZVNEW '-' ZVAUX ;
                      DZVN = 'NOMC' DZVN LEMO ;
                      ZVARN = ZVARN '+' DZVN ;
 
*                     mess 'Erreur retenue', ZERR4;
*                     stockage de l erreur normalisee dans la liste des erreurs
                      LIERR1 = LIERR1 'ET' ('PROG' ZERR4);
 
*                     cas particulier TWL2 Sellier theorie WL2
                      'SI' ('EGA' LEMO 'TWL2') ;
                        ZSMAX = 'EXCO' ZVARF1 'SMAX' 'SMAX' ;
                        ZSMAN = 'MANU' 'CHML' MOD_HELM 'SMAX' ('MAXI' ZSMAX) 'TYPE' 'VARIABLES INTERNES' 'STRESSES' ;
                        ZVARN = ZVARN '+' ZSMAN '-' ZSMAX ;
                      'FINS' ;
*                     ---- cas particulier SER Sellier renfort avec gradient
*                     calcul GRADIENT des contraintes sur les renforts (SERi)
*                     stockage dans la VARI HELM1(&BH,2) du gradient selon l absc curviligne 
                      MOPRE = 'EXTR' LEMO 1 3 ;
                      'SI'  ('EGA' MOPRE 'SHR'); 
                        NUME='EXTR' LEMO 4; 
                        NDIM1='VALE' 'DIME';
                        'SI' ( 'EGA' NDIM1 3 );
*                         valeur SERi du pas precedent dans Helm(&bh,1)                             
                          LEMO1=CHAI 'H' &BH '1'; 
*                         valeur SERi du pas actuel 
                          ZAUX1='EXCO' ZVARN LEMO 'SCAL';                 
                          ZVNEW1='NOMC' ZAUX1 'T' ;          
                          ZVAUX1 = 'CHAN' 'NOEUD' MOD_HELM ZVNEW1;              
                          ZVNEW2='PROI' TAHELM . &BH . 'H_MODELE' ZVAUX1 'MINI';
                          GRA1=GRAD TAHELM . &BH . 'H_MODELE' ZVNEW2;                 
                          GRA2=CHAN 'CHAM' 'NOEUD' GRA1 TAHELM . &BH . 'H_MODELE';
                          GRA3=PROI GRA2 MOD_HELM; 
                          LMO3=EXTR GRA3 COMP; 
                          NOM1=EXTR LMO3 1; 
                          NOM2=EXTR LMO3 2;
                          NOM3=EXTR LMO3 3;
                          GRA11=EXCO NOM1 GRA3 'SCAL';
                          GRA22=EXCO NOM2 GRA3 'SCAL';
                          GRA33=EXCO NOM3 GRA3 'SCAL';                  
*                         recuperation vecteur direction de projection
                          LEMON1=CHAI 'VR' NUME '1';
                          NX1=EXCO LEMON1 ZMATHL 'SCAL';
                          NX1=CHAN 'NOEUD' NX1 MOD_HELM;
                          NX11=CHAN NX1 TYPE 'SCALAIRE';
                          LEMON2=CHAI 'VR' NUME '2';
                          NX2=EXCO LEMON2 ZMATHL 'SCAL';
                          NX2=CHAN 'NOEUD' NX2 MOD_HELM;                  
                          NX22=CHAN NX2 TYPE 'SCALAIRE';
                          LEMON3=CHAI 'VR' NUME '3';
                          NX3=EXCO LEMON3 ZMATHL 'SCAL'; 
                          NX3=CHAN 'NOEUD' NX3 MOD_HELM;
                          NX33=CHAN NX3 TYPE 'SCALAIRE';
*                         produit scalaire                  
                          ZVAUX3=(NX11*GRA11)+(NX22*GRA22)+(NX33*GRA33);
*                         stockage du gradient de epse dans YU2j
                          mot4=CHAI 'H' NUME '2'; 
                          zvaux5= 'NOMC' zvaux3 mot4;
                          zvaux6= 'CHAN' 'CHAM' zvaux5 MOD_HELM 'STRESSES' 'VARIABLES INTERNES';
*                         trac zvaux6 MOD_HELM titre '3D H2';                                             
                          ZVAUX7 = 'EXCO' ZVARF1 MOT4 ;
                          DZVN = ZVAUX6 '-' ZVAUX7;
                          ZVARN = ZVARN '+' DZVN ;        
                        'SINON';
*                         valeur SERi du pas precedent dans Helm(&bh,1)                             
                          LEMO1=CHAI 'H' &BH '1';  
                          ZAUX1='EXCO' ZVARF1 LEMO1 'SCAL';                 
                          ZVNEW1='NOMC' ZAUX1 'T' ;          
                          ZVAUX1 = 'CHAN' 'NOEUD' MOD_HELM ZVNEW1;              
                          ZVNEW2='PROI' TAHELM . &BH . 'H_MODELE' ZVAUX1 'MINI';
                          GRA1=GRAD TAHELM . &BH . 'H_MODELE' ZVNEW2;                 
                          GRA2=CHAN 'CHAM' 'NOEUD' GRA1 TAHELM . &BH . 'H_MODELE';
                          GRA3=PROI GRA2 MOD_HELM; 
                          LMO3=EXTR GRA3 COMP; 
                          NOM1=EXTR LMO3 1; 
                          NOM2=EXTR LMO3 2;
                          GRA11=EXCO NOM1 GRA3 'SCAL';
                          GRA22=EXCO NOM2 GRA3 'SCAL';                 
*                         recuperation vecteur direction de projection
                          LEMON1=CHAI 'VR' NUME '1';
                          NX1=EXCO LEMON1 ZMATHL 'SCAL';
                          NX1=CHAN 'NOEUD' NX1 MOD_HELM;
                          NX11=CHAN NX1 TYPE 'SCALAIRE';
                          LEMON2=CHAI 'VR' NUME '2';
                          NX2=EXCO LEMON2 ZMATHL 'SCAL';
                          NX2=CHAN 'NOEUD' NX2 MOD_HELM;                  
                          NX22=CHAN NX2 TYPE 'SCALAIRE';
*                         produit scalaire                  
                          ZVAUX3=(NX11*GRA11)+(NX22*GRA22);
*                         stockage du gradient de epse dans YU2j
                          mot4=CHAI 'H' NUME '2'; 
                          zvaux5= 'NOMC' zvaux3 mot4;
                          zvaux6= 'CHAN' 'CHAM' zvaux5 MOD_HELM 'STRESSES' 'VARIABLES INTERNES';
*                         trac zvaux6 MOD_HELM titr '2DH2';                                            
                          ZVAUX7 = 'EXCO' ZVARF1 MOT4 ;
                          DZVN = ZVAUX6 '-' ZVAUX7;
                          ZVARN = ZVARN '+' DZVN ;        
                        'FINSI';     
                      'FINSI';
*                        ------fin du calcul des gradients  et Hessiens--------
                    'FIN' BH ;
*                 ------fin boucle sur les formulations de Helmholtz ---
                  'FINS' ;
 
*                 On rend ZVARN // si besoin :
                  'SI' PARALLEL ;
                    ZVARN = 'REDU' ZVARN ZMODL ;
                  'FINS' ;
 
*                 Recherche du residu maxi parmi les formulations non locales 
                  ISTEP0=ISTEP;
*                 residu maxi de tous les Helmholtz actifs
                  ZERR6='MAXI' LIERR1;
 
                  'SI' ('EGA' &BIHL WTAB.'MAXSOUSITERATION') ;
                    'MESSAGE' ' SOUS ITERATION HELMHOLTZ CONVERGENCE FORCEE';
*                   on indique au modele de mecanique que la sous 
*                   iteration de Helmholtz est achevee                 
                    ISTEP = 2 ;
                  'SINON' ;
                    'SI' ((( ZERR6 '<' WTAB.'PRECSOUSITERATION') 'ET' ('NEG' ISTEP 1)) 'OU' LOGHLIN ) ;
*                     on indique au modele de mecanique que la sous 
*                     iteration de Helmholtz est achevee 
                      ISTEP = 2 ;
                    'SINON';
*                     on indique au modele de mecanique que les sous
*                     iterations de Helmholtz continuent
                      ISTEP = 3 ;
                    'FINS';
                  'FINS' ;
*
                  'SI' ('EGA' ISTEP0 1);
                    mess1=chai ' SOUS ITERATION HELMHOLTZ:' &BIHL  ' ISTEP:' ISTEP0 ;
                  'SINON';
                    mess1=chai ' SOUS ITERATION HELMHOLTZ:' &BIHL  ' ISTEP:' ISTEP0 ' CRITERE:' ZERR6;
                  'FINSI';
                  'MESSAGE' mess1 ;           
*
                  chm_z = 'ADET' 'NOUV' ZMODL CHASANST tsodeb 'STEP' ISTEP ;
                  che11 = CHE1A 'ET' chm_z ;
                  che11 = che11 'ET' ZSIG0a 'ET' ZVARN 'ET' ZDEI0a ;
 
                  chm_z = 'ADET' 'NOUV' ZMODL CHASANST tsofin 'STEP' ISTEP ;
                  che22 = CHE2A 'ET' chm_z ;
                  che22 = che22 'ET' ZMATT ;
 
                  'SI' PARTLOCA ;
                    souci 0 ;
                    cho22 = 'ASSIS' 'TOUS' 'COMP' MODRELOC che11 che22 ;
                    isoucomp = souci ;
                    cho2  = 'REDU' ZMODL ('ETG' cho22) ;
                    cho22 = 0 ;
                  'SINON';
                    souci 0 ;
                    cho2  = 'COMP' ZMODL che11 che22 ;
                    isoucomp = souci ;
                  'FINS' ;
*                         
                  'SI' ('EGA'  ISTEP 2 );
                    'QUITTER' BIHL;
                  'FINS';
*                 fin de boucle sous-iterations 
                  ZVARF = 'CHANGER' ('EXCO' cho2 com_var 'NOID') 'TYPE' 'VARIABLES INTERNES' ;
*                 NLOC ne traitant pas des champs //, on reduit tout au modele initial.
                  'SI' PARALLEL ;
                    ZVARF = 'REDU'  ZVARF ZMODLI ;
                  'FINS' ;                
                'FIN' BIHL ; 
*      ------ fin des sous iterations de Helmholtz ------
 
*               Un peu de menage
                CHE1A = 1 ;
                CHE2A = 1 ;
                ZVARN = 1 ;
                ZVARF = 1 ;
                chm_z = 1 ;
                LIERR1 = 1 ;      
                ZMATHL = 1; 
              'FINS' ;
*             FIN MODELE NON LOCAL
*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*
              si isoucomp ;
                mess ' Souci dans le comportement' ;
              fins ;
              ZSIGFa = 'CHANGER' ('EXCO' cho2 com_sig 'NOID')'TYPE' 'CONTRAINTES' ;
              ZVARF  = 'CHANGER' ('EXCO' cho2 com_var 'NOID')'TYPE' 'VARIABLES INTERNES' ;
              ZDEIF  = 'CHANGER' ('EXCO' cho2 com_dei 'NOID')'TYPE' 'DEFORMATIONS INELASTIQUES' ;
*
              'SI' ('EXIS' WTAB 'MOD_LIA') ;
                ZFLIA = 'EXCO' cho2 'FLIA' 'NOID';
              'FINS' ;
*
              'SI' (IDYN 'ET' ('NEG' nsoincr 1)) ;
                'SI' ('EGA' &sousinc nsoincr) ;
                  aaa1=zsigfa/2;
                  aaa2=aaa1+ zsigm;
                  'DETR' aaa1; 'DETR' zsigm;
                  zsigm=aaa2;
                  CONT = ZSIGFa + ZSIGM0 ;
                'SINON';
                  aaa2=zsigfa+ zsigm;
                  'DETR' zsigm;
                  zsigm=aaa2;
                'FINS';
              'FINS';
 
*     le fin d'en dessous est le fin de "si isste ... sinon ... finsi"
            'FINS';
*
*           cas particulier poreux et thermique
            'SI' (ITHER 'ET' WTAB.'POR1');
              ZSIGFa = ZSIGFa - ( (COEPI '/' nsoincr) '*' DMSRT0 ) ;
            'FINS';
*
            'SI' ('NON' WTAB.'MODAL') ;
              DSIGTa = ZSIGFa '-' ZSIG0a ;
            'FINSI' ;
*
*...fin du si ISSTE sinon ...
          'FINSI' ;
*
*         Transporter les contraintes sur GEOREF0
          'SI' ('NON' WTAB.'MODAL') ;
            'SI' IGRD ;
              'FORM' GEOREF0 ;
              DSIGTa = 'CONF' DSIGTa ZMODL  ;
            'FINS';
            DSIGTZ = DSIGT '+' DSIGTa ;
            'DETR' DSIGTa ; 'DETR' DSIGT ;
            DSIGT  = DSIGTZ ;
            ZSIGFa = ZSIG0 '+' DSIGT ;
          'FINS';
*
          'SI' ('NEG' &sousinc 1) ;
            'DETR' ZSIG0a ; 'DETR' ZDEF0a ; 'DETR' ZVAR0a ; 
            'DETR' ZDEI0a ;
          'FINSI' ;
          'DETR' DEPSTa ;
*
          ZSIG0a = ZSIGFa;
          ZDEF0a = ZDEFF ;
          ZVAR0a = ZVARF ;
          ZDEI0a = ZDEIF ;
*
        'FIN' sousinc;
*       ---------------------------------------------------------------
*
*       Exprimer sur la configuration GEOREF0
        ZSIGF  = ZSIGFa ;
        DPSMAX = 'MAXI' 'ABS' DEPST 'AVEC' MLDEFOR ;
**         DPSPRI = DEPST 'PRIN' ZMODL ZVAR0;
**         DPSMAX = 'MAXI' 'ABS' DPSPRI 'AVEC' MLPRINC;
**         'DETR' DPSPRI;
        EPSM   = 'MAXI' 'ABS' ZDEFF 'AVEC' MLDEFOR ;
**        EPSPRI = ZDEFF 'PRIN' ZMODL ZVAR0;
**        EPSM = 'MAXI' 'ABS' EPSPRI 'AVEC' MLPRINC;
**        'DETR' EPSPRI;
*
* pour tenir compte de ce que le travail de la correction est 1/2 FU,
* on la multiplie par 2
        'SI' (IDYN 'ET' ('NEG' nsoincr 1));
          ZSIGM = ZSIGM*(2. /nsoincr) ;
          ZSIGM = ZSIGM - CONT ;
          BZSIGM='BSIG' ZMODL ZSIGM ZMAT22 ;
        'FINS';
 
        'SI' IPLAVI ;
*
*   Matrice tangente par perturbation evaluee pour la derniere iteration calculee
*   A voir : cas grand deplacement ZSIGF et ZMAT2, cas poreux et thermique ZSIGF
          'SI' (IKTAN 'ET' IPERT) ;
            Z1COMP = che11 ;
            Z2COMP = che22 'ET' ZSIGFa ;
          'FINS' ;
*
*         Max de la composante 'MOCA' (='EPSE' par defaut) pendant l'iteration
*         nbre de points qui ont une evolution non lineaire
          'SI' ('NEG' WTAB.'MOVA' 'RIEN') ;
            XXX1 = 'EXCO' (WTAB.'MOVA') ZVARF ;
            ACC = 'ABS' ( XXX1 - ACC0 ) ;
            MMC = 'MASQUE' ACC 'SUPERIEUR' 1.D-10 'SOMME' ;
            'SI' (MMC > MMCMAX) ; MMCMAX = MMC ; 'FINS' ;
          'FINS';
        'FINS';
*
* temporaire verification zsigf
**   zsigtmp = 'ELAS' ZMODL zdeff ZMAT22 ;
**   zdiff = (zsigf - zsigtmp);
**   zdmax = maxi abs zdiff;
**   zdref= maxi abs zsigf;
**   mess 'zdmax:' zdmax ' zdref:' zdref;
**        mess 'iraug rate' ' ' iraug ' ' irate;
**        zsigf = zsigtmp;
*
*----------------------------------------------------------------------*
*          Forces nodales equivalentes au champ de contraintes         *
*----------------------------------------------------------------------*
        'SI' IGRD ;
          'FORM' GEOM1 ;
          'SI' ('NON' HPP_EPS) ;
            'SOUC' 0;
*           Changer de GEOREF0 -> GEOM2
            'SI' IJAUMA; 
              ZSIGG = 'CONF' ZSIGF ZMODL ;
              'FORM' GEOM2;
              ZSIGG = 'CONF' ZSIGG ZMODL ;
            'SINON';
              'FORM' GEOM2 ;
              ZSIGG = 'CONF' ZSIGF ZMODL ;
            'FINSI';
*HHO : Modifications appel a BSIGMA (ajout du champ de deplacements necessaire en HHO)
            FEQU2 = 'BSIG' 'NOER' ZMODL ZSIGG ZMAT22 ZDETOT ;
*
            'SI' IJAUMA;
              'FORM' GEOREF0 ;
              ZSIGF = 'CONF' ZSIGG ZMODL ;
            'FINSI' ;
*
            'FORM' GEOM1;
            'SI' (SOUCI) ;
              ZSIGG = 'CONF' ZSIGG ZMODL  ;
              MESS 'BSIG impossible : tentative sur geom1' ;
*HHO : Modifications appel a BSIGMA (ajout du champ de deplacements necessaire en HHO)
              FEQU2 = 'BSIG' 'NOER' ZMODL ZSIGG ZMAT21 ZDETOT ;
              pastest = vrai;
            'FINSI';
          'SINON';
*HHO : Modifications appel a BSIGMA (ajout du champ de deplacements necessaire en HHO)
*    config debut de pas pour predicteur hpp
            ZSIGG = 'CONF' ZSIGF ZMODL ;
            FEQU2 = 'BSIG' ZMODL ZSIGG ZMAT22 ZDETOT;
          'FINSI';
        'SINON';
*         config unique car HPP (igrd faux)                     
          FEQU2 = 'BSIG' ZMODL ZSIGF ZMAT22 ZDETOT ;
        'FINS';
*
        'SI' IRCON;
          FEQU2 = FEQU2 '+' ('REDU' MAI_CONS (RIG_CONS '*' ZDETOT));
        'FINS';
*----------------------------------------------------------------------*
*          Mise a jour des chargements dependants de la geometrie      *
*----------------------------------------------------------------------*
        iACTU = FAUX ;
        'SI' (LOGPRE 'ET' IGRD 'ET' ('NON' HPP_EPS)) ;
          'FORM' GEOM2 ;
          'SI' (('NON' IPILOT) 'ET' (COEPI '<' 1.)) ;
            ZPEXTF = 'TIRE' WTAB.'CHARGEMENT' 'PRES' TFIN ;
          'FINSI' ;
          'SI' ('EXIS' WTAB 'MAT_PRE') ;
            ZFPEXTF = 'BSIG' WTAB.'MOD_PRE' ZPEXTF WTAB.'MAT_PRE' ;
          'SINON' ;
            ZFPEXTF = 'BSIG' WTAB.'MOD_PRE' ZPEXTF ;
          'FINS' ;
          ZFSUIV = ZFPEXTF ;
          iACTU  = VRAI ;
          'FORM' GEOM1 ;
        'FINS' ;
*
        'SI' WTAB.'PROCEDURE_CHARMECA';
          'SI' (IGRD 'ET' ('NON' HPP_EPS)) ;
            'FORM' GEOM2 ;
          'FINSI' ;
          TFP22    = CHARMECA PRECED TFIN ;
          ADDISEC2 = 'EXIS' TFP22 'ADDI_SECOND' ;
          'SI' ADDISEC2 ;
            FP22   = TFP22.'ADDI_SECOND' ;
            'SI' iACTU ;
              ZFSUIV = ZFSUIV '+' FP22 ;
            'SINON';
              ZFSUIV = FP22 ;
            'FINS';
          'FINS';
          'FORM' GEOM1 ;
        'FINS' ;
*
      'FINS';
*     Finsi de Update or total lagrangian (IFEFP) ---------------------
*
      'SI' ('EXIS' WTAB 'MOD_LIA') ;
        FLIAI = 0.D0 ;
        'SI' ('EGA' ('TYPE' ZFLIA) 'MCHAML  ') ;
          NZLIA = 'EXTR' ZFLIA 'NBZO' ;
          'SI' (NZLIA > 0) ;
            'REPETER' BZLIA NZLIA ;
* un point support par zone - 2010 kich
              FLIAI= FLIAI +
  ('EXCO' ('EXTR' ZFLIA 'FLIA' &BZLIA 1 1) MLDUAL_LIA 'NOID') ;
            'FIN' BZLIA ;
          'FINS' ;
          FEQU2 = FEQU2 + FLIAI ;
        'FINS' ;
      'FINS' ;
*
      'SI' ISOL  ;
        'SI' IPLAVI ;
          ZMATPF = 'CHANGER' ('EXCO' cho2 com_maa 'NOID')
          'TYPE' 'CARACTERISTIQUES';
          MA_POR = 'REDU' ZMATPF WTAB.'MOD_POR' ;
        'FINS';
        XXX1 = 'GRAD' WTAB.'MOD_POR' ZDETOT MA_POR 'CONS'  ;
        XXXS = (1. - WTAB.'TETA' )*GRAP0 + (WTAB.'TETA'*XXX1) ;
        XXX2 = ZDT * ('GNFL' WTAB.'MOD_POR' XXXS) ;
        XXX3 = FEQU2 ;
        FEQU2 = XXX3 - XXX2;'DETR' XXX3; 'DETR' XXX2 ;
      'FINS' ;
*
* si on a fait de la sous incrementation ( a cause du dynamique
      'SI' (IDYN 'ET' ('NEG' nsoincr 1));
        FEQU2 = FEQU2 + BZSIGM ;
      'FINS';
    'FINSI' ;
*FOR_MECA
*----------------------------------------------------------------------*
*   Contributions propres a la dynamique intervenant dans le residu    *
*                       [ 4/h2*m + 2/h *c ] * du                       *
*----------------------------------------------------------------------*
    'SI' IDYN;
      'DETR' FFDYN;
      FFDYN = 'COPIER' FEQU2 ;
      XXX1 = WTAB.'MASSE' * ZDEPT;
      XXX3 = 4. * UNSURH * UNSURH  * XXX1 ;
      'DETR' XXX1;
      'SI' ('EXIS' WTAB 'MOD_LIA');
        'SI' (NZLIA > 0) ;
          XVIT2 = 0.D0;
          'REPETER' BZLIA NZLIA ;
            XVIT2 = XVIT2 +
       ('EXCO' ('EXTR' ZFLIA 'FLIA' &BZLIA 1 1) MVPRIM_LIA MLPRIM_LIA 'NOID') ;
          'FIN' BZLIA ;
        'FINS' ;
      'FINS' ;
      'SI' ( 'NEG' WTAB.'AMORTISSEMENT' 'INCONNU') ;
        XXX1 = WTAB.'AMORTISSEMENT' * ZDEPT;
        XXX2 = 2. * UNSURH * XXX1 ;
        XXX4 = XXX3 + XXX2 ;
        XXX3 = XXX4 ;
      'FINS' ;
      XXX4 = FEQU2 + XXX3;
      'DETR' XXX3 ; 'DETR' FEQU2 ;
      FEQU2 = XXX4 ;
*
*         forces correctrices en cas de liaison persistante :
*  on veut avoir (forces inertielles + forces visqueuses) compatibles
*  avec accelerations et vitesses relatives nulles aux points de contact
*  (pendant le contact). On modifie fequ2 ---> residu et l'iteration sui
*  fournira les bonnes reactions
*
      'SI' IMPLP;
        VADD XXX3 = VITETFOR ZRAID_T WTAB FCORU ;
        XXX4      = FEQU2 - XXX3;
        'DETR' XXX3 ; 'DETR' FEQU2 ;
        FEQU2     = XXX4;
      'FINS' ;
    'FINS';
*----------------------------------------------------------------------*
*             Calcul du nouveau coefficient pour le pilotage           *
*----------------------------------------------------------------------*
    'SI' WTAB.'FOR_MECA' ;
      'SI' IPILOT ;
        'SI' ACCEL ;
*
          COEPI0 = COEPI;
          'SI' ('OU' IVISCO IVIDOM ITHER LOGDEF);
            'MESS' 'ALPHA calcule avec la norme de l increment' ;
            COEPI = 'MINI' (COEPI0 * AL1) 1.D0 ;
          'SINON' ;
            'SI' ((COEPI0 'EGA' 1.D0) 'ET' (AL1 'EGA' 1.D0));
              COEPI = 1.D0;
            'SINON' ;
              XXX1 = ZFEXT '+' (COEPI '*' ZFSUIV) '-' FEQU2 ;
              XXX2 = 'ENLE' RESIDNOR 'FLX' ;          
*
              XX1 = 'XTY' XXX1 ZDEPT MNDUAL MNPRIM;
              XX2 = 'XTY' XXX2 ZDEPT MNDUAL MNPRIM;
              COEINC =  XX1 '/' XX2;
*
              COEPI = 'MINI' (COEPI0 - COEINC) 1.D0 ;
            'FINS';
          'FINS' ;
*
*         Mise a jour des termes force de l'acceleration de convergence
          XXXZ = RESIDNOR ;
          'SI' WTAB.'PROCEDURE_CHARMECA';
            'SI' ADDISEC2 ;
              XXXZ = RESIDNOR '+' (FP22 '-' FP022) ;
            'FINS';
          'FINS';
          COEINC = COEPI0 '-' COEPI ;
          XXX1   = 'ENLE' (COEINC '*' XXXZ) 'FLX' ;
          XXX3   = ACFEP2 - XXX1 ; 'DETR' ACFEP2; 
          ACFEP2 = XXX3;
          XXX4   = ACFEP1 - XXX1 ; 'DETR' ACFEP1; 
          ACFEP1 = XXX4;
          'DETR' XXX1;
*
*         Actualisation des forces/deplacements
          DFEXT = COEPI '*' DFEXT0F ;
          ZFEXT = DFEXT '+' FEXT0 ;
          'SI' IMPO12;
            DUIMP = (COEPI '*' DUIMPO) '+' DUUNIL ;
          'SINON' ;
            DUIMP = COEPI '*' DFEXT0L ;
          'FINS';
          ZFLX1 = DUIMP '+' FLXREC ;
          ZFSUIV = COEPI '*' ZFSUIV ;
*
        'FINS' ;
      'FINS' ;
*----------------------------------------------------------------------*
*          Equilibre (RESIDC) et nouveau second membre (RESIDU)        *
*----------------------------------------------------------------------*
*     Forces exterieures (sans les reactions)
      XFEXT = ZFEXT '+' ZFSUIV ;
*     Modification du chargement
      'SI' LOGPIL ;      
        XFEXT = XFEXT '+' (ETA '*' ZFPILIN);      
      'FINSI' ;
*
*     Forces exterieures - Forces interieures
      XFORC = XFEXT '-' FEQU2;
*
*     Mise a jour des jeux (FLX) qui travaillent en incremental :
*     - pour les depl imposes cela revient a imposer un increment nul
*     - pour les rela unilaterales cela revient a mettre a jour le jeu
      XJEUX  = ZFLX1 '-' FCORU ;
*
      RESIDU = XFORC '+' XJEUX ;
      RESIDC = XFORC '-' FCORF ;
*----------------------------------------------------------------------*
*             Critere pour statuer de la convergence ou non            *
*----------------------------------------------------------------------*
      XNUMF = 'MAXI' RESIDC 'ABS' 'SANS' MXMFLX;
      XNUMG = 'MAXI' ((ZDEPL '-' ZDEPLP) '*' ZCLIM) 'ABS';
      XNUMF = 'MAXI' ('PROG' XNUMF XNUMG);
*
      'SI' IFTOL;
        ZPREC = ZFTOL ;
        XDENO = 1.;
      'FINS';
      'SI' IMTOL;
        ZPRECM = ZMTOL ;
        XDENOM = 1.;
      'FINS';
*
      XCONV  = XNUMF '/' XDENO ;
      XCONVM = XCONV ;
*
      'SI' TSTMOM ;
        ZFAU1 = XFEXT '-' FCORF '-' ZFPLO ;
        XAUXM = 'MAXI' ZFAU1  'ABS' 'AVEC' MXMYMZ;
        'SI' ('<' XAUXM XPETIT) ;
          XCONVM = 0.;
        'SINON';
          XAUXF = 'MAXI' ZFAU1  'ABS' 'SANS' MXMFLX;
          XNUMM = 'MAXI' RESIDC 'ABS' 'AVEC' MXMYMZ;
          'SI' ('<' (XAUXF / XAUXM) 1.D-12) ;
            XCONV = 0.;
          'FINS';
          XCONVM = XNUMM / XDENOM ;
        'FINS';
      'FINS' ;
*
    'FINSI' ;
*FOR_MECA
*
    'SI' WTAB.'NVSTNL'                               ;
        ZFLCONV = FAUX                               ; 
        mofofl = 'EXTR' ZMODL 'FORC   '     ;
          dvit2 = 'XTX'  ('EXCO' ZVIF1 movifl)       ;
          err  = dvit2 '**' 0.5                      ;
          res2 = 'XTX' ('EXCO' Resnsl mofofl 'NOID') ; 
          errres2  = res2 **0.5                      ;
      'MESSAGE' ' ite:' IT*26 ' errdu:' err*52 ' errres:' errres2*78  ; 
       'SI' (err < 1.d-9); ZFLCONV = VRAI ;   'FINSI'; 
       'SI' (ZFLCONV 'OU' (it '>EG' WTAB.'MAXSOUSPAS')) ; 
                             'QUITTER' ETIQ ; 'FINSI' ;
    'FINSI' ;
*----------------------------------------------------------------------*
*                resume de l'iteration                                 *
*----------------------------------------------------------------------*
    'SI' WTAB.'FOR_MECA' ;
      'SI' IPILOT ;
        'MESS' IT*13 MMC*26 XCONV*39 DPSMAX*52 EPSM*65 XCONVM*78 COEPI*91;
      'SINON' ;
        'MESS' IT*13 MMC*26 XCONV*39 DPSMAX*52 EPSM*65 XCONVM*78;
      'FINSI' ;
*       sauver depstp pour le test increment de deformation
      depstps = depstp;
*       acceptation de l'etat courant pour le sous pas suivant
      depstp zdeptp zsigfp fcorfp correcp geom2p zmat22p = depst zdept zsigf fcorf correc geom2 zmat22;
      IRATE = FAUX;
*       premiere iteration prise en compte
      IPREM = FAUX;
*----------------------------------------------------------------------*
*          test de convergence : equilibre de la structure             *
*----------------------------------------------------------------------*
      TABCONV.IT=XCONV;
*
*     A CONSERVER ????
      'SI' IPILOT;
        'SI' (('EGA' MMC 0) 'ET' (MMCMAX > 0) 'ET' (COEPI < 0.));
*       on refuse de converger si on est elastique et en decharge
          COEPI = 'ABS' COEPI; PASTEST = VRAI;
        'FINSI' ;
      'FINSI';
*
*     Si frottement, faire au moins 2 iterations
      'SI' (WTAB.'CAFROTTE' 'ET' IMPO12 'ET' (IT < 2 ));
        PASTEST = VRAI;
      'FINS';
*
      PASUNIL = FAUX;
      'SI' ('EXIS' ZRAID_T 'OK');
        PASUNIL = 'NON' ZRAID_T.'OK' ;
      'FINS';
*
      'SI' ('NON' PASTEST) ;
*
        'FORM' GEOREF0;
*       Variation de Despi entre 2 itérés
        DEPSTD  = DEPST - DEPSTPS ;
        DEPSTDM = 'MAXI' 'ABS' DEPSTD 'AVEC' MLDEFOR;
*
*       Pas de test sur la convergence a la premiere iteration
*       ou apres une initialisation a partir de la solution precedente
        'SI' ((IT > 1) 'OU' ('NON' INIT)) ;
*
*         Si les criteres (deplacements + moments) sont < precision
*         et la variation sur de epsi est < precision souhaitee : on a convergé!
          ZPRECHPP = ZPREC; 
          ZPRECHPPE = ZPRECEPS;
          'SI' HPP_EPS; ZPRECHPP = ZPREC * 1d-1; ZPRECHPPE = ZPRECEPS * 1D-1; 'FINSI';
*
          iCONV = (XCONV '<' ZPRECHPP) 'ET' (DEPSTDM '<' ZPRECHPPE) ;
          'SI' TSTMOM ;
            iCONV = iCONV 'ET' (XCONVM '<' ZPRECM) ;
          'FINSI' ;
*
          'SI' iCONV ;
*
            'SI' PASUNIL ;
              'MESS' ' ****** NON CONVERGENCE DE CONTACT A L''ITERATION' ' ' IT ' SOUS-PAS' ' ' WTAB.'ISOUSPAS';
              ZICONV = FAUX;
              'QUITTER' ETIQ;
            'FINS';
*
            'SI' (('NON' IPILOT) 'ET' (COEPI '<' 1.));
              'MESS' ' ****** NON CONVERGENCE DUE A LA RELAXATION A L''ITERATION' ' ' IT ' SOUS-PAS' ' ' WTAB.'ISOUSPAS';
              ZICONV = FAUX;
              'QUITTER' ETIQ;
            'FINS';
*
            'SI' ((IRAUG 'ET' AUTAUG 'ET' ((DPSMAX > ZPRECEPS) 'OU' (augmult > 0.6e-1)) 'OU'
                  (resmul < 0.99) 'ET' VRAI) 'OU' DE_CNTRL);
              'MESS' ' ****** NON CONVERGENCE DUE A L AUGMENTATION A L''ITERATION' ' ' IT ' SOUS-PAS' ' ' WTAB.'ISOUSPAS';
*              'MESS' 'augmult dpsmax avant modif' ' ' augmult dpsmax;
              ZICONV = FAUX;
              augmult = augmult * 0.1;
              'SI' (DPSMAX < ZPRECEPS) ; IRAUG = FAUX; 'FINSI';
              'QUITTER' ETIQ;
            'FINS';
*
            iRECA = FAUX ;
            'SI' (WTAB.'CONTACT' 'ET' ('NON' HPP_EPS) 'ET' ('NON' NONCONV));
              'SI' ('NON' (WTAB.'MODAL' 'OU' WTAB.'FROCABL')) ;
                MODCON = WTAB.'MODCONTA';
                'FORM' GEOM2 ;
                'SI' ('EXIS' WTAB 'MATCONTA') ;
                  CJEU CRRN = 'RFCO' MODCON WTAB.'CONV' WTAB.'MATCONTA';
                'SINON' ;
                  CJEU CRRN = 'RFCO' MODCON WTAB.'CONV'                ;
                'FINS';
                'FORM' GEOM1 ;
                ZDLX = 'EXCO' ZDEPT 'LX' 'NOID' 'LX' ;
* 
                REA1 = 'REAC' CRR  ZDLX ;
                CMPF = 'EXTR' REA1 'COMP' ;
                PSC1 = 'PSCA' REA1 CMPF REA1 CMPF ;
                NOR1 = 'MAXI' 'ABS' PSC1 ;
*
                'SI' ('>' NOR1 XDENO) ;
                  REA2 = 'REAC' CRRN ZDLX ;
                  DREA = REA2 '-' REA1 ;
                  PSC2 = 'PSCA' DREA CMPF DREA CMPF ;
*
                  RAP4 = ('MAXI' 'ABS' PSC2) '/' NOR1 ;
*
                  'SI' ('>EG' RAP4 ZPREC) ;
                    'MESS' 'Actualiser les conditions de contact-frottement' rap4 nor1;
*                   En gardant les statuts courants
*                   ZRIBLO_M = 'MOT' 'INCONNU' ;
                    iRECA = VRAI ;
                  'FINSI' ;
                'FINSI' ;
              'FINSI' ;
            'FINSI' ;
*
*           Tout semble ok mais il reste encore a verifier que
*           - le resultat est obtenu sans predicteur HPP
*           - les relations de cont-fro sont "coherentes" avec la conf finale
            'SI' (HPP_EPS 'OU' isoucomp 'OU' iRECA) ;
              'SI' (IGRD 'ET' HPP_EPS) ; 
                'MESS' ' passage en grands deplacements'; 
              'FINSI' ;
*
              HPP_EPS = FAUX;
              URG     = VRAI;
              ITACC   = 4;
*
              DLTAIT  = WTAB.'DELTAITER' '+' IT ;
              PASTEST = VRAI;
              'SI' autaug; iraug = faux; de_cntrl = faux ; 'FINSI';
            'SINON';
              'MESS' ' ****** CONVERGENCE A L ITERATION' ' ' IT '  SOUS-PAS' ' ' WTAB.'ISOUSPAS';
              'MESS' ' ';
*             Nombre d'iteration avant convergence
              ZNITE = ZNITE '+' IT;
              'QUITTER' ETIQ ;
            'FINS';
          'FINS';
        'FINS';
      'FINS' ;
*----------------------------------------------------------------------*
* test de non convergence
*----------------------------------------------------------------------*
      XCONVREF=1E50;
      'SI' (IT > DLTAIT);
        XCONVREF=TABCONV.(IT-DLTAIT) * 0.99;
      'FINS';
*
*   si on a depasse le nombre max d'iterations ou si le residu augmente
*   ou si on aurait du converger et que cela n'est pas le cas :
*   => non convergence detectee !
      'SI' (('NON' PASTEST) 'ET' ('NON' NONCONV) 'ET' (IT > 1) 'ET' 
            ('NON' IPILOT)  'ET' ('NON' PASUNIL)) ;
        'SI' ((IT '>EG' ZMAXIT) 'OU' (XCONV '>' XCONVREF) 'OU' iCONV);
          'MESS' ' non convergence detectee 2' ;
          ZMAXIT  = 3 * IT;
          NONCONV = VRAI;
          PASTEST = VRAI;
*  test pv
*
          HPP_EPS = FAUX;
          'SI' (ITURG > 3); URG     = VRAI; ITACC = 4; 'FINSI';
*
          COEPI = WTAB.'RELAXATION_NONCONV';
*
*         Modification de l'objectif
          'SI' (COEPI '<' 1.) ;
            ITACC = 4 ;
*
*           - Nouvel instant objectif
            TFIN = TDEBUT '+' (ZDT '*' COEPI) ;
*
*           - Reevaluer le chargement final
            DFEXT = COEPI '*' DFEXT0F ;
            ZFEXT = DFEXT '+' FEXT0 ;
            'SI' IMPO12;
              DUIMP = (COEPI '*' DUIMPO) '+' DUUNIL ;
            'SINON' ;
              DUIMP = COEPI '*' DFEXT0L ;
            'FINS';
            ZFLX1 = DUIMP '+' FLXREC ;
*
*           - Reevaluer le residu
            XFEXT  = ZFEXT '+' (COEPI '*' ZFSUIV);
            XFORC  = XFEXT '-' FEQU2;
            XJEUX  = ZFLX1 '-' FCORU ;
            RESIDU = XFORC '+' XJEUX ;
*
*           - Reevaluer le materiau
            'SI' WTAB.'MATVAR' ;
              'SI' IGRD ; 'FORM' GEOREF0 ; 'FINSI' ;
              ETATI  = PAS_ETAT PRECED TFIN ;
              ZMAT2I = PAS_MATE PRECED ETATI;
              ZMAT21 = 'REDU' ZMAT2I ZMODL ;
              'SI' IGRD ; 
                'FORM' GEOM1 ;
                ZMAT21 = 'CONF' ZMODL ZMAT21 ;
              'FINSI' ;
            'FINSI' ;
          'FINSI' ;
*
        'FINS';
      'FINS';
*
*     Changement de la precision en non convergence
      'SI' nonconv ;
*       'SI' (it > (zmaxit *2 /3));
        'SI' (it > 15); 
          zprecnc=zprecnc*2;
        'FINS';
      'FINS';
*
      'SI' (('NON' PASTEST) 'ET' (dpsmax < zprecnc));
        'SI' ((IT '>EG' ZMAXIT) 'OU' 
              (NONCONV 
               'ET' (XCONV '>' XCONVREF)
               'ET' (DPSMAX '>EG' (DPSMAXP /4 ))) );
          'MESS' ' ****** NON CONVERGENCE A L''ITERATION' ' ' IT ' SOUS-PAS' ' ' WTAB.'ISOUSPAS';
          ZICONV = FAUX;
          'QUITTER' ETIQ;
        'FINS' ;
      'FINS' ;
*----------------------------------------------------------------------*
*     Menage d'objets temporaires
*----------------------------------------------------------------------*
      'DETR' FEQU2; 'DETR' ZDEFF;
*
*     Supprimer les configurations intermediaires
      'SI' (IGRD 'ET' ('NON' HPP_EPS));
        'DETR' GEOM2 ;
      'FINS' ;
*
      'SI' IFEFP ; 
        'DETR' ZRIKTA ; 
        'SI' IGRD; 'DETR' GEOM2;  FINSI ;
      'FINS' ;
*
    'FINSI' ;
*FOR_MECA
*
    dpsmaxp = dpsmax;
*
  'FIN' ETIQ ;
*----------------------------------------------------------------------*
*           Fin de la boucle de convergence ETIQ                       *
*----------------------------------------------------------------------*
*######################################################################*
  'SI' IDYN  ;
    VITI = UNSURH * 2. * ZDEPT - conti. 'VITESSES' ;
    VITI = 'ENLE' VITI 'LX' ;
    ZFP  = ZDYFEXT '+' ZFSUIV '-' FFDYN;
*
    'SI' ('EXIS' WTAB 'MOD_LIA');
      'SI' ('EXIS' WTAB.'MOD_LIA' 'MATE' 'NEWMARK') ;
        'SI' ('EGA' 'CHPOINT  ' ('TYPE' FLIAI) ) ;
          MODNEW = 'EXTR' WTAB.'MOD_LIA' 'MATE' 'NEWMARK' ;
          MAUNEW = 'EXTR' MODNEW 'MAIL' ;
          'SI' ('EXIS' ZFLIA MODNEW) ;
            ZFLIANEW = 'REDU' ZFLIA MODNEW ;
            NZNEW = 'EXTR' ZFLIANEW 'NBZO' ;
            ZFPNW = 0. ; XVITW = 0. ;
            'REPETER' BZLIAW NZNEW ;
*             un point support par zone - 2010 kich
              ZFPNW = ZFPNW +
                  ('EXCO' ('EXTR' ZFLIANEW 'FLIA' &BZLIAW 1 1) MLDUAL_LIA 'NOID') ;
              XVITW = XVITW +
                  ('EXCO' ('EXTR' ZFLIANEW 'FLIA' &BZLIAW 1 1) MVPRIM_LIA MLPRIM_LIA 'NOID') ;
            'FIN' BZLIAW ;
            ZFP = ZFP +  ZFPNW ;
          'FINS' ;
        'FINS' ;
      'FINS' ;
 
      'SI' (NZLIA > 0) ;
        XVIT2 = 0.D0;
        'REPETER' BZLIA NZLIA ;
          XVIT2 = XVIT2 +
                  ('EXCO' ('EXTR' ZFLIA 'FLIA' &BZLIA 1 1) MVPRIM_LIA MLPRIM_LIA 'NOID') ;
        'FIN' BZLIA ;
        'SI' ('EGA' ('TYPE' XVIT2) 'CHPOINT  ') ;
          MAAUT2 = 'DIFF' ('EXTR' VITI 'MAIL') ('EXTR' XVIT2 'MAIL') ;
          VITI = ('REDU' VITI MAAUT2) + XVIT2 ; 'DETR' XVIT2 ;
          'SI' ('EGA' ('TYPE' XVITW) 'CHPOINT  ') ;
            MAAUT3 = 'DIFF' ('EXTR' VITI 'MAIL') ('EXTR' XVITW 'MAIL');
            VITI = ('REDU' VITI MAAUT3) + XVITW ; 'DETR' XVITW ;
          'FINS' ;
        'FINS' ;
        'DETR' FLIAI ;'DETR' ZFLIA ;
      'FINS' ;
    'FINS' ;
*DMODI_NB
    ZVITET = VITI - conti. 'VITESSES';
    ACCEI  = UNSURH * 2. * ZVITET - conti.'ACCELERATIONS';
    ACCEI  = 'ENLE' ACCEI 'LX' ;
*FMODI_NB
  'FINS' ;
*
  'SI' WTAB.'NVSTNL' ; 
    estim.'VITESSES_FLUIDE_0' = estim.'VITESSES_FLUIDE' ;
    estim.'VITESSES_FLUIDE' = ZVIFL  ;
    'SI' ('NON' ('EXISTE' WTAB 'QNSL')) ;
      WTAB.'QNSL' = 'MASSE' 'PRES' ZMODL WTAB.'MAT_NSL' ;
    'FINSI' ;
    QNSL = WTAB.'QNSL' ;
    FPREFL = 'PROI' ZMODL ZVIFL ;
    ZPREFL = 'RESOUT' QNSL FPREFL ;
    estim.'PRESSION_FLUIDE' = ZPREFL ;
    WTAB.'CONV' = ZFLCONV            ; 
    'SI' ('NON' ZFLCONV) ; 'MESS' 'Pas de convergence NVST NLIN'; 'FINSI' ;
    'QUITTER' BONOCONV ;
  'FINSI' ;
*
* Redefinir les variables associees a l'instant TDEBUT, pour soit :
* -> ZICONV = VRAI : initialiser et preparer le pas suivant
* -> ZICONV = FAUX : redemarrer le pas courant en convergence forcee
  'SI' WTAB.'FOR_MECA' ;
*
    augmult = augmult * 0.55 ;
    si (augmult < 1d-3); augmult = 1d-3; finsi;
*
*   Conserver la partie Forces nonlineaires trouvee pendant ce pas pour 
*   s'en servir pour l'initialisation du pas suivant
*   DFNL = K*DU  - DF  - residu
*   Ktot*dutot ne contient pas les forces de reactions mais les forces
*   internes dues a un champ de deformation initiales ( thermique)
*   pour les forces suiveuses on fait delta FP
*   pour etre plus precis on fait  aussi intervenir le residu
*
    XXX1 = ZRAID * ZDEPT;
    XXX3 = (COEPI '*' ZDFINI) '+' ZFSUIV ;
    XXX4 = XXX1 - XXX3;
    XXX3 = XXX4 + RESIDC;
    XXX1 = XXX3 - FREAP ;
    XXX5 = XXX1 'ENLEVER' 'FLX ' ;
    'DETR' XXX1 ; 'DETR' XXX3 ; 'DETR' XXX4 ;
* en cas de non convergence on cumule les forces non lineaires
    'SI' WTAB.'CONV';
      ZFNONL = XXX5 ;
    'SINON';
      ZFNONL = XXX5 + ZFNONL ;
      'DETR' XXX5 ;
    'FINS';
*
*   Garder la derniere matrice KTAN calculee dans ETIQ si necessaire
    'SI' IKTAN ;
      'SI' IFEFP ;
*       'MESS' 'FEFP: Last KTAN is kept for next increment' ;
        ZLASTKTAN = ZRIKTA ;
      'SINON' ;
        'SI' (IKT_SAUV 'ET' IKT) ;
*         'MESS' 'KTAN : La matrice est conservee pour le pas suivant' ;
          ZLASTKTAN = ZRIKTA ;
        'FINS' ;
      'FINS' ;
    'FINS' ;
*
    'DETR' ZFEXT ; 'DETR' FREAP ;
*
    'SI' (IGRD 'OU' IFEFPUL) ;
      'SI' (knoconv '>' 1) ;
        'DETR' GEOM1 ;
      'FINS';
    'FINS';
*
*   Utile pour les calculs d'usure
    'SI' WTAB.'CAFROTTE' ;
      ZFLX = 'EXCO' (ZDEPT '*' ZCLIM) 'FLX' 'NOID' 'FLX' ;
      WTAB.'POST_COFR'.'RIGI_UNILA' = CRR ;
      WTAB.'POST_COFR'.'MULT_SOLUT' = 'EXCO' ZDEPT 'LX' 'NOID' 'LX' ;
      WTAB.'POST_COFR'.'GLISSEMENT' = ZFLX ;
    'FINSI' ;
*
    ZXDENO  = XDENO;
    ZXDENOM = XDENOM;
*
*   Faut-il partir en convergence forcee ou ok?
    WTAB.'CONV' = ZICONV ;
    'SI' ('NON' ZICONV) ;
*
*     -----------------------------------------------------------------
*     Preparation des donnees pour le nouveau sous-pas
*
*     Recalcul du pas a partir de TDEBUT + ZDT
      TDEBUT = TDEBUT '+' ZDT ;
      PASDT  = TI '-' TDEBUT ;
      'SI' (('&lt;EG' PASDT XPETIT) 'ET' ('OU' IVISCO IVIDOM IVIEXT IDYN ISOL)) ;
        'ERRE' 1155 ;
      'FINSI' ;
*
      'SI' (IGRD 'OU' IFEFPUL) ;
        GEOM1 = GEOM2;
        'FORM' GEOM2 ;
      'FINS';
*
      ZDEP0 = ZDETOT ;
      ZDEF0 = ZDEFF ;
      ZSIG0 = ZSIGF ;
      'SI' IPLAVI ;
        ZDEI0 = ZDEIF ;
        ZVAR0 = ZVARF ;
        'SI' ('NEG' WTAB.'MOVA' 'RIEN') ;
          ACC0 = 'EXCO' ZVAR0 WTAB.'MOVA' ;
        'FINS' ;
      'FINS';
      ZGRDU0 = ZGRDUF;
*  
      'SI' ITHER ;
        ZTET1  = ZTET2 ;
        ZTEMP1 = ZTEMP2 ;
        ETT0   = 'REDU' ETT ZMODL ;
        'SI' WTAB.'POR1' ;
          MSRTT0 = 'REDU' MSRTT ZMODL ;
        'FINSI' ;
      'FINSI' ;
*
      'SI' LOGDEF;
        ZDEFOR1 = ZDEFOR2 ;
      'FINS';
*     -----------------------------------------------------------------
    'FINS' ;
*
    'SI' (ZICONV 'OU' ('NON' WTAB.'CONVERGENCE_FORCEE')) ;
      'QUITTER' BONOCONV;
    'FINS';
*
*   Convergence forcee -> on reinitialise ZCLIM
    ZCLIM = ZCLIM0 ;
  'FINSI' ;
*FOR_MECA
*
* Nombre maximum de sous-pas atteint? si oui, arret de pasapas
  WTAB.'ISOUSPAS' = WTAB.'ISOUSPAS' + 1;
  'SI' (WTAB.'ISOUSPAS' >EG WTAB.'MAXSOUSPAS');
    'ERREUR' 996 ;
  'FINS';
*
'FIN' BONOCONV;
*----------------------------------------------------------------------*
*           Fin de la boucle de non convergence BONOCONV               *
*----------------------------------------------------------------------*
*######################################################################*
*----------------------------------------------------------------------*
*                     Preparation du pas suivant                       *
*----------------------------------------------------------------------*
WTAB.'KNOCONV' = KNOCONV;
*
'SI' WTAB.'FOR_MECA' ;
*
  'SI' ('NON' IFEFP) ; 'FORM' GEOREF0; 'FINSI' ;
  estim.'DEPLACEMENTS' = ZDETOT ;
  estim.'REACTIONS'    = 'REAC' estim.'DEPLACEMENTS' ZCLIM ;
  estim.'CONTRAINTES'  = 'REDU' ZSIGF ZMODLI;
  estim.'DEFORMATIONS' = 'REDU' ZDEFF ZMODLI ;
*
  'SI'  IPLAVI ;
    estim.'VARIABLES_INTERNES'        = 'REDU' ZVARF ZMODLI ;
    estim.'DEFORMATIONS_INELASTIQUES' = 'REDU' ZDEIF ZMODLI ;
  'FINS';
*
  'SI'  IDYN  ;
    'SI' ('NEG' WTAB.'REAPREC' 'INCONNU');
      maireac = 'EXTR' estim.'REACTIONS' 'MAIL' ;
      'SI' ('EGA' ('TYPE' maireac) 'MAILLAGE');
        reacdif = 'REDU' WTAB.'REAPREC' maireac ;
        estim.'REACTIONS' = estim.'REACTIONS' - REACDIF;
      'FINSI' ;
    'FINSI';
    estim.'VITESSES' = VITI;
*
    'SI' IMPLP; 
*     Correction des vitesses pour avoir des vitesses relatives nulles
*     aux points qui sont en contact
      estim.'VITESSES' = estim.'VITESSES' + VADD ;
    'SINON' ;
*     Appuis unilateraux + choc elastique : essai de corriger les 
*     vitesses fournies par le schema
      'SI' ('NEG' WTAB.'ZRAIDV' 'INCONNU') ;
        estim.'VITESSES' = VITEUNIL WTAB.'ZRAIDV' WTAB.'MASSE' 
                           estim.'VITESSES' ZDEPT conti.'DEPLACEMENTS' 
                           ZDT ZSDMBR WTAB ;
        'SI' ('EXIS' WTAB 'RATE_VITEUNIL');
          'ERREUR' 'Probleme dans viteunil' ; 
        'FINSI';
      'FINSI';
    'FINS';
*
    estim.'ACCELERATIONS' = ACCEI ;
    WTAB.'FOPL'           = ZFP ;
  'FINS' ;
*
  'SI' LOGPIL ; 
    PRECED.'COEFFICIENT_DE_PILOTAGE' = PRECED.'COEFFICIENT_DE_PILOTAGE' 
          'ET' ('PROG'ETA) ;                                                                         
  'FINSI' ;
*----------------------------------------------------------------------
*
  WTAB.'CLIM' = ZCLIM ;
  'SI' (IGRD 'OU' IFEFPUL) ;
    WTAB.'FOR' = GEOM2 ;
  'FINS';
  WTAB.'MAT1' = ZMAT22 ;
*
* Pour initialisation a partir du pas precedent
  WTAB.'DTPREC'    = WTAB.'DT';
  WTAB.'FNONL'     = ZFNONL ;
  WTAB.'INCREMENT' = ZINCREMENT ;
  WTAB.'RESIDU'    = RESIDC ;
*
* Contact-frottement
  'SI' ('NEG' ZRIBLO_M 'INCONNU') ;
    WTAB.'RIBLO_M' = ZRIBLO_M ;
    WTAB.'LISEA_M' = ZLISEA_M ;
  'FINSI' ;
*
* Pilotage
  'SI' IPILOT ;
    WTAB.'AUTOCOEF' = COEPI ;
  'FINSI' ;
*
  WTAB.'XDENO'             = ZXDENO ;
  WTAB.'XDENOM'            = ZXDENOM ;
  WTAB.'LASTKTAN'          = ZLASTKTAN ;
  WTAB.'NOMBRE_ITERATIONS' = ZNITE ;
*
  'SI' ITHER ;
    WTAB.'ETHER2' = 'REDU' ETT ZMODLI ;
    'SI' WTAB.'POR1' ;
      WTAB.'MSRTHER2' = 'REDU' MSRTT ZMODLI;
    'FINSI' ;
  'FINSI' ;
*
* Mise a jour du MODELE du PAS precedent
  WTAB.'MO_TOT_PREC' = ZMODLI ;
*
  WTAB.'DPSMAX' = DPSMAX ;
*
'FINSI' ;
*FOR_MECA 
*
'OPTI' 'PARA' FAUX;
*
IERRMEC = ('NON' WTAB.'CONV') 'ET' ('NON' WTAB.'CONVERGENCE_FORCEE') ;
*
'FINPROC' IERRMEC ;