c******************************************************************************
c  Programme cmf - UNIX
c******************************************************************************
        subroutine cmf(cmfdat,npo,nw1,nw2,win,kdecal,gna,cna,fna,fl,
     1                  tab,tabb) 
        dimension x(2000),y(2000),
     1            cr(10),profm(2000,8),prof100(101,11),coef(11),
     1            cf(11,2),
     1            kz(512),profmm(2000)
        dimension xr(11,3,2),yr(11,3,2),tabp(128),tabq(128),
     1            yline(1024,11,8),grad(11),ylinem(1024,11),
     2            trj(8),ncw(8),kmin(8),kdecal(8)
        integer*4 displ,nrej(90),smoothi,smoothj,day,time,curv,
     1            win(8),winp, som1, som2
        integer*2 tab(400,60,32),tabb(400,2,8),
     1            rec(1024),fl(400,60,8),pm(2000)
        character COMM*80,nom*15,fname*20,fnameps*23,
     1            cmfdat*6,flat*6,end*3,titre*64,canal*2,
     2            ftype*9,hdeb*6,jdate*9
        character*20 gna(8),cna(8),fna(8)
c        write(czfdat(1:6),'(a6)')'czfdat'  !!!!!!!!!!!!!!puis character czfdat
        write(flat(1:6),'(a6)')'flat  '
        write(end(1:3),'(a3)')'end'

      type*,' cmf1: kdecal ',kdecal

c Input: .DAT and C-file
c Output: F-file = Header = 15 values (6,im,jm,nm, im,jm,1,
c                                      7,day, time, lbda,
c                                      jt100, ja100, jb100, jz100)
c                  Table Flat field = im*jm*nm
c                  Mean Profile = im values
C
C     8=max number channels
C-----------------------------
c	CALL HORODAT(JHD1,JHD2,JHD3,JHD4)
C-----------------------------
C LECTURE SUR .DAT
1	FORMAT(A)
2	FORMAT(2x,A)
3       format(1x)
4       format(1x,a6)
12	FORMAT(12I6)
13      format(11f7.1)
14	FORMAT(2X,12I6)
16	FORMAT(' donnees sur .par')
22	FORMAT(/)
24	FORMAT(1H1)
c
      rewind(8)
      call SNAME(8,cmfdat, ier) 
      write(3,16)
      READ(8,1)COMM
      write(3,2)COMM
      READ(8,12)nm,mupris,mustep,smoothi,smoothj,il1p,il2p,isym
      write(3,14)nm,mupris,mustep,smoothi,smoothj,il1p,il2p,isym

         n1=1
         n2=nm

          ncm=(1+nm)/2
c         ---             canal de reference pour  cmd.f

c  smoothi=1     constant zero function versus i
c  smoothi=2     linear zero function versus i
c  smoothj=1     constant zero function versus j
c  smoothj=2     linear zero function versus j
c  if mupris.ne.0
c                 read file g....  before c....
c                 mupris= translation by PRISms (microns)
c                 mustep= STEP between slits before prisms
c  isym   =1      le profil moyen est remplace par le profil symetrise

c                                                     boucle nw
c                                                     ---------
      do 260 nw=nw1,nw2

      rewind(8)
      call sname(8,cmfdat,ier)
      call sname(8,'cmf1  ',ier)
        read(8,1)comm
        write(3,2)comm
                       
          do n=1,8
        read(8,12)winp, nc,inclin,curv,iliss,jparliss,lissprofm,
     1             jt100,ja100,jb100,jz100
        write(3,14)winp,nc,inclin,curv,iliss,jparliss,lissprofm,
     1             jt100,ja100,jb100,jz100
          if(winp.eq.win(nw))goto30
          enddo

30        ncw(nw)=nc
c         -------    canal ou sera recherche le centre raie (option inclin=1)

c        nc=canal ou sera cherche le centre raie (option incli=1)
c        inclin=0   forme des raies lue par ja100,jb100,jz100
c                      jz100 = 100 fois le j  de la raie en (1+im)/2
c                      ja100 = 100 fois le J  de la raie en i=1   
c                      jb100 = 100 fois le j  de la raie en i=im 
c        jt100 = 100 fois
c                    la translation en j correspondant a la meme difference
c                    en lambda qu'entre deux canaux successifs (initiaux)
c        curv = 0/1  (prise en compte de la courbure de raie)
c        iliss    lissage en i sur 2*iliss+1 avant detection centre raie
c        jparliss     "      j parabolique        "   "        "     "
c        lissprofm    "      k parabolique  sur le profil moyen qui sert
c                                          pour calculer fl

      do n=1,8
      do j=1,128
      do i=1,400
       tab(i,j,n)=0
       fl(i,j,n)=0
      enddo
      enddo
      enddo
