c   Programme dme.f - UNIX
c**********************************************************************
      subroutine dme(dmedat,nfm,npo,dnanfp,iue,        iuj,
c                ---   etiq            in  e(q ou p)in  in 
     1               lineref,sundec, ideld)
c                     in      in      out
c vmoy lu  dans fichier l, dam/s, >0 blueshift
c le jname sert pour
c     lire le fichier j  si lcor.lt.0  (-lcor est le tableau de reference)
c     creer le fichier j si lcor=0 (recouvrements a partir de pixstep)
c                        ou lcor.gt.0 (correl a partir du tableau lcor)
c le dessin continue le q....ps
c            le plot de ond (2 derniers tableaux)
c            est toujours recale sur les "tabmoy" correspondants

c                              iud  iue     iuj
c les fichiers ouverts seront   51   52      53
c                             les d  q ou p   j

      dimension corm(150),kedge(150),cr(10),ond(1024,128,3)
      dimension xdes(2),ydes(2),poids(1024,128,3),tabmoy(30),sigma(30),
     1          sig(128,30,6),sigmoy(30,6)
      integer*4 head(512),ll(30),i1(150),jr(150),nrej(90),
     1          cmoy(100,2,30),dlb(30)
      integer*4 pixstep,pixstepp,corstep,fac100,dlbd
      integer icor(150),jcor(150),xcor(150),ycor(150)
      integer*2 lec(2048),tab(400,4000)

      character*6 dmedat
      character*10 jdate
      character*20 dname,jname,jnameref,ename,pname,dnanfp
      character*1 let1,let2
      character comm*80,data*6,end*3
      dimension dnanfp(2000,2)


      type *,'etape0'!!!!!!!!!
      call date(jdate)
      type *,'jname :',jname !!!!!!!!!!!
       type *,'dmedat : ',dmedat
c      jname=dname !!!!!!!!!!
c      write(jname(1:1),'(a)')'j' !!!!!!!!
c                  voir plus loin lineref
c                     pour le cas de plusieurs raies simeltanees
c                     reference correlations
    
