Sources Esope | Cast3M

cflu3d
C CFLU3D    SOURCE    FD218221  24/02/07    21:15:03     11834          
c     A.Sellier sept. 2010 11:55:47 CEST 
c     revision  sept. 2022        
      SUBROUTINE cflu3d(WRK52,WRK53,WRK54,MWRKXE,WR14,
     # nbnb,idimb,teta13d,teta23d,nvarib,nstrsb,ifourb,dtb,trefb)
 
C     WRK52,53,54 segments déclarés dans le common DECHE
 
C     XMAT(NCOMAT)  =  COMPOSANTES DE MATERIAU
C     IVAL(NCOMAT)  =  INDICE DES COMPOSANTES DE MATERIAU
C     NCOMAT        = NOMBRE DE COMPOSANTES DE MATERIAU
C     XCAR(ICARA)   =  CARACTERISTIQUES
C     MFR           =  NUMERO DE LA FORMULATION DE L'ELEMENT FINI
C                   = 1       MASSIF
C                   = 3       COQUE MINCE ( COQ2 , COQ3  ET DKT )
C                   = 5       COQUE EPAISSE ( COQ6 , COQ8 )
C                   = 7       POUTRE
C                   = 9       COQUE MINCE AVEC CISAILLEMENT TRANSVERSE  ( COQ4 ET DST )
C                   = 11      LIQUIDE
C                   = 13      TUYAU
C                   = 15      LINESPRING
C                   = 17      TUYAU FISSURE
C                   = 19      RACCORD MASSIF
C                   = 21      RACCORD COQUE
C                   = 23      SURFACE LIBRE
C                   = 25      MEMBRANE
C                   = 27      UNIAXIALE
C                   = 29      THERMIQUE
C                   = 31      INCOMPRESSIBLES
C                   = 33      POREUX
C                   = 35      JOINT
C                   = 37      HOMOGENEISE
C                   = 39      TUYO
C                   = 41      TUYAU ACOUSTIQUE PURE
C                   = 43      RACCORD TUYAU FLUIDE
c*    DDAUX = MATRICE DE HOOKE ELASTIQUE
c*    NSTRS = NBRE DE COMPOSANTES DES DEFORMATIONS
c*    CMATE = NOM DU MATERIAU
c*    VALMAT= TABLEAU DE CARACTERISTIQUES DU MATERIAU
c*    VALCAR= TABLEAU DE CARACTERISTIQUES GEOMETRIQUES
c*    N2EL  = NBRE D ELEMENTS DANS SEGMENT DE HOOKE
c*    N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE
c*    MFR   = NUMERO DE LA FORMULATION
c*    IFOU  = type de formulation
c*    IB    = NUMERO DE L ELEMENT COURANT
c*    IGAU  = NUMERO DU POINT COURANT
c*    EPAIST= EPAISSEUR
c*    NBPGAU= NBRE DE POINTS DE GAUSS
c*    MELE  = NUMERO DE L ELEMENT FINI
c*    NPINT = NBRE DE POINTS D INTEGRATION
c*    NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES
c*    NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES
c*    SECT  = SECTION
c*    LHOOK = TAILLE DE LA MATRICE DE HOOKE
c*    TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI = TABLEAUX UTILISES
c*    POUR LE CALCUL DE LA MATRICE DE HOOKE
C-----------------------------------------------------------------------    
C      VARIABLES PASSEES PAR LES COMMONS COPTIO , ECOU  ET NECOU             
C                                                                           
C  IFOUR    INDICE DU TYPE DE PROBLEME
C            -3  DEFORMATIONS PLANES GENERALISEES
C            -2  CONTRAINTES PLANES                                         
C            -1  DEFORMATIONS PLANES                                       
C             0  AXISYMETRIQUE                                             
C             1  SERIE DE FOURIER                                          
C             2  TRIDIMENSIONNEL                                      
C  ITYP      TYPE DE FORMULATION MECANIQUE                              
C  --------------- ATTENTION ---------------
C     IL EST ACTIF APRES L APPEL DE VISAVI
C  -----------------------------------------
C       ITYP=1   CAS DES ELEMENTS MASSIFS                          
C       ITYP=2   CAS DES COQUES
C       ITYP=3   CAS DES MEMBRANES                               
C       ITYP=4   CAS DES CABLES ET DES BARRES                       
C       ITYP=5   CAS QUELCONQUE                                          
C       ITYP=6   CAS DES CONTRAINTES PLANES                      
C       ITYP=7   CAS DES COQUES A NU=0. OU CONTRAINTES PLANES        
C       ITYP=8   CAS DES MEMBRANES A NU=0. OU CONTRAINTES PLANES   
C       ITYP=9   CAS DES COQUES EPAISSES                           
C       ITYP=10  CAS DES JOINTS                                          
C       ITYP=11  CAS DES POUTRES                                            
C       ITYP=12  CAS DES TUYAUX                                              
C       ITYP=13  CAS DES COQUES AVEC CISAILLEMENT TRANSVERSE    
C            
C     ISTEP      flag utilise pour separer les etapes dans un calcul non local
C                ISTEP=0 -----> calcul local
C                ISTEP=1 -----> calcul non local etape 1 on calcule les seuils
C                ISTEP=2 -----> calcul non local etape 2 on continue le calcul
C                               a partir des seuils moyennes
C                ISTEP=3 ----->sous iteration de Helmholtz
C                ISTEP est défini par DEFINI.ESO      
C                                                                 
C     EPIN0      déformations inélastiques initiales
C     EPINF      déformations inélastiques finales
C     epst0      déformation totale initiale
C     EPSTF      déformation totale finale
C
C-----------------------------------------------------------------------    
C  SORTIES                                                                   
C     SIGF(NSTRS)   = CONTRAINTES FINALES                                    
C     VARF(NVARI)   = VARIABLES INTERNES FINALES
C     DEFP          = DEFORMATIONS PLASTIQUES                               
C     KERRE         = 0    TOUT OK
C-----------------------------------------------------------------------
 
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
c     ajout sellier 2022
-INC PPARAM
c     fin ajout sellier 2022
-INC CCOPTIO
-INC DECHE
*
      SEGMENT IECOU
*      COMMON/IECOU/NYOG,NYNU,NYALFA,NYSMAX,NYN,NYM,NYKK,
       INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7,
C      INTEGER NYOG, NYNU, NYALFA,NYSMAX,NYN, NYM,  NYKK,
     1  icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16,
C    .  NYALF1,NYBET1,NYR,  NYA,   NYRHO,NSIGY,  NNKX,  NYKX,  IND,
     2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,
C    .  NSOM,  NINV, NINCMA,NCOMP, JELEM, LEGAUS,INAT,  NCXMAT,
     3 icow25,icow26,icow27,icow28,icow29,icow30,ICARA,
C    .  LTRAC, MFR,  IELE,  NHRM,   NBNN, NBELEM,ICARA,
     4  icow32,icow33,NSTRS1,MFR1,icow36,icow37,icow38,
C    .  LW2,   NDEF,  NSTRSS,MFR1,  NBGMAT,NELMAT,MSOUPA,
     5  icow39,icow40,icow41,icow42,icow43,icow44
C    .  NUMAT1,LENDO, NBBB,  NNVARI,KERR1, MELEME
        INTEGER icow45,icow46,icow47,icow48,icow49,icow50,
     .          icow51,icow52,icow53,icow54,icow55,icow56
     .          icow57,icow58
      ENDSEGMENT
      SEGMENT XECOU
*      COMMON/XECOU/DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00
      REAL*8 xcow1, xcow2,xcow3,xcow4,DT,xcow6, xcow7
C     REAL*8 DTOPTI,TSOM, TCAR, DTT,  DT,TREFA,TEMP00
      ENDSEGMENT
 
c     segment de coordonnees des noeuds de l element
      integer nbnb,insb,idimb
      SEGMENT MWRKXE
        REAL*8 XE(3,nbnb)
      ENDSEGMENT  
 
*     
*     sellier/millard segment pour les Helmholtz 04 04 2020
      integer NBVIA
      SEGMENT WR14
        INTEGER INLVIA(NBVIA)
      ENDSEGMENT     
 
c-----------------------------------------------------------------------
c     nombre de noeuds max et table locale des coordonnees des 
c     noeuds de l element fini
 
c     ******************************************************************
c     NBNMAX3D
c     XE3D
c      include './nombre_noeuds.h'
-INC HNBRNEU
c     ******************************************************************
 
c     ********  taille du pseudo vecteur des contraintes ***************
c     tjrs 6 en raison de son utilisation dans as3d 
      integer nstrs3d 
      parameter (nstrs3d=6)       
      real*8 sig03d(nstrs3d),sigf3d(nstrs3d)
      real*8 epst03d(nstrs3d),epstf3d(nstrs3d)
      real*8 depst3d(nstrs3d)
      real*8 dtb,tref3d,alpha
 
c     variable logique : elasticite initialement isotrope
      logical iso1
 
c     nombre de parametres de tailles suivant la formulation
      integer ifou11,ntail3d
 
c     dimension de l espace de travail
      integer idimb3d 
c     temperatures debut et fin de pas , moyenne, pas de temps, volume rgi     
      real*8 teta13d,teta23d,temp3d,dt3d
c     variables de transfert de donnees ( a declarer suivant idvar4 et idvisc)
      integer ierr1,mfr11
c     variable pour passer le numero de l etape non locale      
      integer istep3d,nvarib,nstrsb,ifourb
 
c-----------------------------------------------------------------------
 
c     ******************************************************************
c     *      dimension des tableaux parametres materiau                *
c     ******************************************************************   
 
c     *****nombre de renforts et de parametres par renfort**************
c     NB_RENF
c     NB_PARA_PAR_RENF
c     NB_PARA_RENF
c     NB_VARI_PAR_RENF
c     NB_VARI_RENF
c      include './nombre_renforts.h'
-INC HNBRREN
c     ******************************************************************
 
c     ********Nombre de parametres pour variables non locales***********
c     NB_HELM: Nombre de variable d etat principales 
c     NB_PARA_PAR_HELM: Nombre de parametres par variable d etat 
c     NB_VARI_PAR_HELM: Nombre de variable internes par variables d etat 
c     NB_PARA_HELM: nombre de parametre total resultant 
c     NB_VARI_HELM: nombre de variables internes totales resultantes 
c     INLVIA3D(NB_HELM):copie du pointeur sur les numeros des vari traitee par helmholtz 
c     NBVIA3D: recuperation du nombre de Helmholtz actives  
c      include './nombre_helmholtz.h'
-INC HNBRHEL
c     ******************************************************************
 
c     parametres commun aux modele mecaniques elastiques isotropes
      integer NBELAS3D,NB_DECAL,NMAT3D
c     en isotrope le nombre de parametre de base est 4 = YOUN NU ALP RHO
      parameter (NBELAS3D=4)
c     decalage entre les parametres obligatoires (le modele)
c     et les parametres facultatifs (les parametres pour Helm holtz)
      parameter (NB_DECAL=3)
c     nombre de parametre ajoute par idmatr en fin de tableau des facultative
      parameter (NB_FIN=0)
 
c     ****************Nombre de parametres partie non lineaire de
c     la loi de comportement
c     NB_PARA_FLUENDO3D
c     NB_VARI_FLUENDO3D
c     NB_PARA_SUPP_FLUENDO3D
c     NB_VARI_SUPP_FLUENDO3D      
c      include './nombre_fluendo3d.h'
-INC HNBRF3D  
      integer NMAT0,NMAT1,NMAT2,NMAT3,NMAT4
c     nombre de parametres propre au modele sans les renforts ni le non local
      parameter (NMAT0=NBELAS3D)
      parameter (NMAT1=NMAT0+NB_PARA_FLUENDO3D)
      parameter (NMAT2=NMAT1+NB_PARA_SUPP_FLUENDO3D)
      parameter (NMAT3=NMAT2+NB_PARA_RENF+NB_DECAL)
      parameter (NMAT3D=NMAT3+NB_PARA_HELM+NB_FIN)
 
c     ******************************************************************
c     dimension maxi du tableau des parametres y compris non local   
      real*8 XMAT3D(NMAT3D)
c     ******************************************************************
c
 
c     ******************************************************************      
c     *         dimension des tableaux variables internes              *
c     ****************************************************************** 
 
c     **** nombre de variables internes loi de comportement non lineaire
 
 
c     nbr de variables pour la loi de comportement non lineaire
      integer NVAR0,NVAR1,NVAR2,NVAR3,NVARI3D
      parameter (NVAR0=0)
      parameter (NVAR1=NVAR0+NB_VARI_FLUENDO3D)
      parameter (NVAR2=NVAR1+NB_VARI_SUPP_FLUENDO3D)
      parameter (NVAR3=NVAR2+NB_VARI_RENF)
      parameter (NVARI3D=NVAR3+NB_VARI_HELM)
c     ******************************************************************
 
c     **** nombre total de variables internes **************************            
c     tableau local des variables internes      
      real*8 VAR03D(NVARI3D),VARF3D(NVARI3D)
c     ******************************************************************
 
c-----------------------------------------------------------------------
 
c***********************************************************************
c  remarque : parametres elastiques en formulation poreux massif (mfr=33)
c    en debut de xmat :
c   'YOUN'  :  module d'Young
c   'NU  '  :  coefficient de poisson
c   'RHO '  :  masse volumique (a la fin en non poreux)
c   'ALPH'  :  coefficient de dilatation thermique 
c    à la fin de xmat :
c   'MOB '  :  module de Biot
c   'COB '  :  coefficient de Biot
c   'PERM'  :  perméabilité intrinsèque
c   'VISC'  :  viscosité dynamique du fluide
c   'ALPM'  :  coefficient de couplage pression - température
c   attention L ORDRE N EST PAS LE MEME SI LE MATERIAU EST ORTHOTROPE 
c   cf deche.inc pour recuperer les noms exacts des variables
c   transfert des variables dans les tableaux du modele
c*********************************************************************** 
 
 
c*****debut du transfert des donnees vers les modeles (sellier 2022)****
 
c      ** on verifie si on a du non local via Helmholtz ****************
c      le decalage est RHO ALPH VISQ avant les parmetres non locaux si
c      il y a des non locales, alors qu'il est de 1 si on est en local
       if(WR14.EQ.0) then
         NBVIA = 0
       else
         NBVIA=INLVIA(/1)
c         print*,'NBVIA = ',NBVIA
c         do i=1,NBVIA
c            print*, 'I' ,i, 'INLVIA ' ,INLVIA(i)
c         end do
       endif  
 
c     suivant la formulation
      if(mfr.ne.33) then
c       ****** formulation non poreux **********************************       
c       print*,'ds cfluendo3d istep=',istep
        if(istep.eq.0) then
c            ***** formulation locale***********************************   
             if ((NMAT3D-NB_PARA_HELM).ne.nmatt) then
                 print*,'pb dimension de xmat cfluendo3d  local'
                 print*,'NBELAS3D',NBELAS3D
                 print*,'NMAT0',NMAT0
                 print*,'NB_PARA_FLUENDO3D',NB_PARA_FLUENDO3D
                 print*,'NMAT1',NMAT1
                 print*,'NB_PARA_SUPP_FLUENDO3D',NB_PARA_SUPP_FLUENDO3D
                 print*,'NMAT2',NMAT2
                 print*,'NB_DECAL',NB_DECAL
                 print*,'NB_PARA_RENF',NB_PARA_RENF
                 print*,'NMAT3',NMAT3
                 print*,'NB_PARA_HELM',NB_PARA_HELM
                 print*,'NB_FIN',NB_FIN
                 print*,'NMAT3D-NB_PARA_HELM',NMAT3D-NB_PARA_HELM
                 print*,'NMATT',NMATT
                 do i=1,nmatt
                     print*,'xmat(',i,')=',xmat(i)
                 enddo
                 read*
                 kerre=1
                 return
             else
                do i=1,(NMAT3D-NB_PARA_HELM)
                    XMAT3D(i)=xmat(i)
                enddo
c               on met a zero les Helmholtz
                do i=(NMAT3D-NB_PARA_HELM+1), nmatt 
                    XMAT3D(i)=0.
                enddo
             end if         
        else
c            **** formulation non locale *******************************
c            ces parametres facultatifs sont definies dans IDMATR        
             if(NMAT3D.ne.nmatt) then
                 print*,'pb dimension de xmat cfluendo3d  local'
                 print*,'NBELAS3D',NBELAS3D
                 print*,'NMAT0',NMAT0
                 print*,'NB_PARA_FLUENDO3D',NB_PARA_FLUENDO3D
                 print*,'NMAT1',NMAT1
                 print*,'NB_PARA_SUPP_FLUENDO3D',NB_PARA_SUPP_FLUENDO3D
                 print*,'NMAT2',NMAT2
                 print*,'NB_DECAL',NB_DECAL
                 print*,'NB_PARA_RENF',NB_PARA_RENF
                 print*,'NMAT3',NMAT3
                 print*,'NB_PARA_HELM',NB_PARA_HELM
                 print*,'NB_FIN',NB_FIN
                 print*,'NMAT3D',NMAT3D
                 print*,'NMATT',NMATT
                 do i=1,nmatt
                     print*,'xmat(',i,')=',xmat(i)
                 enddo
                 read*
                 kerre=1
                 return
             else
                do i=1,nmatt
                    XMAT3D(i)=xmat(i)
                enddo
             end if 
        end if       
       else       
c        *****formulation poreux ****************************************     
         print*,'on est en poreux dans cfluendo3d'
         print*,'Pb affectation des caracteristiques dans cas3d'
         do i=1,nmatt
             print*,'xmat(',i,')=',xmat(i)
         end do
         print*,'voir idmatr et deche.inc pour parametres en poreux'
         kerre=1
         return
       end if              
 
c********** recuperation des variables internes debut de pas ***********       
       if(NVARI3D.ne.nvarib) then
             print*,'Pb dimension de var0 dans cas3d'
             print*,'Verifier dimension table des variables internes'
             print*,'et CMPATBLTCR avec idvar4 et as3d'
             read*
             kerre=1
             return
       end if
c      transfert vers le tableau de fluendo3d           
       do i=1,NVARI3D
             VAR03D(i)=var0(i)            
       end do
 
c*********** recuperation des contraintes totales debut de pas *********
       if(mfr.ne.33) then 
c           on est pas en poreux       
            if(nstrsb.lt.nstrs3d) then
                do i=1,nstrsb
                    sig03d(i)=sig0(i)
                    sigf3d(i)=sigf(i)
c                   ATTENTION les depst 3-6 doivent etre des gamas                
                    depst3d(i)=depst(i)
                    epst03d(i)=epst0(i)
                    epstf3d(i)=epstf(i)
                end do
                do i=nstrsb+1,nstrs3d
                    sig03d(i)=0.d0
                    sigf3d(i)=0.d0
c                   ATTENTION les depst 3-6 doivent etre des gamas                 
                    depst3d(i)=0.d0
                    epst03d(i)=0.d0
                    epstf3d(i)=0.d0
                end do
            else
                do i=1,nstrs3d
                    sig03d(i)=sig0(i)
                    sigf3d(i)=sigf(i)
c                   ATTENTION les depst 3-6 doivent etre des gamas                 
                    depst3d(i)=depst(i)
                    epst03d(i)=epst0(i)
                    epstf3d(i)=epstf(i)
                end do
            end if
       else
            print*,'on est en poreux'
            print*,'Pb affectation des contraintes dans cfluendo3d'
            print*,'a terminer'
            do i=1,nstrsb
                print*,'sig(',i,')=',sig0(i)
            end do
            kerre=1
            return
       end if         
 
 
c**************   autres parametres a renseigner ***********************
 
c      initialisation indicateur d erreur       
       ierr1=0
c      indicateur isostropie elastique et de resistance
       iso1=.true.
c      numero de la formulation (33 pour poreux)
       mfr11=mfr 
c      type de formulation
       ifour11=ifourb         
c      pas de temps
       dt3d=dtb     
c      numero  pour le traitement non local eventuel       
       istep3d=istep
 
c************* recuperation des coordonnees de noeuds ******************
c      boucle sur les noeuds
       if(nbnb.le.NBNMAX3D) then
           NBNB3D=nbnb
           if(idimb.le.3) then
              idimb3d=idimb
           else
              print*,'pb table de coord des neouds ds cfluendo3d'
              kerre=1
              return
           end if
       else
          print*,'Element avec + de noeuds que capacite de cfluendo3d'
          kerre=1
          return
       end if
c      chargement des coordonnees des noeuds        
       do insb=1,NBNB3D
          do j=1,idimb3d       
c             print*,'xe(',j,insb,')=',xe(j,insb)   
             XE3D(j,insb)=xe(j,insb)
          end do
          if (idimb.lt.3) then           
            do j=idimb+1,3
              XE3D(j,insb)=0.d0
            end do
          end if            
       end do
c      on double les neouds si le probleme est 2D       
       if ((ifourb.eq.0).or.(ifourb.eq.-1)) then
c        on double les noeuds avec dimension 3 recuperee ds xmat
c        toujours placé en NMAT1-1 cf def de namat1 plus haut + idvisc
         dimension3=XMAT3D(NMAT1-1)
         NBNB3D=2*nbnb
         do insb=nbnb+1,NBNB3D
            do j=1,2       
c               print*,'xe(',j,insb,')=',xe(j,insb)   
                XE3D(j,insb)=XE3D(j,(insb-nbnb))
            end do
            XE3D(3,insb)=dimension3
         end do
       end if
 
c      ****** recuperation de la temperature de reference **************
       tref3d=trefb
c       print*,'verif temperature de reference fluendo3d',trefb
 
c      * intialisation des variables de sortie  a leur valeur initiale * 
c      en cas de 1ere etape non locale
       if(istep.eq.1) then
             do i=1,nvari3d
                 varf3d(i)=var03d(i)
             end do
       else   
             do i=1,nvari3d
                 varf3d(i)=0.d0        
             end do
       endif
c****************** fin du transferts des donnees **********************
 
c*************recuperation des numeros des variables Helmholtz *********
 
c      transfert dans un tableau de fluendo3d dimensionne au max des
c      facultative helmholtz de IDMATR (actuellement 10)
       if(NBVIA.GT.NB_HELM) then
          print*,'Pb de dimensionnement avec les variables de Helmholtz'
          print*,'dans cfluendo3d:',NBVIA,' > ',NB_HELM
          print*,'Verifier NB_HELM et nbr de non locales'
          print*,'declarees nombre_helmholtz.h (NBRHEL.INC)'
       else
c         transfert vers les variables locales
          if (NBVIA.gt.0) then
             NBVIA3D=NBVIA
             do i=1,NB_HELM
                 if(i.le.NBVIA3D) then
                     INLVIA3D(i)=INLVIA(i)
                 else
                     INLVIA3D(i)=0
                 end if
             end do 
          else
             NBVIA3D=NBVIA
             do i=1,NB_HELM
                 INLVIA3D(i)=0
             end do
          end if
       end if
c************* fin de recuperation des numeros des variables Helmholtz *        
 
c************* appel au modele  fluendo3d*******************************  
 
      call fldo3d(XMAT3D,NMAT3D,sig03d,sigf3d,depst3d,nstrs3d,
     #VAR03D,VARF3D,NVARI3D,teta13d,teta23d,dt3d,ierr1,iso1,mfr11,
     #ifour11,istep3d,epst03d,epstf3d,tref3d,NBELAS3D,NB_PARA_FLUENDO3D,
     #NB_PARA_SUPP_FLUENDO3D,NB_VARI_FLUENDO3D,NB_VARI_SUPP_FLUENDO3D,
     #NB_DECAL,NBNMAX3D,NBNB3D,IDIMB3D,XE3D,NBVIA3D,INLVIA3D,
     #NMAT0,NMAT1,NMAT2,NMAT3,NVAR0,NVAR1,NVAR2,NVAR3,NB_PARA_FIBRE_U)
 
 
c************ re-affectation des tableaux de variables internes ********
       do i=1,nvarib
            varf(i)=VARF3D(i)
*           print*,'cas3d f(',i,')=',varf(i),' i',var0(i)
       end do
c       read*
c      ecriture des contraintes totales fin de pas       
       do i=1,nstrsb
          sigf(i)=sigf3d(i)
       end do        
 
c***********************************************************************
c       traitement de l erreur d ecoulement
        if(ierr1.eq.0)then
            kerre=0 
        else
            kerre=1
        end if
      return
      end
 
c***********************************************************************