c
      do i=1,512
      kz(i)=0
      enddo
c      type*,' cmf3: kdecal ',kdecal
c*****************
c                                                     images TAB
c                    (les backgrounds sont deja soustraits dans bmc)
c      type*,' cmf3:nw,cna ',nw,cna(nw)

      iuc=26
      call openold20(cna(nw),sundec,26)
c          ---------
c      type*,' cmf3a ',kdecal


      rewind(iuc) !!!!!!!!!!!!
      read(iuc)na3,im,jm,nm,imm,jmm,nmm,naxis3,naxis4,nd,day,time,
     1 iline,lbda,(kz(i),i=15,512) 
      type *,im,jm,nm,day,time,iline,lbda
 
      jcentre=(1+jm)/2
c     -------

110   format(1xi6,8x)
120   format(9xi6)
125   format(' lecc: header z ou c'/6i8/6i8/6i8) !!!!!!!!!!!!!!!10 -> 11
      write(3,125)na3,im,jm,nm,imm,jmm,nmm,naxis3,naxis4,nd,day,time,lbda
c
      do n=1,nm
        do j=1,jm
        read(iuc)(rec(i),i=1,im) !!!!!!!!!!!!!!!
          do i=1,im
          tab(i,j,n)=rec(i)
c         ----------
          enddo
        enddo

        write(9,*)' cmf lit z: n, tab(187,25,n) ',n,tab(187,25,n)
      enddo

c          type *,' lecc2 ',kdecal

      im2=im/2
      jm2=jm/2
      write(3,*)' tab ',(tab(im2,jm2,n),n=1,nm)
      type    *,' tab ',(tab(im2,jm2,n),n=1,nm)  

c      type *,' nm,im ',nm,im
c      do jgd=1,2
c        do n=1,nm
c        read(iuc)(rec(i),i=1,im)  !!!!!!!!!!!!!!!!!
c          do i=1,im
c          tabb(i,jgd,n)=rec(i)
c      if(jgd.eq.1.and.n.eq.1)write(3,*)jgd,n,i,kdecal(1),tabb(i,jgd,n)
c          enddo
c        enddo
c      enddo
c        type*,' lecc3', kdecal

c      write(3,*)' backgrounds 1 ',(tabb(im2,1,n),n=1,nm)
c      write(3,*)' backgrounds 2 ',(tabb(im2,2,n),n=1,nm)

      close(unit=26)

c*****************
c                                     Calcul de TRJ
c                                     *************
         write(3,*)' jt100',jt100
      if(jt100.ne.0)then
      trj(nw)=jt100*0.01
      else

c lecture de vecjb et chajb
c            transl et largeur canaux  en unites fichier b
c
      iug=25
      open(unit=25,status='old',file=gna(nw)) !!!!!!!!!!!!!!!!!!
c                               -----------
         rewind(iug)    !!!!!!!!!!!!!!!
        do n=1,19
        read(iug,3) !!!!!!!!!!!!!!!
        enddo
      read(iug,13)xr !!!!!!!!!!!!!!
      read(iug,13)yr !!!!!!!!!!!!!!
      close(unit=iug) !!!!!!!!!!!!!!
c
      chajb=0.5*(yr(nc,1,2)-yr(nc,1,1)+yr(nc,3,2)-yr(nc,3,1))
      vecjb=0.125*(yr(nc+1,1,1)+yr(nc+1,3,1)+yr(nc+1,1,2)+yr(nc+1,3,2)
     1            -yr(nc-1,1,1)-yr(nc-1,3,1)-yr(nc-1,1,2)-yr(nc-1,3,2))
      type *,chajb,vecjb
c
      trj(nw)=float(jm-1)*mustep*vecjb/(float(mupris)*chajb)
c     -------
c      write(3,*)' nw, trj ',nw,trj
      endif

c******************
      type *,'calcul de yline100'    !!!!!!!!!!!!!!
c Calculation of YLINE100(1024,11,8)   + plot
c ----------------------------------     ----
c           seules seront utiles la pente (et la courbure) de YLINE

      i1=1
      i2=im
              im16=im-16
              im8=im-8
 
        xl=im-1
        yl=jm-1

      il1=im*float(il1p)/100.
      il2=1+im*float(il2p)/100.
                                write(3,*)' il1,il2 ',il1,il2
      iliss2=2*iliss+1
      iil1=il1+iliss
      iil2=il2-iliss
      iilm=iil2-iil1+1

      jmarge=jm/10
c                                      <---------   modifiable


      if(inclin.eq.0)then
c     ---------------

        xxl=xl*0.5
       al=0.01*(ja100+jb100)/2.