c      
c attention 150 cliches maximum
c          128 lignes maximum par tableau
c         1024 pixels maximum par ligne
c         4000 lignes maximum dans l'image finale
c         4000 fichiers maximum
c
c        iminus = intervalles en i elimines en debut et fin    inutilise  
c        jminus = intervalles en j elimines en debut et fin    inutilise  
c        jlap2 = demi-intervalle pour reduction du recouvrement
c                     (cas des grands delta lambdas non couverts dans tout le
c                      champ)
c        lcor = indice du tableau servant de reference pour les correlations
c             si lcor = -1, on lit un nom de fichier "correlations" type  J
c                            de la raie lineref
c             si lcor = 0, on utilise pixstep sans calcul de correlations
c             si lcor.gt.0, on calcule pour chaque ny
c                         a partir des tableaux l=lcor,npo=nporef
c        jrsom # jrmoyen .le. 128 = recouvrement sur (1,jm) # jm-pixstep
c        ijcormax .le. 50
c        milcor = 1000 * correlation minimum admissible
c        corstep=pas pour les differences premieres en i dans calcul correl
c               =pas pour l'integration
c        jdel=   pas pour les differences premieres en j dans calcul correl
c              (pas 1 en j pour l'integration)
c     milsig  moyennes et
c             corrections calculees a partir pts < moy+/-milsig*sigma/1000
c       pas de calcul de corrections si moins de 3 images
c     les corrections sont lissees sur im/10
c     corrections des defauts d'interpolation: icorrec=1
c                                  les moy sont calculees sur les pts dont 
c                                        l'int est > suplec/4
c                                       (suplec=max 1-51/im-50,im)
c         (vitesses:   ll=-1 tab=vit-moy  en dam/s   ??????)
c         (protus: en dehors du disque, ajouter 2000, sauf aux points 0 ??)
c      icorrec=0   pas de correction des defauts d'interpolation
c              1   correction avec repliement en fonction de trj
c              2   correction sans repliement

c     iplot =1  plots de controle pour 2 dernieres cartes
c                   
c---------     
c
      pi=3.1415927
      nm=nfm

      end='end'
1     format(a)
2     format(1x,a)
3     format(12i6)
4     format(4i6,f6.3)
5     format(i7)
6     format(2(5x,a))
7     format(10i6)
9     format(' ltrj=',i6)


      type *,'dans dme step1' !!!!!!!!!!!!!!!!!
      rewind(8)
      type *,'dmedat(1:6)',dmedat(1:6)
      call SNAME(8,dmedat)
      type *,'apres sname'
      read(8,1)comm
      write(3,1)comm
      read(8,7)ll
      write(3,7)ll
c               1 pour chaque tableau a traiter du fichier d
c               0 pour les tableaux non traites
      read(8,1)comm
      write(3,1)comm
      read(8,3)iminus,jminus,jlap2,fac100,lcor,ijcormax,corstep,
     1         milcor,pixstepp,maxdecal,jdel
      write(3,3)iminus,jminus,jlap2,fac100,lcor,ijcormax,corstep,
     1         milcor,pixstepp,maxdecal,jdel
      if(jdel.eq.0)jdel=corstep

      type *,'step 2'
c  Si pixstepp=0, on prend pour pixstep la valeur lue dans les
c                 headers des fichiers b
c                                  a programmer ensuite!!!
      type *,'pixstepp',pixstepp
        if(pixstepp.ne.0)pixstep=pixstepp
      type *,'corstep',corstep
      if(corstep.eq.0)corstep=1
      read(8,1)comm
      write(3,1)comm
      read(8,3)lbda,dlbd,icorrec,milsig,mintab,
     1         iplot,jsor,ilissond,nporef,ideld
      write(3,3)lbda,dlbd,icorrec,milsig,mintab,
c                          0/1/2
     1         iplot,jsor,ilissond,nporef,ideld
c                          0/1  0/1   0/1        ~im/20  1/2

        if(iplot.ne.0)call pgadvance

      type *,'step 3'
      read(8,1)comm
      write(3,1)comm
      read(8,3)iz1,iz2,iz3,iz4
      write(3,3)iz1,iz2,iz3,iz4
      type *,'step4' 
      read(8,1)comm
      write(3,1)comm  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
      read(8,6)let1,let2
      write(3,6)let1,let2
c        if(jrsom.eq.0)jrsom=40
c        if(jrsom.gt.128)jrsom=128
c        if(ijcormax.eq.0)ijcormax=20
c        if(ijcormax.gt.50)ijcormax=50
c
c            x=v*convervx  (x in units ltrj)
c--------------------------------------
c                            lecture du header 1er fichier d
        iud=51
      dname=dnanfp(1,npo)
      call openold20(dname,sundec,51)
c          ---------                 
       type *,'dname:',dname

      jname=dname
c     -----
      write(jname(1:1),'(a)')'j'
      jnameref=jname


      rewind(iud)
      read(iud)head
                     type*,' dme:head de d ',(head(k),k=1,20)
      im=head(2)
      jm=head(3)    
      lm=head(4)
      npo=head(5)
      nxx=1
      nyy=1

c      day=head(6)
c      time=head(7)
c      iline=head(8)
c      lbda=head(9)

      ncm=head(11)
      ltrj=head(12)
          if(icorrec.eq.2)ltrj=10000
         type*,' dme: head(12) ',ltrj
c                             ceci remplace les moyennes par periode
c                                       par les moyennes sur tout le champ
c       convervx=-float(lbda)*float(ltrj)/(3000000.*float(dlbd))

       do l=1,lm
       dlb(l)=head(l+12)
       enddo
c      center=head(lm+10)
c      sum=head(lm+11)
c      lmp=head(lm+12)....
c      lbdl=head(lm+13)....
c      lbdstep=head(lm+14)....
c                les dimensions du fichier d sont:  IJ l po sto XY
c                        "           "     e seront:UV l po sto  Y
c les boucles qui suivent respectent l'ordre de creation de e:
c                                           Y  sto po  l   X  J   I
c                                          150    140 100 80 50  3interval.
      type*,' dme: lecture header d ',(head(i),i=1,lm+12) 
     
c--------------------------------------------
c                 calcul des correlations si lcor.gt.0    sur npo
      if(lcor.gt.0)then

c       if(i1i2.eq.0)then
      type *,'avant mulcor' !!!!!!!!!!!!!!!
       call mulcor
c           ------
     1  (nfm,npo,dnanfp,sundec,pixstep,lcor,jrsom,ijcormax,ijc,corstep,
     2  im,jm,nm,i1,jr,iminus,jminus,i1i2,icor,jcor,corm,milcor,
     3     kedge,xcor,ycor,maxdecal,jdel)
c                        nm est reduit au dernier de correl > milcor/1000
c       endif
c
       type *,'apres mulcor' !!!!!!!!!!!!!!!!!!!!!
      is=i1i2+1
c
      js=jm-2*jminus
        do n=2,nm
        js=js+jm-jr(n)
        enddo
c
      jsp=jm-2*jminus+(nm-1)*pixstep
215   format(' js,jsp=',2i6)
      write(3,215)js,jsp
c      if(js.ne.jsp)jr(nm)=jr(nm)+js-jsp
c      js=jsp
c
c      nf2=nf1+nm-1
        type *,' Fin de mulcor'
      goto18

      endif

c--------------------------------------
c pas pixstep si lcor=0          calculer a partir de im,jm,pixstep
      if(lcor.eq.0)then
            is=im-2*iminus
            js=jm-2*jminus+(nm-1)*pixstep
          do n=1,nm
          corm(n)=0
          enddo
        jr(1)=0
        i1(1)=1
          do n=2,nm
          jr(n)=jm-pixstep
          i1(n)=1
          enddo
                   write(3,*)' dme: jr ',(jr(n),n=1,nm)
         goto 21
      endif
c--------------------------------------------
c calcul precedent de correlations lu dans fichier "j" si lcor=-1
      if(lcor.eq.-1)then
       jnameref=dnanfp(1,nporef)
       write(jnameref(1:1),'(a)')'j'
c       write(jnameref(20:20),'(i1)')lineref
c                                              a modifier quand on connaitra
c                                              les codes de caracteres
       call LECJ(jnameref,sundec,is,js,nm,jminus,jm,vmoy,
c           ----
     1          i1,jr,corm)
      goto 21
      endif
c**********************************************
c  On dispose maintenant de 
c          iminus, i1(...nm), i1i2
c          jminus, jlap2
c  qui decrivent les jonctions du balayage nxx.
18    continue
      if(iplot.ne.0)
     1         call sdescor(ijc,icor,jcor,corm,nm,xcor,ycor)
c                   -------
21    continue
c------
      if(icorrec.eq.3)then
      call scorond3(nfm,npo,dnanfp,sundec,jm,im,lm,lec,kpiv,sig,sigmoy,
c          -------
     1   iminus,jminus,milsig,ll,ilissond)

      endif
c-------
          do i=1,is
          do j=1,js

          tab(i,j)=32000
          enddo
          enddo
c
c            cliche 1   JA = JB
c                        .
c                       JC
c                          .                   (commun)  JA du 2
c                          +JRecouvrement      (commun)  JD du 1
c            cliche 2   JB
c                        .
c                       JC
c                          .                   (commun)  JA du nm
c                          +JR                 (commun)  JD du 2
c            cliche nm  JB
c                        .
c                       JC = JD
c
c
22    format(' is=' i5,'  js=' i5)
      write(3,22)is,js
        ldeb=0
c                       suivent boucle tableaux l
c                                 boucle X (nf) - ouvertures d
c                                   boucle ltz (recherche tableau l dans nf)


      type *,'avant boucle sur les tableaux : is,ij',is,js !!!!!!!!!!!!
c-------------------------
c                                                  grande boucle ltz
      do100 ltz=1,30
            ides=3
            if(ltz.ge.lm-1) ides=ltz-lm+2    
      if(ll(ltz).eq.0)goto100
      type *,' ltz ',ltz
      ldeb=ldeb+1
      if(ldeb.gt.1.and.ltz.gt.lm)goto100
c
      type *,' nm=',nm
c      if(iplot.ne.0)then
c        if(iplot.gt.1)then
c          ides=ides+1
c        else
c          if(nf1.eq.1)ides=ides+1
c        endif
c      endif
c
c      if(ides.gt.3)ides=3
c
225   format(/' ltz,milsig,ides=',3i6)
      write(3,225)ltz,milsig,ides
        lpiv=ll(ltz)
        kpiv=ltrj
c------------------                 calcul    corrections ondulations

c                                 icorrec=0   pas de correction
c                                         1   moyennes repliees
c                                         2   moyennes non repliees
c                                         3   histogrammes non replies
c                                                 deja fait
        if(icorrec.eq.0)goto228

c------           cas icorrec=1 ou 2
      if(icorrec.eq.1.or.icorrec.eq.2)then

      if(icorrec.eq.2)kpiv=100*(jm+1)
      call scorond(nfm,npo,dnanfp,sundec,jm,im,lec,ltz,kpiv,ond,poids,
c          -------
     1   iminus,jminus,milsig,ides,tabmoy,  sigma,lpiv,ilissond)!!!!!!!!!!
c                  tabmoy: pour normaliser les dessins d'intensites
      type *,' Fin scorond'
c------                                dessin ond et poids
c vmoy pour ltz=2
c

      if((ides.eq.1).and.(iplot.ne.0))then
      piv=float(js)
        call pgvport(0.1,0.9,0.1,0.8)
        call pgsch(1.5)
        call pglabel('correl /   x / y      /        poids ', 
     1              ' ',dname//'     '//jdate)
      call pgvport(0.7,0.9,0.1,0.8)
      call pgwindow(0.,1.,1.,piv)
      call pgbox('bcnts',0.5,5,'bcmts',0.,0)
      call pgsch(1.)
      endif
c
      endif
c----
228          jta=1

        write(3,*)'avant boucle X, jta=',jta 
c -----------------------------------------  boucle nfm
        do80 ncli=1,nfm
c                   ---
        iud=51
        dname=dnanfp(ncli,npo)
            type*,dname
        call openold20(dname,sundec,51)
c            ---------
      rewind(iud)
      read(iud)head

        if(ncli.ne.1)jta=jta+pixstep
c                        -----------
        if(ldeb.eq.1)then
c        type 1,dname
c        write(3,1)dname
        endif

      lmax=0
      do llt=1,30
      if(ll(llt).ne.0.and.llt.le.lm)lmax=lmax+1
      enddo
      if(lm.lt.ltz)goto70
c
          lapa=jlap2
          lapb=jlap2
c          jr2a=jr(ncli)/2
c             if(ncli.ne.nfm) jr2b=jr(ncli+1)/2
c          if(jlap2.lt.0)then
c            lapa=jr2a-1
c            lapb=jr2b-1
c          endif
        ia=i1(ncli)
c             write(3,*)' dme: ncli,jr2a,jr2b ',ncli,jr2a,jr2b

        ja=1
        if(ncli.ne.1) ja= 1
        jb=ja
        if(ncli.ne.1) jb= 1+ jr(ncli)-lapa
        jd=jm
        if(ncli.ne.nm)jd= jm
        jc=jd 
        if(ncli.ne.nm)jc= jm-jr(ncli+1)+lapb
c       
      if(ldeb.eq.1)then      
23    format(' ncli,ja,jb,jc,jd,jta,head ='i4,4i3,i4,
     1         i2,3i4,i2,2i6)
      write(3,23)ncli,ja,jb,jc,jd, jta,(head(ih),ih=1,7)
c      type 23, ncli,ja,jb,jc,jd, jta
24    format(' jminus or jlap2 too large')
26    format(' jlap2 too large')
          if(ja.lt.jminus+1)then
           type *,' ja.lt.jminus+1'
          write(3,24)
          close(unit=iud)
          close(iue)
          close(3)
          close(8)
          call pgend
              stop
          endif

          if(jd.gt.jm-jminus)then
           type *,' jd.gt.jm-jminus'
          write(3,24)
          close(unit=iud)
          close(iue)
          close(3)
          close(8)
          call pgend
              stop
          endif


          endif

c          if(jc.lt.jb)then
c          write(3,26)
c          goto100
c          endif
  
c-----------------------------------------   ecriture header q ou p
        if(ldeb.eq.1.and.ncli.eq.1)then
        head(4)=lm
        head(2)=is
        head(3)=js
         ipiv=8+lmax
27    format(' head file q ou p ',10i6)
        write(3,27)(head(i),i=1,ipiv)
        write(iue)head
c             ---
        endif
c-----------------------------------------------------------
c                                                 boucle lm  de lecture
          do60 l=1,lm
c               type*,'    saut avant ltz, lecture tableau ltz de d'
            do50 j=1,jm
            read(iud)(lec(i),i=1,im)
c                ---
            if(l.ne.ltz)goto50
c--------
             ipiv=im/2
             jpiv=jm/2
        if(j.eq.jpiv)type*,' dme: ncli,l  lec(im/2,jm/2) ',
     1                            ncli,l, lec(ipiv)  

            
            jt=jta+j-ja
c-----------------------------------------------1er intervalle
            if(j.ge.ja.and.j.le.jb)then
            fac=1.
            if(jb.ne.ja)then
c              if(jlap2.lt.0)then
c              fac=float(fac100)/100.
c              else
              piv=pi*float(j-ja)/float(jb-ja)
              fac=(1.-cos(piv))/2.
c              endif
      if(ides.eq.1.and.iplot.ne.0)then
      yd=jt
      xd=1.-fac
      call pgpoint(1,xd,yd,2)
c          -------        
      xd=fac
      call pgpoint(1,xd,yd,2)
c          -------
      endif
c
            endif
c
              do30 it=1,is
              i=it+ia-1
                if(tab(it,jt).eq.32000.and.ncli.ne.1)then
                goto30
                endif
                if(lec(i).eq.32000)then
                 tab(it,jt)=32000
                 goto30
                endif
c
                piv1=tab(it,jt)
c                piv2=lec(i)-tabmoy(ltz)
                piv2=lec(i)
                  if(piv2.eq.32000.)then
                  tab(it,jt)=32000
                  goto30
                  endif                
           if (icorrec.eq.1.or.icorrec.eq.2)piv2=piv2-ond(i,j,ides)
           if(icorrec.eq.3)piv2=cormom(piv2,sig,sigmoy,j,l)
c                               ------
c
c                  corvit=convervx*ond(i,j,ides)
c                  jcorvit=j+corvit
c                  if(jcorvit.gt.jm)jcorvit=jcorvit-ltrj
c                  if(jcorvit.lt.1)jcorvit=jcorvit+ltrj
c                  piv2=lec(i)+ond(i,jcorvit,ides)
28    format(' ltz,jt,j,fac,ancien,nouveau,tab',
     1           3i5,      f6.3, 3f7.0, i7)
              tab(it,jt)=(1.-fac)*piv1+fac*piv2
c             ----------
      if(it.eq.is/2)then
      piv3=piv2-piv1
c      write(3,28)ltz,jt,j,fac,piv1,piv2,piv3,tab(it,jt)
      endif
30            continue
            endif
c----------------------------------------------2eme intervalle
            if(j.ge.jb.and.j.le.jc)then
      yd=jt
      xd=1.
      if(ides.eq.1.and.iplot.ne.0)call pgpoint(1,xd,yd,1)
c                                      -------
c
              do it=1,is
              i=it+ia-1
c
              if(lec(i).ne.32000)then
c               piv=lec(i)-tabmoy(ltz)
                piv=lec(i)
           if(icorrec.eq.1.or.icorrec.eq.2)piv=piv-ond(i,j,ides)
           if(icorrec.eq.3)piv=cormom(piv,sig,sigmoy,j,l)
c                              ------
c
c                if (imuladd.eq.0) then
c                  tab(it,jt)=lec(i)*ond(i,j,ides)
c                else
c                  corvit=convervx*ond(i,j,ides)
c                  jcorvit=j+corvit
c                  if(jcorvit.gt.jm)jcorvit=jcorvit-ltrj
c                  if(jcorvit.lt.1)jcorvit=jcorvit+ltrj
c                  tab(it,jt)=lec(i)+ond(i,jcorvit,ides)
 
             tab(it,jt)=piv
c            ----------
              else
                tab(it,jt)=32000
              endif
              enddo
            endif
c-------------------------------------------------3eme intervalle
            if(j.ge.jc.and.j.le.jd)then
      yd=jt
      xd=1.
      if(ides.eq.1.and.iplot.ne.0)call pgpoint(1,xd,yd,1)
c                                      -------
c
              do it=1,is
                i=it+ia-1
                if(lec(i).ne.32000)then
c                  piv=lec(i)-tabmoy(ltz)
                   piv=lec(i)
c                  corvit=convervx*ond(i,j,ides)
c                  jcorvit=j+corvit
c                  if(jcorvit.gt.jm)jcorvit=jcorvit-ltrj
c                  if(jcorvit.lt.1)jcorvit=jcorvit+ltrj
c                  tab(it,jt)=lec(i)+ond(i,jcorvit,ides)
            if(icorrec.eq.1.or.icorrec.eq.2)piv=piv-ond(i,j,ides)
            if(icorrec.eq.3)piv=cormom(piv,sig,sigmoy,j,l)
c                               ------
c

                  tab(it,jt)=piv
c                 ----------
                else
                  tab(it,jt)=32000
                endif
40            enddo
            endif
c                                   fin 3eme intervalle
50          continue
c                                        fin boucle j
60        continue
c-----------------------------------------    fin boucle lm de lecture 
70        continue

        close(unit=iud)
c       -----
80    continue
       type *,'fin boucle X'
c-------------------------------------------------------- fin boucle nfm
      if(iz1.ne.0)then
        if(iz4.eq.0)then
        iz4=is
        iz2=iz4/2
        iz3=iz2+1
        endif
      endif
c
      type*,' dme: is,js,ltz,tabmoy(ltz) ',is,js,ltz,tabmoy(ltz)

        do j=1,js
          do i=1,is
          lec(i)=tab(i,j)
            if(lec(i).eq.32000)then
             lec(i)=tabmoy(ltz)
            else
c              if(ll(ltz).eq.1)lec(i)=2000.*lec(i)/tabmoy(ltz)
c                                 intensites non normalisees
               if(ll(ltz).eq.-1)lec(i)=lec(i)-tabmoy(ltz)
c si ll=-1, vitesses recalees 
            endif
          enddo
c
      
        call zero3(iz1,iz2,iz3,iz4,0,0,0,lec,1,is,1,
c            ----- si 0, inefficace
     1             0,ar,br,cr,item,nrjm,nrej,piv1,piv2,ifin)
c               type *,j,ar,br,item,nrjm,piv1,piv2
c       sans effet si iz1=0
c
          if(ifin.ne.0)then
              do i=1,is
              tab(i,j)=0
              enddo
c          type 85
85        format(' Sortie en erreur dans zero3')
c          write(3,85)
c
           else
          do i=1,is
          tab(i,j)=lec(i)
c         --------
          enddo
           endif
c
        enddo
c
c      if(ltz.eq.2)call srejec(tab,is,js,3000,vmoy)
c                      ------
c
      type *,'is,js=',is,js
          do j=1,js
c            if(iz1.ne.0)then
c            do i=1,is
c            kpiv=tab(i,j)
c            if(kpiv.ne.32000)tab(i,j)=kpiv+ymoyij
c            enddo
c            endif

          write(iue)(tab(i,j),i=1,is)
c               -----                         ecriture sur q ou p
90    format(' ltz, tab='10i6)
            jtest=j-100*(j/100)
            if(jtest.eq.0)write(3,90)ltz,(tab(i,j),i=1,is,200)
          enddo
c
100   continue
      type *,'fin boucle ltz'
c-------------------------------                             fin boucle ltz

      type *,'nrepetmax : ',nrepetmax,'jm : ',jm,'ltrj : ',ltrj
      nrepetmax=100*jm/ltrj+1
      type *,'nrepetmax : ',nrepetmax !!!!!!!!!!!!!!!!!!!
      den=nm*nrepetmax
c-------
c                              dessin des 2 derniers tableaux calcules
       type *,'iplot : ',iplot 
      if(iplot.ne.0.and.(icorrec.eq.1.or.icorrec.eq.2))
     1    call sdesond(ond,poids,im,jm,den,tabmoy)
c              -------
c-------------------------------

c                                Creation du fichier j (vmoy et corr)
c        type *,'jname: ',jname 
c        call system('rm '//jname)
c        call opennew20sf(jname,sundec,54)!opennew20 -> opennew20sf seq format
c            ---------
107   format(' moyenne dopplershift           ')
108   format(f12.1)
109   format(a)
110   format(5i6)
111   format('    is    js    nm jminus   jm  ')
112   format(2i6,f6.3)
113   format('    i1    jr  corm              ')
114   format(3i4)
      write(iuj,109)jname
      write(3,109)jname
c      write(iuj,114)2,nyy,npol
c      write(3,114)2,nny,npol
c        endif

      write(iuj,107)
      write(3,107)
        type 107
c        vmoy(npo,ny)=tabmoy(2)
c        type 108,ny,npo,vmoy(npo,ny)
      write(iuj,108)tabmoy(2)
      write(3,108)tabmoy(2)
      write(iuj,111)
      write(iuj,110)is,js,nm,jminus,jm
c                                      attention nm
      write(3,111)
      write(3,110)is,js,nm,jminus,jm
      type 111
      type 110,is,js,nm,jminus,jm
      write(iuj,113)
      write(3,113)
      type 113
       do n=1,nm
       write(iuj,112)i1(n),jr(n),corm(n)
       write(3,112)i1(n),jr(n),corm(n)
       type 112, i1(n),jr(n),corm(n)
       enddo

c      close(unit=iuj)

      return
      end
c***************************************************************************
      subroutine mulcor
c                ------
     1       (nfm,npo,dnanfp,sundec,pixstep,lcor,jrsom,ijcormax,ijc,
     2        corstep,im,jm,nm,i1,jr,iminus,jminus,i1i2,
     3        icor,jcor,corm,milcor,kedge,xcor,ycor,maxdecal,jdel)
c
      integer*4 i1(150),jr(150),kedge(150),head(512),
     1          icor(150),jcor(150),pixstep,corstep,sundec
      integer*2 tab(400,100,150),         !!!   attention place memoire!!!!!
c                                  sans rapport avec le tab du principal
     1          tab1(1024,128),tab2(1024,128)
      dimension cor(101,101),corm(150),kcor(101),xcor(150),ycor(150)
      integer*2 lec(1024)

      character*20 dnanfp(2000,2)

c
c      if(lcor.le.0)return
c
1     format(a)
2     format(12i6)
c                                      boucle x      (nfm)
      iud=51

      do 20 nf=1,nfm
      call openold20(dnanfp(nf,npo),sundec,iud)
c          ---------      
      rewind(iud)
      read(iud)head
      im=head(2)
      jm=head(3)
      lm=head(4)
      npol=head(5)
c      nsto=head(6)
c                             !!!!!!!!!!!!! verifier
      nxx=1
      nyy=1
    
          jrsom=jm-pixstep
          if(jrsom.gt.128)jrsom=128
          ijc=jm-pixstep
          if(ijc.gt.ijcormax)ijc=ijcormax
5         format(15x,' im,jm,lm,jrsom,ijc ' 5i5)
c          type 5,im,jm,lm,jrsom,ijc
c          write(3,5)im,jm,lm,jrsom,ijc
  
      nm=nfm
      nmp=nfm-1
         icor(1)=ijc
         jcor(1)=ijc
         corm(1)=1.
     
      do l=1,lcor
        do j=1,jm
          if(l.eq.lcor)then
          read(iud)(tab(i,j,ncli),i=1,im)
c                   -------------
          else
          read(iud)(lec(i),i=1,im)
          endif
        enddo
      enddo 
        close(unit=iud)

20    continue       
c                        fin boucle X de lecture de tab
c---------------------------------------
c                           boucle X pour calcul des correlations

      do60 ncli=1,nfm-1
      type *,'Dans boucle60 ncli : ',ncli

        do jp=1,jrsom
        jp1=jp-jrsom+jm
        jp2=jp
          do i=1,im
          tab1(i,jp)=tab(i,jp1,ncli)
          tab2(i,jp)=tab(i,jp2,ncli+1)
          enddo
        enddo
c
      type *,'avant cortab' !!!!!!!!!!!!!!!!!!!!!!!!!!!
      call CORTAB(im,jrsom,ijc,corstep,jdel,iminus,jminus,
c          ------
     1 tab1,tab2,cor,icor(ncli+1),jcor(ncli+1),corm(ncli+1),milcor,
     2 kedge(ncli+1))
         type *,'apres cortab' !!!!!!!!!!!!!
           type *,'ncli,corm(ncli+1)',ncli,corm(ncli+1)
        if(corm(ncli+1).eq.0.)then
        nm=ncli
        endif
c
40    format(12f6.3)
      ijcm=2*ijc+1
        if(kedge(ncli+1).ne.0)then
          nm=ncli
        write(3,52)ncli
        type 52
        endif

50    format(' ncli,  icpas,  icor, jcor, corm',4i6,f8.3)
51    format(12i6)
52    format(' ncli=' i3,
     1       ' kedge=1 : jminus too large  or  ijc too small')
53    format(' correlations')
      ica=1
      icb=2*ijc+1
      icpas=1+(icb-1)/10
      write(3,50)ncli,icpas,icor(ncli+1),jcor(ncli+1),corm(ncli+1)
      type 50,ncli,icpas,icor(ncli+1),jcor(ncli+1),corm(ncli+1)

      write(3,53)
      write(3,51)0,(i,i=ica,icb,icpas)
      do j=ica,icb,icpas
        do i=ica,icb,icpas
        kcor(i)=1000.*cor(i,j)
        enddo
      write(3,51)j,(kcor(i),i=ica,icb,icpas)
      enddo
60    continue

61    format(//' Nombre de cliches correles ='i3//)
      type *,nm !!!!!!!!!!!!!!61->*
      write(3,61)nm
c---------------------------------
c                                   calcul de i1(n) et jr(n)
62     if(nm.eq.1)then
       i1(1)=1
       i1i2=im-1
       jr(1)=0
        type *,'i1,i1i2,jr', i1(1),i1i2,jr(1)
       return
      endif
c
      kpiv=nm-1
        do n=2,nm
        xcor(n)=icor(n)-ijc
        ycor(n)=jcor(n)-ijc
        enddo

      if(maxdecal.ne.0) then
      call zerocorr(ijc,icor,jcor,kpiv,xcor,ycor)
c          --------
      endif

        sr=0.
        jr(1)=0
        sjr=0.
      do n=2,nm
        sr=sr+jrsom+ycor(n)
        jr(n)=sr-sjr+0.5
        sjr=sjr+jr(n)
c           jrsom=jm-pixstep
      enddo
c
      i1(1)=iminus+1+maxdecal
      sx=i1(1)
      si=i1(1)
      do n=2,nm
        sx=sx+xcor(n)
        i1(n)=i1(n-1)+sx-si+0.5
        si=si+i1(n)-i1(n-1)
      enddo
c
      imin=1000
      imax=0
      do n=1,nm
      imin=min0(imin,i1(n))
      imax=max0(imax,i1(n))
      enddo             
c
c      ipiv=maxdecal+iminus+1
c      if(ipiv.gt.imin)imin=ipiv
c      ipiv=im-maxdecal-iminus
c      if(imax.lt.ipiv)imax=ipiv
c
c      i1i2=imax-imin
       i1i2=im-2*iminus-1-imax+imin
c      i1i2=i1i2-imax+imin
       if(imin.lt.iminus+1)then
        idell=iminus+1-imin
        do n=1,nm
        i1(n)=i1(n)+idell
        enddo
       endif
c
c
70    format('ncli,  i1,i1i2,jr,kedge, corm(ncli) =' i4,2x4i4,f8.3)
       do ncli=1,nm
       write(3,70)ncli,i1(ncli),i1i2,jr(ncli),kedge(ncli),corm(ncli)
       enddo
      return
      end
c-------------
      subroutine CORTAB(im,jrsom,ijc,corstep,jdel,iminus,jminus,
c                ------
     1                  tab1,tab2,cor,icc,jcc,corm,milcor,kedge)
      integer*2 tab1(1024,128),tab2(1024,128)
      dimension cor(101,101)
      integer corstep
c
1     format(10i5)
2     format(/)
c
c   differences premieres sur i
      idel=corstep
c
4     format(20i4)
5     format(' tab1: premiere et derniere lignes')
6     format(' tab2: premiere et derniere lignes')
  
c      type *,'im,jrsom,ijcormax',im,jrsom,ijc
c

      den=im-idel
c
      do j=1,jrsom-jdel
        si1=0.
        si2=0.
c
        do i=1,im-idel
        kt1=tab1(i,j)
        kt2=tab1(i+idel,j+jdel)
          kpiv=kt2-kt1
          if(kt1.eq.32000.or.kt2.eq.32000)kpiv=0
          tab1(i,j)=kpiv
            si1=si1+kpiv
        kt1=tab2(i,j)
        kt2=tab2(i+idel,j+jdel)
          kpiv=kt2-kt1
          if(kt1.eq.32000.or.kt2.eq.32000)kpiv=0
          tab2(i,j)=kpiv
            si2=si2+kpiv
        enddo

          si1=si1/den
          si2=si2/den
        do i=1,im-idel
        tab1(i,j)=tab1(i,j)-si1
        tab2(i,j)=tab2(i,j)-si2
        enddo
c
      enddo
c
        kj1=1+jminus
        kj2=jrsom-jminus-jdel
        kj12=kj2-kj1
      write(3,5)
      do j=kj1,kj2,kj12
      write(3,4)(tab1(i,j),i=1,im-idel,50)
      enddo     

      write(3,6)
      do j=kj1,kj2,kj12
      write(3,4)(tab2(i,j),i=1,im-idel,50)
      enddo

c
      do i=1,101
      do j=1,101
      cor(i,j)=0.
      enddo
      enddo
c
c      jcopas=1+(jrsom-2*jminus)/20
       jcopas=1
      ijcm=2*ijc+1
        do jc=1,ijcm
        do ic=1,ijcm
        som=0.
        s1=0.
        s2=0.
        nt=0
          do12 j=kj1,kj2,jcopas
          jp=j+jc-ijc
          if(jp.lt.kj1.or.jp.gt.kj2)goto12
          somj=0.
          sj1=0.
          sj2=0.
            do10 i=iminus+1,im-idel-iminus,corstep
            ip=i+ic-ijc
            if(ip.lt.iminus+1.or.ip.gt.im-idel-iminus)goto10
            t1=tab1(i,j)
            if(t1.ge.32000.)goto10
            t2=tab2(ip,jp)
            if(t2.ge.32000.)goto10
            nt=nt+1
c              if(nt.gt.2000)goto10
            somj=somj+t1*t2
            sj1=sj1+t1**2
            sj2=sj2+t2**2
10          continue
          som=som+somj
          s1=s1+sj1
          s2=s2+sj2
12        continue
        piv=s1*s2
        if(piv.ne.0.)cor(ic,jc)=som/sqrt(piv)
      enddo
      enddo
c
      corm=0.
      do ic=1,ijcm
      do jc=1,ijcm
      corm=amax1(corm,cor(ic,jc))
        if(corm.eq.cor(ic,jc))then
        icc=ic
        jcc=jc
        endif
      enddo
      enddo
15    format(12f6.3)
c
      kedge=0
      if((icc-1)*(icc-ijcm)*(jcc-1)*(jcc-ijcm).eq.0)kedge=1
      if((corm.lt.float(milcor)/1000.).or.(kedge.eq.1))then
        icc=ijc
        jcc=ijc
        corm=0.
      endif
      return
      end
c-----------------------
      subroutine LECJ(jname,sundec,is,js,nm,jminus,jm,vmoy,
c                ----
     1                i1,jr,corm)
c                                           vmoy pour centrages q/p
      dimension corm(150)
      integer*4 i1(150),jr(150)
      character*20 jname
      character com32*32,com20*20
      integer sundec
1     format(a)
2     format(1x,a)
3     format(5i6)
4     format(2i6,f6.3)
5     format(f12.1)


      call openold20sf(jname,sundec,60)
c          -----------
      read(60,1)com20
      write(3,1)com20

      read(60,1)com32
      write(3,1)com32
      read(60,5)vmoy
      write(3,5)vmoy
      read(60,1)com32
      write(3,1)com32
      read(60,3)is,js,nm,jminus,jm
      write(3,3)is,js,nm,jminus,jm
      read(60,1)com32
      write(3,1)com32
          type *,com32
        do n=1,nm
        read(60,4)i1(n),jr(n),corm(n)
          type 4,i1(n),jr(n),corm(n)
          write(3,4)i1(n),jr(n),corm(n)
        enddo

10      close(unit=60)
c
      return
      end
c----------------------------
c      subroutine CALMOY(fname,nf1,n1,n2,ic1,ic2,ic3,ic4,cmoy)
c                                                                NON UTILISE
c      integer*4 cmoy(100,2,30),norm(100,2,30),head(512)
c      integer*2 lec(2048)
c
c      nfa=nf1+n1-1
c      nfb=nf1+n2-1
c      do l=1,30
c      do k=1,2
c      do j=1,100
c        cmoy(j,k,l)=0
c        norm(j,k,l)=0
c      enddo
c      enddo
c      enddo
c
c1     format(a)
c      do 10 nf=nfa,nfb

c      rewind(iud)
c      read(iud)head
c      im=head(2)
c      jm=head(3)
c      lm=head(4)
c
c      if(ic1.eq.0)then
c        ic1=1
c        ic2=im/2
c        ic3=ic2+1
c        ic4=im
c      endif
c
c      do l=1,lm
c      do j=1,jm
c        read(10)(lec(i),i=1,im)
c        do i=ic1,ic2
c        if(lec(i).ne.32000)then
c           cmoy(j,1,l)=lec(i)+cmoy(j,1,l)
c           norm(j,1,l)=1+norm(j,1,l)
c        endif
c        enddo
c        do i=ic3,ic4
c        if(lec(i).ne.32000)then
c           cmoy(j,2,l)=lec(i)+cmoy(j,2,l)
c           norm(j,2,l)=1+norm(j,2,l)
c        endif
c        enddo
c      enddo
c      enddo
c
c10    continue
c
c      do l=1,lm
c      do k=1,2
c      do j=1,jm
c        kpiv=norm(j,k,l)
c        if(kpiv.ne.0)then
c          cmoy(j,k,l)=float(cmoy(j,k,l))/float(kpiv) 
c        else
c          cmoy(j,k,l)=32000
c        endif
c      enddo
c      enddo
c      enddo
c
c      return
c      end
c----------------------------------------------------------------
      subroutine scorond(nfm,npo,dnanfp,sundec,jm,im,
c                -------
     1      lec,ltz,ltrjp,ond,poids,
     1      iminus,jminus,milsig,ides,tabmoy,sigma,ll,   ilissond) 
c                                                  ici sans dim
      integer*4 head(512)
      integer*2 lec(2048)
      dimension ond(1024,128,3),poids(1024,128,3),tabmoy(30),
     1          sigma(30)
      character*20 dnanfp(2000,2)
c
      i1=1+iminus
      i2=im-iminus
      j1=1+jminus
      j2=jm-jminus
      ltrj=iabs(ltrjp)
          type*,' scorond: i1,i2,j1,j2,ltrj ',i1,i2,j1,j2,ltrj

        ablec1=-100000.
        ablec2=100000.
c 
      do j=1,128
      do i=1,1024
      ond(i,j,ides)=0.
      poids(i,j,ides)=0.
      enddo
      enddo
c                              1   calcul ablec
c                              2   calcul tabmoy
c                              3   calcul ond
c----------------------------------------------------   boucle moyond 
       do90 moyond=1,3 !!!!!!!!!!!!
      type *,'                        valeur de moyond : ',moyond 
5     format(' moyond,milsig=',2i6)
      write(3,5)moyond,milsig
c--------------------------------------------------
c                                        etapes 2,3
            if(moyond.ne.1)then
            tabmoy(ltz)=sijly/sijl
            sigma(ltz)=sijly2/sijl-tabmoy(ltz)**2
            sigma(ltz)=sqrt(sigma(ltz))
              ablec1=-sigma(ltz)*float(milsig)/1000.
              ablec2=-ablec1
          type *,'ltz,tabmoy,sigma,milsig,ablec1,ablec2'
          type *, ltz,tabmoy(ltz),sigma(ltz),milsig,ablec1,ablec2
10    format(' ltz,tabmoy,sigma,milsig,ablec1,ablec2'/
     1        i3,2f8.0,i5,2f8.0)
      write(3,10)ltz,tabmoy(ltz),sigma(ltz),milsig,ablec1,ablec2
            endif
c----------------------------------------------------
c                                        etapes 1,2,3
c----------                       boucle X
           nm=nfm
         type *,'avant do80 dans scorond nm,npo',nm,npo

      do 80 ncli=1,nm

      iud=51
      call openold20(dnanfp(ncli,npo),sundec,iud)
c          ---------
      rewind(iud)
      read(iud)head

      lm=head(4)
c      type *,'lm :',lm

c        type *,'ncli',ncli
12    format(' ncli=',i4)
        write(3,12)ncli
c-------------------------------------------------
      if(moyond.eq.3)goto40
c                                       etapes 1,2
         sijl=0.
         sijly=0.
         sijly2=0.

        do35 l=1,ltz
                  sijy=0.
                  sijy2=0.
                  sij=0.

        do30 j=1,jm
           read(iud)(lec(i),i=1,im)
c               ---
c      write(3,*)' dme/scorond: l,j,lec(10) ',l,j,lec(10)
           if(l.ne.ltz)goto30
           if(j.lt.j1.or.j.gt.j2)goto30
                  si=0.
                  siy=0.
                  siy2=0.
                     if(ll.gt.0)then
                     suplec=0
                     do i=i1,51
                       piv=lec(i)
                       suplec=amax1(piv,suplec)
                     enddo
                     do i=i2-51,i2
                       piv=lec(i)
                       suplec=amax1(suplec,piv)
                     enddo
                     endif
c                             if(j.eq.jm/2)type *,'suplec',suplec
15    format(' suplec,ll',f8.0,i5)
c                             if(j.eq.jm/2)write(3,15)suplec,ll
          do i=i1,i2
             if(lec(i).ne.32000)then
             piv=lec(i)
             test=-1.
c   -------------------------
               if(moyond.ne.1)then
               piv2=piv-tabmoy(ltz)
               test=(piv2-ablec1)*(piv2-ablec2)
               endif
c   -------------------------                          protus
c        if(ll.lt.0.or.(ll.gt.0.and.piv.gt.suplec/4..and.test.le.0.))
c     1        then
                   siy=siy+piv
                   siy2=siy2+piv*piv
                   si=si+1.
c              endif
c----------------------------
             endif
          enddo

                   sij=sij+si
                   sijy=sijy+siy
                   sijy2=sijy2+siy2 
30     continue
              if(l.ne.ltz)goto35

                   sijl=sijl+sij
                   sijly=sijly+sijy
                   sijly2=sijly2+sijy2
35    continue
c
38    format(' sijl,sijly,sijly2=',3e12.3)
      write(3,38)sijl,sijly,sijly2

      goto70     !!!!!             ??????
c------------------------------------------------
c                                         etape 3   calcul ond
40     type *,'passage a 40 avec ncli :',ncli
       continue
c
        do l=1,ltz  
              jpm=1
          do j=1,jm
            read(iud)(lec(i),i=1,im)
c                ---
            if(l.ne.ltz)goto50
c
               jtest=0
               if(j.ge.j1.and.j.le.j2)jtest=1

c              type*,' scorond: ltrj ',ltrj
            kpiv=(100*(j-j1))/ltrj
            jp=j-kpiv*float(ltrj)/100.
            jpm=max0(jp,jpm)
c
             do i=i1,i2
             if(jtest.eq.1)then
               piv=lec(i)-tabmoy(ltz)
45    format(' lec,tabmoy,ablec1,ablec2='i6,3f8.0)
      if(i.eq.im/2.and.j.eq.jm/2)write(3,45)lec(i),tabmoy(ltz),
     1                           ablec1,ablec2
               if((piv-ablec1)*(piv-ablec2).le.0.)then
               poids(i,jp,ides)=poids(i,jp,ides)+1.
               ond(i,jp,ides)=ond(i,jp,ides)+piv
               endif
             endif
             enddo
50         enddo
         enddo           
      type *,'avant 70'

70      close(unit=iud)

80    continue
       type *,'apres 80 continue'
c--------------------------------------------                                   
      if(moyond.ne.3)goto90
c                                                 etape 3
        type*,' dme: moyond,jpm,ond(1,1,1),poids(1,1,1)',
     1               moyond,jpm,ond(1,1,1),poids(1,1,1)
      do i=i1,i2
        do jp=j1,jpm
         if(poids(i,jp,ides).eq.0.)then
         ond(i,jp,ides)=0.
         else
         ond(i,jp,ides)=ond(i,jp,ides)/poids(i,jp,ides)
c        --------------
         endif
        enddo
      enddo
c
      if(j2.gt.jpm)then
      do i=i1,i2
        do j=jpm+1,j2
        kpiv=(100*(j-j1))/ltrj
        jpp=j-kpiv*float(ltrj)/100.
        ond(i,j,ides)=ond(i,jpp,ides)
        poids(i,j,ides)=poids(i,jpp,ides)
        enddo
      enddo
      endif
c
90    continue
c                                           fin boucle moyond
c-------------------------------------------------
c                                   Moy glissantes 
      call sgliss(ond,i1,i2,j1,j2,ides,ilissond)
c          ------
      type *, ' Fin sgliss'
100   format(/' ltz,j1,j2,ltrj,jpm,ides='6i4,
     1      / ' j  ond(1,j) ond(im/2,j) ond(im,j)   poids...')
101   format(i4, 3f8.3, 6x,3f6.0,2f8.0)
      write(3,100)ltz,j1,j2,ltrj,jpm,ides
        ims2=im/2
      do j=1,jm,50
      write(3,101)j,ond(1,j,ides),ond(ims2,j,ides),ond(im,j,ides),
     1        poids(1,j,ides),poids(ims2,j,ides),poids(im,j,ides)
      enddo
c
      return
      end
c************************************************************
      subroutine scorond3(nfm,npo,dnanfp,sundec,jm,im,lm,
c                -------
     1      lec,ltrjp,sig,sigmoy,
     1      iminus,jminus,milsig,ll,ilissond) 
      integer*4 head(512)
      integer*2 lec(2048)
      dimension sig(128,30,6),sigmoy(30,6),pds(128,30),ll(30),
     1          xdes(128),ydes(128)
      character*20 dnanfp(2000,2)
      character*1 index(9)
      double precision s0,sx,sx2,sx3,sx4,sx5,sx6
      data index/'1','2','3','4','5','6','7','8','9'/

      write(3,*)' scorond3: jm,im,lm ',jm,im,lm

        mm=6
      i1=1+iminus
      i2=im-iminus
      j1=1+jminus
      j2=jm-jminus
      ltrj=iabs(ltrjp)
          
      do l=1,30
        do m=1,mm
          sigmoy(l,m)=0.
        enddo
      do j=1,128
        pds(j,l)=0.
        do m=1,mm
          sig(j,l,m)=0.
        enddo
      enddo
      enddo
 

c-------------------------------   calculs des moments centres
      do86 moyond=1,2

c-------------------------------------                boucle nfm

      do 80 ncli=1,nfm

      iud=51
      call openold20(dnanfp(ncli,npo),sundec,iud)
c          ---------
      rewind(iud)
      read(iud)head

      lm=head(4)
c      type *,'lm :',lm

c        type *,'ncli',ncli
12    format(' ncli=',i4)
c        write(3,12)ncli
c-------------------------------------------------


        do35 l=1,lm
c               write(3,*)' scorond3: moyond,l,xmoy(1) ',
c     1                              moyond,l,sig(1,l,1)

        do30 j=1,jm
           read(iud)(lec(i),i=1,im)
c               ---
           if(ll(l).ne.1)goto30
           if(j.lt.j1.or.j.gt.j2)goto30

                if(moyond.eq.2)xmoy=sig(j,l,1)
           s0=0
           sx=0
           sx2=0
           sx3=0
           sx4=0
           sx5=0.
           sx6=0.

         do25 i=i1,i2
             if(lec(i).ne.32000)then
             x=lec(i)            

            if(moyond.eq.1)then
              s0=s0+1.
              sx=sx+x

            else
              x=x-xmoy
                      x2=x*x
              sx2=sx2+x2
                      x3=x2*x
              sx3=sx3+x3
                      x4=x3*x
              sx4=sx4+x4
                      x5=x4*x
              sx5=sx5+x5
                      x6=x5*x
              sx6=sx6+x6

            endif
            endif

25       continue

        if(moyond.eq.1)then
        pds(j,l)=pds(j,l)+s0
        sig(j,l,1)=sig(j,l,1)+sx

        else
          sig(j,l,2)=sig(j,l,2)+sx2
          sig(j,l,3)=sig(j,l,3)+sx3
          sig(j,l,4)=sig(j,l,4)+sx4
          sig(j,l,5)=sig(j,l,5)+sx5
          sig(j,l,6)=sig(j,l,6)+sx6
        endif


30     continue

35    continue

70      close(unit=iud)

80    continue
       type *,'apres 80 continue'
c-----------
      do85 l=1,lm
        if(ll(l).ne.1)goto85

      do j=1,jm
           if(pds(j,l).eq.0)then
                     pds(j,l)=1.
                     type*,' attention pds=0'
                     write(3,*),' attention pds=0'
            endif

        if(moyond.eq.1)then
            if(pds(j,l).eq.0)then
                     pds(j,l)=1.
                     type*,' attention pds=0'
                     write(3,*),' attention pds=0'
            endif
        sig(j,l,1)=sig(j,l,1)/pds(j,l)
        
        else
          do84 m=2,mm
          piv=sig(j,l,m)/pds(j,l)
          expo=1./float(m)
            if(piv.lt.0.)then
            write(3,*)' attention moment.lt.0 '
            type*,' attention moment.lt.0 '
            sig(j,l,m)=1.
            goto84
            endif
          sig(j,l,m)=piv**expo
84        continue
        endif
      enddo

85    continue
86    continue
c-----------------   fin boucle moyond

87    format(/'     l     m  sigmoy   sig(j=1,jm,10)')
88    format(2i6,           f8.0, 2x,    10f6.0)

      do l=1,lm
      do m=1,mm
        piv=0.
          do j=5,jm-4
          piv=piv+sig(j,l,m)
          enddo
        sigmoy(l,m)=piv/float(jm-8)
      enddo
      enddo
c-------
        do l=1,lm
                  write(3,87)
        do m=1,mm
        write(3,88)l,m,sigmoy(l,m),(sig(j,l,m),j=1,jm,10)
        enddo
        enddo
c--------
      call pgadvance
      call pgvport(0.1,0.9,0.1,0.8)
      call pglabel('coupes/j pour l=1,lm','moments',
     1             'moyennes (1) et moments centres (2,3,4) ')

      jm2=jm/2
      piv1=sig(jm2,lm,1)
      piv2=sig(jm2,lm,mm)
      ymax=1.2*amax1(piv1,piv2)

c      piv1=sig(jm2,1,1)
c      piv2=sig(jm2,lm,3)
c      ymin=amin1(piv1,piv2)

      do j=1,jm
      xdes(j)=j
      enddo

      xjm=float(jm)

      do l=1,lm
        xmin=0.1+(l-1)*0.8/float(lm)
        xmax=0.1+l*0.8/float(lm)
        call pgvport(xmin,xmax,0.1,0.8)

        call pgwindow(1.,xjm,0.,ymax)
c      if(l.eq.1)then
        call pgbox('bcnts',20.,2,'bcnts',100.,10)
c      else
c        call pgbox('bcnts',20.,2,'bcnts',100.,10)
c      endif

        do m=1,mm
                          xind=jm2
                          yind=sig(jm2,l,m)+4.
                          call pgtext(xind,yind,index(m))
        call pgsls(m)
          do j=1,jm
          ydes(j)=sig(j,l,m)
          enddo
        call pgline(jm,xdes,ydes)
        enddo
          call pgsls(1)

      enddo
c-------------
          
      return
      end
c***************************************
      function cormom(y,sig,sigmoy,j,l)
c              ------

      dimension sig(128,30,6),sigmoy(30,6)
      mcor=2

      cormom=sigmoy(l,1)+(y-sig(j,l,1))
     1                   *(sigmoy(l,mcor)/sig(j,l,mcor))

c      cormom=y*sigmoy(l,1)/sig(j,l,1)

      return
      end

c*************************************************************
      subroutine sgliss(x,i1,i2,j1,j2,ides,ilissond)
c                ------
      dimension x(1024,128,3),xx(1024)
      liss=ilissond
      liss2=2*liss+1
      fliss2=float(liss2)
c
      type *,'i1,i2,j1,j2,ides',i1,i2,j1,j2,ides
        do j=j1,j2
          do i=i1,i2
          xx(i)=x(i,j,ides)
          if(xx(i).eq.0..and.i.ne.i1)xx(i)=xx(i-1)
          enddo
        s=0.
          do i=i1,i1+liss2-1
          s=s+xx(i)
          enddo
          x(liss+i1,j,ides)=s/fliss2
             if(j.eq.1)type *,' i,somme',i,s
c
        do i=i1+liss+1,i2-liss
        s=s+xx(i+liss)-xx(i-liss-1)
        x(i,j,ides)=s/fliss2
        enddo
c
        do i=i2-liss+1,i2
        x(i,j,ides)=x(i2-liss,j,ides)
        enddo

      enddo
      return
      end
c*****************************************************************
      subroutine sdescor(ijc,icor,jcor,corm,nm,xcor,ycor)
c                -------
      dimension x(150),y(150),icor(150),jcor(150),corm(150),
     1          xcor(150),ycor(150),x1(150)

      call pgvport(0.1,0.2,0.1,0.8)
        ynm=nm
        if(nm.lt.2)ynm=2.
      call pgwindow(0.5,1.2,1.,ynm)
      call pgbox('bcnst',0.5,5,'bcnst',0.,0)
c
      if(nm.lt.2)return
      nm1=nm-1
      do n=1,nm1
      y(n)=n+1
      x(n)=corm(n+1)
      enddo
      call pgline(nm1,x,y)
c
      xmax=ijc-5
      call pgvport(0.2,0.4,0.1,0.8)
      call pgwindow(-xmax,xmax,1.,ynm)
      call pgbox('abcnst',10.,5,'abcst',0.,0)
c
      do n=1,nm
      x(n)=icor(n)-ijc
      x1(n)=xcor(n)
      y(n)=n
      enddo
      call pgsls(2)
      call pgline(nm,x1,y)
      call pgsls(1)
      call pgline(nm,x,y)
c
      call pgvport(0.4,0.6,0.1,0.8)
      call pgwindow(-xmax,xmax,1.,ynm)
      call pgbox('abcnst',10.,5,'abcmst',0.,0)
c
      do n=1,nm
      x(n)=jcor(n)-ijc
      x1(n)=ycor(n)
      enddo
      call pgsls(2)
      call pgline(nm,x1,y)      
      call pgsls(1)
      call pgline(nm,x,y)
c
      return
      end
c********************************************************
      subroutine sdesond(ond,poids,im,jm,den,tabmoy)
c                -------   (tabmoy pour normalisation intensites dessins)
c dessin des 2 derniers  tableaux l
c    3 coupes en i de chacun
c          
      dimension icoupe(3),x(128),y(128),ond(1024,128,3),
     1          poids(1024,128,3),xp(128),tabmoy(30)
      call pgadvance
c
          xjm=jm
          do j=1,jm
          y(j)=j
          enddo
c
      type *,'avant boucle50' !!!!!!!!!!!!!
      type *,'den : ',den !!!!!!!!!!

      call pgvport(0.1,0.9,0.1,0.8)
      call pglabel('coupes .1 .5 .9*im',
     1    'Tableau 1 (int)         Tableau 2 (vit)','Corrections')
      do50 idess=1,2
      ypiv1=0.1+0.35*(idess-1)
      ypiv2=ypiv1+0.35
c
      do icoup=1,3
        xpiv1=0.1+(icoup-1)*0.8/3.
        xpiv2=xpiv1+0.8/3.
      call pgvport(xpiv1,xpiv2,ypiv1,ypiv2)
        if(idess.eq.1)then
        piv1=-0.30
        piv2=0.30
        else
        piv1=-150.
        piv2=150.
        endif
      call pgwindow(piv1,piv2,1.,xjm)
      if(idess.eq.1)then
      call pgbox('bcnts',0.,0,'bcnts',0.,0)

      else
      call pgbox('bcmts',0.,0,'bcnts',0.,0)

      endif
c
        i=im*(0.1+0.4*(icoup-1))
          do j=1,jm
          x(j)=ond(i,j,idess)
c             type *,'tabmoy(idess) :',tabmoy(idess) !!!!!!!!!!!!!
            if(idess.eq.1)x(j)=x(j)/tabmoy(idess)
          xp(j)=poids(i,j,idess)/den
          enddo
        call pgline(jm,x,y)
c
      call pgwindow(0.,2.,1.,xjm)
      call pgsls(2)
      call pgline(jm,xp,y)
      call pgsls(1)
c
      enddo
50    continue
      return
      end
c------------------------------------------
      subroutine zerocorr(ijc,icor,jcor,nm1,xcor,ycor)
c                --------
      dimension icor(150),jcor(150),xcor(150),ycor(150)
c
        nm=nm1+1
      if(nm1.le.2)then
        do n=2,nm
        xcor(n)=icor(n)-ijc
        ycor(n)=jcor(n)-ijc
        enddo
      xcor(1)=0.
      ycor(1)=0.
      icor(1)=ijc
      jcor(1)=ijc
      return
      endif
c
      na=2
      nb=nm
      ab=nb-na+1.
      siab=0.
      sjab=0.
c
      icor(1)=ijc
      jcor(1)=ijc
      xcor(1)=0.
      ycor(1)=0.
        do n=1,nm1
        xcor(n+1)=icor(n+1)-ijc
        ycor(n+1)=jcor(n+1)-ijc
        enddo
c
      do n=na,nb
      siab=siab+xcor(n)
      sjab=sjab+ycor(n)
      enddo
        siab=siab/ab
        sjab=sjab/ab
c
      do n=2,nm
      xcor(n)=xcor(n)-siab
      ycor(n)=ycor(n)-sjab
      enddo
c
10    format('   icor-ijc/ xcor/ jcor-ijc/ ycor')
      write(3,10)
11    format(i10,f6.2,i10,f6.2)
      do n=1,nm
        kpiv1=icor(n)-ijc
        kpiv2=jcor(n)-ijc
      write(3,11)kpiv1,xcor(n),kpiv2,ycor(n)
      enddo
c      
      return
      end
c-------------------------------------------
      subroutine srejec(tab,is,js,milsig,ymoy)
c                ------
      integer*2 tab(400,4000)
c
      double precision sy,sy2,pds
c
        ymoy=0.
        ec=0.
      do n=1,5
      sy=0.
      sy2=0.
      pds=0.
      do i=1,is,5
      do j=1,js,5
        y=tab(i,j)
        if(y.ne.32000.)then
          if(n.eq.1.or.abs(y-ymoy).lt.ec) then
          sy=sy+y
          sy2=sy2+y*y
          pds=pds+1.
          endif
        endif
      enddo
      enddo
c
        if(pds.eq.0.)then
          type *,' Rejection impossible pour le calcul de moyennes'
          return
        endif
c
        ymoy=sy/pds
        var=sy2/pds-ymoy*ymoy
        var=sqrt(var)
        ec=(var/1000.)*float(milsig)
        type *,' rejec: n,ymoy,pds,var,ec ='
        type *,         n,ymoy,pds,var,ec
      enddo
c
      return
      end
c************************************************
 