c      al=0.01*(ja100+jb100)/(2.*xxl**2)
       bl=0.01*(jb100-ja100)/(2.*xxl)
        do i=1,im
        xi=i-xxl
        yline(i,nc,nw)=al+bl*xi
        rec(i)=al+bl*xi
c        yline(i,nc,nw)=xi*(bl+xi*al)
        enddo
      goto200
c              le dessin sera le meme que dans l'option inclin=1
      endif

      if(inclin.eq.1)then
c     ---------------
c      type*,'kdecal ',kdecal            
      do i=il1,il2
             iip=i-il1+1
             y(iip)=i

      do j=1+jmarge,jm-jmarge
                ipiv1=i-iliss
                ipiv2=i+iliss 
                piv=0.
            do ip=ipiv1,ipiv2
             piv=piv+tab(ip,j,nc)
c                    ---
            enddo
           tabp(j)=piv/float(iliss2)
      enddo

          call parafit(tabp,1+jmarge,jm-jmarge,jparliss,tabq)
c              -------

      if(i.eq.im/2)type *,'parafit: i, tabp,tabq',tabp(20),tabq(20)

      piv1=100000.
            do j=jmarge+1,jm-jmarge
c              piv2=tab(i,j,nc)
              piv2=tabq(j)
            if(piv2.lt.piv1)jmin=j
              piv1=amin1(piv1,piv2)
            enddo
            piv1= tab(i,jmin-1,nc)
            piv2= tab(i,jmin,nc)
            piv3= tab(i,jmin+1,nc)
c       type *,'avant calcul',piv1,piv2,piv3                    
       rec(i)=jmin + 0.5*(piv1-piv3) / (piv1 - 2*piv2 + piv3)
       yline(i,nc,nw)= rec(i)
c       type *,'apres calcul'
       enddo
           write(3,*)' rec ',(rec(i),i=1,im,100)  
      endif
c--------------                    lissage yline

        im8=im-8
        im16=im-16
        jm2=jm/2
      iic=0.5*(il1+il2)
      iic1=iic+1
      type *,' iil1,iil2 ',iil1,iil2
      if(curv.eq.0)then
c        ---------
      call ZERO3(iil1,iic,iic1,iil2,0,0,0,rec,1,im,0,0,
     1           ar,br,cr,item,nrjm,nrej,y1,y2,ifin)
      cr(1)=br
      cr(2)=ar
      cr(3)=0.
      else
c
      call ZERO3(iil1,0,0,iil2,3000,5,3,rec,1,im,0,0,
     1           ar,br,cr,item,nrjm,nrej,y1,y2,ifin)     
      endif

      do i=1,im
      yline(i,nc,nw)=cr(1)+i*(cr(2)+i*cr(3))
      enddo
c                            rec=raie    yline=raie lissee

c-----------------------------
c                                      dessin A1 options 0 et 1

200   continue
      if(inclin.eq.0.or.inclin.eq.1)then
c        --------------------------
      call pgadvance
c          ---------
      call pgvport(0.1,0.5,0.1,0.8)
      call pgwindow(1.,yl,1.,xl)
        write(canal(1:2),'(i2)')nc
      call pgbox('bcnts',20.,2,'bcnts',0.,0)
      call pglabel(' ',' ',canal)

      do i=iil1,iil2
      iip=i-iil1+1
      x(iip)=rec(i)
      y(iip)=i
      enddo
      call pgline(iilm,x,y)
c          ------

      do i=1,im
      x(i)=yline(i,nc,nw)
      y(i)=i
      enddo
      call pgline(im,x,y)
c          ------
      endif
c---------------------------------


      if(inclin.eq.2)then
c     ---------------
      call pgadvance
c                                 option non testee
c   Recherche des canaux a abs(pente) maximum au centre du canal
      ic=(1+im)/2
      jc=(1+jm)/2
      jpiv1=jc-jm/10
      jpiv2=jc+jm/10
c
      gradp=-100000.
      gradm=100000.    
      do n=1,nm
      grad(n)=tab(ic,jpiv2,n)-tab(ic,jpiv1,n)
      type *,'n: ',n
      type *,' tab(ic,jpiv2,n): ',tab(ic,jpiv2,n)
      type *,' tab(ic,jpiv1,n): ',tab(ic,jpiv1,n)
      gradm=amin1(grad(n),gradm)
      gradp=amax1(grad(n),gradp)
      enddo

      do n=1,nm
      if(grad(n).eq.gradm)nfte2=n
      if(grad(n).eq.gradp)nfte1=n
      enddo

      type *,' nfte1,nfte2 ',nfte1,nfte2

c                                  boucle 2 canaux
            ndiff=nfte2-nfte1
      do n=nfte1,nfte2,ndiff
        xa=0.1+(n-nfte1)*0.267/ndiff
        xb=xa+0.267
        call pgvport(xa,xb,.1,.8)
        call pgwindow(1.,yl,1.,xl)
          write(canal(1:2),'(i2)')n
c          call pgsch(0.6)
        call pgbox('bcnts',20.,2,'bcnts',200.,2)
        call pglabel(' ',' ',canal)
c          call pgsch(1.)

      iliss2=2*iliss+1
      iil1=il1+iliss
      iil2=il2-iliss
      iilm=iil2-iil1+1

      ic=im/2
      jc=jm/2
      zlevel=tab(ic,jc,n)

          do210 i=iil1,iil2   !!!!!--------- boucle en i
          iip=i-iil1+1
          y(iip)=i
c
c      jmarge=4
c      jmarge2=8
c          rec(i)=jmarge+1
c          do j=1,jm-jmarge2
c                                 
              jmax=jc+(jm-jc)*0.25
c             ----                          valeur modifiable
          do j=jc,jmax
           do igd=-1,1
            jj1=jc+(j-jc)*igd
            jj2=jj1+igd
                        z1=0.
                        z2=0.
              do id=-iliss,iliss
              z1=z1+tab(i+id,jj1,n)
              z2=z2+tab(i+id,jj2,n)
              enddo
                    deni=2*iliss+1.
                    z1=z1/deni
                    z2=z2/deni

          z1=tab(i,jj1,n)
          z2=tab(i,jj2,n)
            test=(z1-zlevel)*(z2-zlevel)
                rec(i)=jc
            if(test.le.0.)then
               rec(i)=j
c               rec(i)=jc
               x(iip)=j
               goto210
             endif
            enddo
           enddo
      if(iip.ne.1)rec(i)=rec(i-1)
210   continue
c                        -----------      fin boucle en i
c        den=float(nfte2-nfte1+1)
c      do i=1,im
c      rec(i)=rec(i)/den
c      enddo
c
c       type*,' cmf: n, rec(i)',n,(rec(i),i=iil1,iil2,10)
        im8=im-8
        im16=im-16
        jm2=jm/2
      iic=0.5*(il1+il2)
      iic1=iic+1
      type *,il1,il2,iil1,iil2
      if(curv.eq.0)then
c        ---------
      call ZERO3(iil1,iic,iic1,iil2,0,0,0,rec,1,im,0,0,
     1           ar,br,cr,item,nrjm,nrej,y1,y2,ifin)
      cr(1)=br
      cr(2)=ar
      cr(3)=0.
      else
c
      call ZERO3(iil1,0,0,iil2,3000,5,3,rec,1,im,0,0,
     1           ar,br,cr,item,nrjm,nrej,y1,y2,ifin)     
      endif

220   continue
      do i=1,im
      y(i)=i
      x(i)=cr(1)+i*(cr(2)+i*cr(3))
c      if(n.lt.nc)ylinea(i)=x(i)
      yline(i,n,nw)=x(i)
c      if(n.gt.nc)ylineb(i)=x(i)
      enddo
      call pgline(im,x,y)
c          ------
      type *,(x(i),i=1,im,20)
c
      enddo
C                                    fin boucle canaux
      den=2.
      do i=1,im
        piv=0.
      do n=nfte1,nfte2,ndiff
        piv=piv+yline(i,n,nw)
      enddo
      yline(i,nc,nw)=piv/den
      enddo
c                                  dessin A1 pour option 2
      endif
c                                    fin option inclin=2
c
c---------------------------------
c                                    calcul de jz,ja,jb, jt pour toutes options
      ic=(1+im)/2
      jz100=100.*yline(ic,nc,nw)
      jt100=100.*trj(nw)
      ja100=yline(1,nc,nw)*100.
      jb100=yline(im,nc,nw)*100.
      type *, jt100,ja100,jb100,jz100
240   format(/' jt100,ja100,jb100,jz100 = '4i8) 
      write(3,240)jt100,ja100,jb100,jz100
250   continue
      type*,' cmf: kdecal ',kdecal
c***************
c dessin                                        titre et A2
      write(nom(1:1),'(a1)')'f'
        write(titre(1:20),'(a)')fna(nw) !!!!!!!!!!!!!!!!
        write(titre(21:25),'(a5)')'     ' !!!!!!!!!!!!!
        write(titre(26:31),'(i6)')jt100 !!!!!!!!!!!!!!!
        write(titre(32:37),'(i6)')ja100 !!!!!!!!!!!!!!!
        write(titre(38:43),'(i6)')jb100 !!!!!!!!!!!!!!!
        write(titre(44:49),'(i6)')jz100 !!!!!!!!!!!!!!
        write(titre(50:55),'(a)')'      '
          call date(jdate)
        write(titre(56:64),'(a)')jdate
            type *,' titre ',titre
      call pgvport(0.1,0.9,0.1,0.85)
c     call pgsch(1.)
      call pglabel(' ',' ',titre)
c     call pgsch(0.7)
c
      call pgvport(0.632,0.90,.1,.8)
      call pgwindow(1.,yl,1.,xl)
c        call pgsch(0.8)
      call pgbox('bcnts',20.,2,'bcnts',200.,2)
c        call pgsch(1.)
c        do k=1,3
c        piv=trj(nw)*(k-2)
c          do i=1,im
c          x(i)=yline(i,nc,nw)+piv
c          y(i)=i
c          enddo
c        call pgline(im,x,y)
c        enddo
c***************
      call pgsls(2)
      piv=0.5*yl-yline(ic,nc,nw)+1.
      x(1)=piv+yline(iil1,nc,nw)-trj(nw)/2.
      x(2)=piv+yline(iil1,nc,nw)+trj(nw)/2.
      y(1)=iil1
      y(2)=iil1
      call pgline(2,x,y)
      x(1)=piv+yline(iil2,nc,nw)-trj(nw)/2.
      x(2)=piv+yline(iil2,nc,nw)+trj(nw)/2.
      y(1)=iil2
      y(2)=iil2
      call pgline(2,x,y)
c
      do i=iil1,iil2
      iip=i-iil1+1
      x(iip)=yline(i,nc,nw)+piv-trj(nw)/2.
      y(iip)=i
      enddo
      call pgline(iilm,x,y)
c
      do i=iil1,iil2
      x(i-iil1+1)=yline(i,nc,nw)+piv+trj(nw)/2.
      enddo
      call pgline(iilm,x,y)
      call pgsls(1)

c*******************************
c     Calcul des coefficients par la moyenne sur le canal
      cf(n1,1) = 1.0
      cf(n1,2) = 1.0
      idec = trj(nw)+0.5
c                             calcul sur trj et trj+1
      do ib=1,2
         type *,ib
         id= idec + ib - 1
      do n= n1+1, n2
         som1=0
         do j= 9, jm-id -8   
            do i= 9, im-8
               som1= som1 + tab(i,j,n-1)
            enddo
         enddo
         som2=0
         do j= 9+id, jm-8
            do i= 9, im-8
               som2= som2 + tab(i,j,n)
            enddo
         enddo
         cn= float(som1)/float(som2)
         cf(n,ib) = cf(n-1,ib) * cn
      enddo
      enddo

      piv1=cf(nc,1)
      piv2=cf(nc,2)
      rap= trj(nw) - idec
c                                 interpolation entre trj et trj+1
      do n=n1,n2
      coef(n)=cf(n,1)/piv1 + ((cf(n,2)/piv2-cf(n,1)/piv1) * rap)
      enddo
c***************
      type *,'calcul de prof100'      !!!!!!!!!!!!!!!!!!!!!!!
c calcul de PROF100(101,11)
c           ---------------  
      profmax=0.
c                             boucle n
        do n=n1,n2
c                                 boucle l
        do l=1,101
        piv=-trj(nw)*(l-51)/100.+1.+0.5*yl
c           ------------------------------
        sum=0.
c                           moyennes sur i
          do i=9,im8
          yj=yline(i,nc,nw)+piv-yline(ic,nc,nw)

c            j= 1+ (jm-1)/2  -(l-51)*trj/100
c bord canal j= 1+ (jm-1)/2  -49    *trj/100

          j=yj
          alp=yj-j
          sum=sum+tab(i,j,n)*(1.-alp)+tab(i,j+1,n)*alp
c                       -                   ---
          enddo
        prof100(l,n)=sum/float(im16)
        enddo
      enddo

c******************
       type *,'calcul de coef(11)'   !!!!!!!!!!!!!!!!!!!!!!
c calcul de COEF(11) pour raccord entre bandes
c      coef(n1)=1.
c      do n=n1+1,n2
c      type *,'n: ',n                     !!!!!!!!!!!!!!!!!!
c      type *,'prof100(101,n-1): ',prof100(101,n-1) !!!!!!!!!!!!!!
c      type *,'prof100(1,n): ',prof100(1,n) !!!!!!!!!!!!!!!!!!!!!!
c      coef(n)=coef(n-1)*prof100(101,n-1)/prof100(1,n)
c      enddo
c
c      piv=coef(nc)
c      do n=n1,n2
c      coef(n)=coef(n)/piv
c      enddo
c****************
      type*,' cmf: kdecal ',kdecal

      type *,'calcul de profm' !!!!!!!!!!!!!!!
c calcul de PROFM(2000,8)
c -----------------------
      do n=n1,n2
        do l=1,101
        k=1000+(n-ncm)*100+l-51
c       ---------------------- 
c                 rapporte au centre de ncm (=(n1+n2)/2)
        profm(k,nw)=prof100(l,n)*coef(n)
        enddo
      enddo
c
c lissage, extrapolation et dessin de PROFM(2000)
      i1=1000+(n1-ncm)*100-50
      i2=1000+(n2-ncm)*100+50
        do i=1,2000
        x(i)=profm(i,nw)
        enddo
      call parafit(x,i1,i2,lissprofm,y)
c          -------
      piv=(y(i1+100)-y(i1))/100.
        do i=1,i1
        y(i)=y(i1)+(i-i1)*piv
        enddo
      piv=(y(i2)-y(i2-100))/100.
        do i=i2,2000
        y(i)=y(i2)+(i-i2)*piv
        enddo
c
        do i=1,2000
          profm(i,nw)=y(i)
        enddo
c------------------------------- 
c                                 dessin  B
      profmax=0.
c
      do n=n1,n2
      do l=1,101
      profmax=amax1(profmax,prof100(l,n))
      enddo
      enddo
c
      do i=i1,i2
      profmax=amax1(profmax,profm(i,nw))
      enddo
c
      call pgadvance
c          ---------
      call pgvport(0.1,0.9,.1,.8)
      xa=n1-2.
      xb=n2+1.
      dmax=profmax*1.2
      call pgwindow(xa,xb,0.,dmax)
      piv=500.
      npiv=5
        if(dmax.lt.500.)then
        piv=100.
        npiv=2
c       call pgsch(0.4)
        endif
      call pgbox('bcgint',1.,0,'bcints',piv,npiv)
c     call pgsch(0.7)
        do n=n1,n2
          do l=1,101
          x(l)=n-1.+(l-1.)/100.
          y(l)=prof100(l,n)
          enddo
        call pgline(101,x,y)
        enddo
c
      xa=i1-100.
      xb=i2+100.
      call pgwindow(xa,xb,0.,dmax)
      call pgslw(2)

        do i=1,2000   
        x(i)=i
        y(i)=profm(i,nw)
        enddo

          profmin=100000.
        do i=700,1300
        profmin=amin1(profmin,y(i))
         if(profmin.eq.y(i))kmin(nw)=i   !!     minimum profm
        enddo
      call pgsls(2)
      call pgline(2000,x,y)

      call pgsls(1)
      call pgslw(1)
260   continue
c                      fin boucle nw
c  On a seulement calcule profm et yline pour chaque nw

c****************************************************************
      type*,' cmf: kdecal ',kdecal
c                      moyennes des trj, profm et yline
c                      --------------------------------

c                       attention le signe de trj depend de invern

      trjm=0.
      do i=1,2000
      profmm(i)=0.
      enddo
      do i=1,1024
      ylinem(i,ncm)=0.
      enddo
c                             moyenne trj   
         den=nw2-nw1+1.
                                         type*,' den ',den
      do  nw=nw1,nw2
      trjm=trjm+trj(nw)/den
      enddo 
c                             moyenne-somme profm

c  le point k=1000 de profmm correspond a la longueur d'onde au point 
c  ic  du centre du canal ncm, dans la fenetre nw1
c  --     ------          ---                  ---
c
c autre i                                          deltak= ylinem(i)*100/trj  
c         autre j                                          -(j-jcentre)*100/trj
c                   autre n                                (n-ncm)*100
c                                       autre  nw          kdecal(nw)

      do k=1,2000
        piv=1.
      do nw=nw1,nw2
          kp=k+kdecal(nw)
c         ---------------
c                        kdecal = decalage theorique en l
c                                             moyenne geom valable pour 2!!!!
          if(kp.ge.1.and.kp.le.2000)
c     1      profmm(k)=profmm(k)+profm(kp,nw)/den
     1      piv=piv*profm(kp,nw)
      enddo
         if(nw1.eq.nw2)profmm(k)=piv
         if(nw2-nw1.ge.1)profmm(k)=sqrt(piv)
         if(nw2-nw1.eq.3)profmm(k)=sqrt(profmm(k))
         pm(k)=profmm(k)
      enddo

c                        symetrie
      if(isym.ne.0)then
      profmin=100000.
        do i=700,1300
        profmin=amin1(profmm(i),profmin)
        if(profmm(i).eq.profmin)imin=i
c          ------               ----
        enddo
      write(3,*)'symetrie: imin, profmin ',imin,profmin
 
      do i=1,2000
      ip=2*imin-i
      if(ip.le.2000)x(i)=0.5*(profmm(i)+profmm(ip))
      enddo

      do i=1,2000
      profmm(i)=x(i)
      enddo
      endif
c                             moyenne yline (recentree sur zero)
c                                                          ----                   
      do i=1,im
      do nw=nw1,nw2
        nc=ncw(nw)
        ylinem(i,ncm)=ylinem(i,ncm)+
     1              (yline(i,nc,nw)-yline(ic,nc,nw))/den
      enddo
      enddo

      jt100m=100.*trjm             !!!  pour sortie
      jz100m=100.*ylinem(ic,ncm)
      ja100m=100.*ylinem(1,ncm)
      jb100m=100.*ylinem(im,ncm) 

      type*,' cmf: jt100m,jz100m,ja100m,jb100m ',
     1             jt100m,jz100m,ja100m,jb100m

c------------------------
c                              Dessin C
      call pgadvance
      call pgvport(0.3,0.7,0.1,0.8)
      call pgwindow(-150.,+150.,700.,1300.)
      call pgbox('bcnts',100.,4,'bcnts',100.,4)
      call pglabel('kdecal','kmin','Entrelac')

        call pgsch(2.)
      do nw=nw1,nw2
      xp=kdecal(nw)
      yp=kmin(nw)
          type*,' kdecal,kmin ',xp,yp
      call pgpoint(1,xp,yp,2)
c          -------                      Graphique Entrelac
      enddo
        call pgsch(1.)

c        type*,' cmf: apres dessin entrelac'

c                                        bissectrice
        xgr=0.
        ygr=0.
      do nw=nw1,nw2
      xgr=xgr+kdecal(nw)/den
      ygr=ygr+kmin(nw)/den
      enddo

        call pgsls(2)
      do ip=1,2000
         i=ip-1000
      x(ip)=i
      y(ip)=ygr+i-xgr
      enddo
        call pgline(2000,x,y)
c            ------
        call pgsls(1)

c-------------------------
c                              Dessin D

      xl=(n2-n1+3)*100.
      
      xk1=1000-xl/2.
      xk2=1000+xl/2.
      k1=xk1
      k2=xk2

           pmin=100000.
           pmax=0.
        do k=k1,k2
        piv=profmm(k)
        pmax=amax1(pmax,piv)
        pmin=amin1(pmin,piv)
            if(piv.eq.pmin)kminm=k
c                          -----
          do nw=nw1,nw2
          piv=profm(k,nw)
          pmax=amax1(pmax,piv)
          enddo
        enddo

        piv=500.
        npiv=5
        if(pmax.lt.500)then
          piv=100.
          npiv=5
        endif

        call pgadvance
      call pgvport(.1,.9,.1,.8)
      call pgwindow(xk1,xk2,0.,pmax)
      call pgbox('bcnts',100.,4,'bcnts',piv,npiv)
c                                                   courbes profm
      do nw=nw1,nw2
      call pgsls(nw+1)
 
        do k=1,2000
        x(k)=k
          kp=k+kdecal(nw)
c         ---------------
          if(kp.ge.1.and.kp.le.2000)
     1      y(k)=profm(kp,nw)
        enddo

      call pgline(2000,x,y)
c          ------
      enddo
c                                                  courbe moyenne profmm
      call pgsls(1)
      do k=1,2000
      y(k)=profmm(k)
      enddo

      call pgline(2000,x,y)
c          ------

c******************               Boucle nw pour creation fichiers F
c                                 ----------------------------------

c        type*,' cmf: nw1.nw2 ',nw1,nw2

      do500 nw=nw1,nw2

c                                              calcul de FL(im,jm,nm)
c               Fl est calcule pour
c               reproduire le profil moyen profmm
c               quand la longueur d'onde du point i,j,n
c               est comptee par reference a ncm et kdecal(nw)
c               (c'est a dire dans les conditions de cmd.f)

          type *,' cmf: trjm,cna(nw),kdecal(nw)',trjm,cna(nw),kdecal(nw)

      iuc=26
      call openold20(cna(nw),sundec,26)
c          ---------
      rewind(iuc)
      read(iuc)kz

c------------------
      do400  n=1,nm
        write(3,*)' nm,im,jm ',nm,im,jm
c      type *,' nw, n ',nw,n

        do j=1,jm
            
        read(iuc)(rec(i),i=1,im)
c                 ------

        yk=1000+(n-ncm)*100-(j-(1+jm)*0.5)*100./trjm-kdecal(nw)
c       --       -           -                      -----------
        if(j.eq.jm/2)write(9,*)' cmf: n,ncm,kdecal(nw) ',
     1                                n,ncm,kdecal(nw)

          do i=1,im
c             if(trj(nw).eq.0.)return
          yyk=yk+(ylinem(i,ncm)-ylinem(ic,ncm))*100./trjm
c         --- --         -             --            ----
          k=yyk
          alp=yyk-k
          piv=profmm(k)*(1.-alp)+profmm(k+1)*alp
c             ---------

          if(i.eq.100.and.j.eq.10)write(3,*)
     1       ' cmf: nw,n, i,j ',
     2              nw,n,i,j

          piv=10000.*rec(i)/piv
c             -----------------

      if(i.eq.187.and.j.eq.25)write(9,*)
     1  ' cmf calcule fl: nw,n, k,rec(187),profmm(k) ',
     2                    nw,n, k,rec(i),profmm(k)

          if(piv.gt.32000.)piv=32000.
          if(piv.lt.0.)piv=0.
          rec(i)=piv+0.5
c         ------
          enddo

        if(smoothi.ne.0)then
        call zero3(9,im2,im2p,im8,0,0,0,rec,1,im,0,0,
c            -----
     1             ar,br,cr,item,nrjm,nrej,y1,y2,ifin)
            if(smoothi.eq.1)then
            y1=0.5*(y1+y2)
            y2=y1
            endif
          do i=1,im
          piv=(i-1)/float(im-1)
          rec(i)=y1*(1.-piv)+y2*piv+0.5
          enddo
        endif

          do i=1,im
          fl(i,j,n)=rec(i)
c         ---------
          enddo
        enddo
c                                      fin des 2 boucles j,i
      do300 i=1,im
        if(smoothj.eq.0)goto300
          do j=1,jm
          rec(j)=fl(i,j,n)
          enddo
        jc=jm/2
        jc1=jc+1
        jm8=jm-8
        call zero3(9,jc,jc1,jm8,0,0,0,rec,1,jm,0,0,
c            -----
     1             ar,br,cr,item,nrjm,nrej,y1,y2,ifin)
        if(smoothj.eq.1)then
        y1=0.5*(y1+y2)
        y2=y1
        endif
c
            do j=1,jm
            piv=(j-1)/float(jm-1)
            fl(i,j,n)=y1*(1.-piv)+y2*piv
            enddo

300     continue
c                                   fin des 2 boucles i,j
400     enddo
c                                         boucle n
      close(unit=26)

c                           dessin de FL(i,jm2,n) et FL(im2,j,n)
      xa=n1-1
      xb=n2
      xl=im-1
      yl=jm-1
c
      call pgadvance
      call pgvport(0.1,0.9,0.1,0.8)
      call pgwindow(xa,xb,0.,2.)
      call pgbox('bcgint',1.,0,'bcints',0.5,5)
      call pglabel('F(i,jm/2,n)',' ',' ')
        do n=n1,n2
          do i=1,im
          x(i)=n-n1+(i-1)/xl
          y(i)=0.0001*fl(i,jm2,n)
          enddo
        call pgline(im,x,y)
c            ------
        enddo
c
      call pgadvance
      call pgvport(0.1,0.9,0.1,0.8)
      call pgwindow(xa,xb,0.,2.)
      call pgbox('bcgint',1.,0,'bcints',0.5,5)
      call pglabel('F(im/2,j,n)',' ',' ')
      im2=im/2
        do n=n1,n2
          do j=1,jm
          x(j)=n-n1+(j-1)/yl
          y(j)=0.0001*fl(im2,j,n)
          enddo
c             type *,(y(j),j=1,jm)
        call pgline(jm,x,y)
c            ------
        enddo

C*******************
c                                       ecriture de FL(im,jm,11)  
c                                       ------------------------
      call system('rm '//fna(nw))
      iuf=27
      call opennew20(fna(nw),sundec,iuf)
c                    -------
      type *,'dans cmf.f avant write(iuf) day,time,lbda',day,time,lbda!!!!!
      write(iuf)6,im,jm,nm,2000,1,1,10,day,time,iline,lbda,ncm,
     1         jt100m,ja100m,jb100m,jz100m,kdecal(nw),kminm,
     2         (kz(i),i=20,512) !!!!!!!!!!!!
c
      do n=1,nm
      do j=1,jm
      write(iuf)(fl(i,j,n),i=1,im) !!!!!!!!!!!!!!!!!!
      enddo
      enddo
c
        do i=1,2000
        pm(i)=profmm(i)+0.5
        if(pm(i).gt.32000.or.pm(i).lt.0)pm(i)=32000   !!!!!!!!!!!!!!!
        enddo
      write(iuf)pm !!!!!!!!!!!!!!!!!!!
980   format(10i8)
c      type 980,(pm(i),i=1,2000,10)
c      write(3,980)(pm(i),i=1,2000,10)
c
      close(unit=27)
500   continue
      return
      end

c***************************************************
