h47037
s 00004/00004/00836
d D 1.2 99/06/11 13:22:28 bonmarti 2 1
c Version du 11 Juin 99
e
s 00840/00000/00000
d D 1.1 99/06/11 12:52:50 bonmarti 1 0
c date and time created 99/06/11 12:52:50 by bonmarti
e
u
U
t
T
I 1

c************************************************
c Sous-programme  cmd.f  -  UNIX
c************************************************

      subroutine cmd(cmddat,nw1,nw2,nwinp,npo,cna,dname,iud, iuj,
c                ---                              pour noms  j pour corv=1
     1               corv,   tab,    ides)
c                     in    economie dessin
c   800= max im
c   100= max jm
c    32= nbre max de canaux au total
c     9= nbre max de canaux par fenetre
c     8= nbre max de fenetres

c    32.ge. 1+2*(lmp(1)+lmp(2)+lmp(3)) nbre de tableaux en sortie
c    32.ge. nm (pour chaque polarisation)
c     arrays "3"    velocity unit = 10 m/s

      integer*4 lmp(3),lbd1(3),lbdstep(3),corv,
     1          day,time,dayf,timef,center,sum,dlbd,dlb(32)
c                                                       lm<=32
                   
      dimension p(32),xp(125),v(32), ylinef(800),kz(512),
c                 nm     4*nm-3  lm<=32
     1          xdes(2),ydes(2)

      integer*2          tab(400,60,32)
c                           nm puis lm
c                        entree    sortie

      character comm*80,cmddat*6,jjname*20,dname*20,
     1          dnanfp*20,cna*20,titre*28
      dimension dnanfp(2000,2),cna(26)

      sqpi=sqrt(3.1415926535)

c Lecture sur .dat

1     format(a)
2     format(2x,a)
3     format(1x,a)
4     format(10i4)
5     format(60i1)
6     format(2x,60i1)
10    format(10f6.2)
12    format(14i6)  !!!!!!!!!!!!!!!
14    format(2x,12i6)
16    format(' Donnees sur .dat')
22    format(/)
24    format(1h1)

      type*,' entree cmd: ides ',ides


      rewind(8)
      type *,'cmddat:',cmddat
      call sname(8,cmddat)
D 2
      type *,'apres sname'
E 2

      read(8,1)comm
      write(3,2)comm
      read(8,12)inveri,inverj,inverl,i1,i2,j1,j2,dlbd, 
     1           nlis,ilis,jlis,idelc
      write(3,14)inveri,inverj,inverl,i1,i2,j1,j2,dlbd,
     1           nlis,ilis,jlis,idelc
c                                        dlbd = entre canaux finaux
c                                                            ------
c                si corv=0, pas de correction de lambdas en entree
c                              (en particulier si longw et lat.ne.0)
c                           on cree un fichier l
c                       =1, on lit le fichier l existant 
c                               on corrige le zero des lambdas
c                           on ne cree pas de fichier l

      read(8,1)comm
      write(3,2)comm
      read(8,12)center,sum,longw,lat
      write(3,14)center,sum,longw,lat
      read(8,1)comm
      write(3,2)comm
        do m=1,3
        read(8,12)lmp(m),lbd1(m),lbdstep(m)
        write(3,14)lmp(m),lbd1(m),lbdstep(m)
        enddo

      lm=center+sum+2*(lmp(1)+lmp(2)+lmp(3))
      If (lmp(3).gt.1) lm = lm +2
c        lm=lm*npol !!!!!!!!!!!!
      if(center.ne.0)dlb(1)=0
      if(sum.ne.0)dlb(1+center)=10000
      kpiv=center+sum
        do m=1,3
          if(lmp(m).ne.0)then
            do mm=1,lmp(m)
            kpiv=kpiv+1
            dlb(kpiv)=(lbd1(m)+lbdstep(m)*(mm-1))*float(dlbd)/4.
            kpiv=kpiv+1
            dlb(kpiv)=-dlb(kpiv-1)
            enddo
          endif
        enddo

        write(3,22)
        piv=3.1415927/180.
          vrot=0.
          if(longw.ne.0.and.lat.ne.0)then
        xlongw=float(longw)*piv
        xlat=float(lat)*piv
        slong=sin(xlongw)
        clat=cos(xlat)
        slat=sin(xlat)
        slat2=slat**2
        vrot=-(2.919-slat2*(0.312+0.212*slat2))*0.696*slong*clat
        vrot=100.*vrot
          endif
c                         dam/s,  + pour blueshift
c---------------------
c                                  ouvertures c 
c                            (d et eventuellement j deja ouverts)
        write(3,*)' cna ',cna

       do nw=nw1,nw2

        iuc=30+nw
      
        call openold20(cna(nw),sundec,iuc)
c            ---------
       enddo
                   
c -------------               Lecture header premier fichier c
c                                pour im,jm,nm,yline, kdec
c                                pour header d 
      iucw=30+nw1

      rewind(iucw)
      read(iucw)na3,im,jm,nmp,imm,nmm,jmm,nxx,nyy,
     1             nd,day,time,iline,lbda,ncmp,
     2             jt100,ja100,jb100,jz100,kdec,           kminm
c                  ----------------------- kdecal pour nw1  moyen
      type *,'na3:',na3,'im:',im,'jm:',jm,'nmp:',nmp
      type *,'imm:',imm,'jmm:',jmm,'nmm:',nmm
      type *,'nxx:',nxx,'nyy:',nyy,'nd:',nd
      type *,'day:',day,'time:',time,'iline:',iline
      type *,'apres read iucw lbda :',lbda

c     nmp= nombre de canaux initiaux
      nm=nmp*nwinp
      ncm=1+(ncmp-1)*nwinp  !!  canal final (sur 32) de reference pour yline
c     ---
      nnm=4*nm-3  !!   nombre de canaux apres interpolation (sur 125)
      ltrj=jt100

D 2

E 2
c      yminm=
c                                      au centre du canal ncm

      nclim=1

c-----------------     debut du header fichier d
      do i=1,512
      kz(i)=0
      enddo

      nd1=4
      nd2=9+lm
c      iz1=11+lm+1+11 !!!!!!!!!
       iz1=18+lm+1
      if(i1.eq.0)i1=1
      if(i2.eq.0)i2=im
      if(j1.eq.0)j1=1
      if(j2.eq.0)j2=jm
        iim=i2-i1+1
        jjm=j2-j1+1

30    format(/' Header fichier d')
      write(iud)nd1,iim,jjm,lm,npo,
     1          nd2,day,time,iline,lbda,ncm,ltrj,(dlb(l),l=1,lm),
c                                           ----
     2          center,sum,lmp,lbdl,lbdstep,
     3          (kz(i),i=iz1,512) !!!!!!!!!!!!!!!
      write(3,30)
      write(3,14)nd1,iim,jjm,lm,npo,
     1           nd2,day,time,iline,lbda,ncm,ltrj,(dlb(l),l=1,lm),
     2           center,sum,lmp,lbdl,lbdstep !!!!!
c----------------------------------------------------
c      if(corv.eq.0)then
c header fichier l
c      rewind(iuj)
c      izl=4+4+lm+2+3*3+1
c      write(iul)3,npol,nxx,nyy,
c     1          nd2,day,time,iline,lbda,(dlb(l),l=1,lm),
c     2          center,sum,lmp,lbdl,lbdstep
c      write(iul)(kz(i),i=iz1,512)
c      endif

c------------------------------         calcul ylinef
c             average line-center = n=ncf and
c                                   j=(ja100+jz100)/100 in i=1
c                                            jz100/100     i=(1+imf)/2
c                                     (jb100+jz100)/100    i=imf

      imf=im
      jmf=jm

      trj=float(jt100)/100.

      xic=(imf+1)/2
         iic=xic

      do i=1,imf
      xpiv=(i-xic)/(xic-1.)
      piv=jz100+xpiv*0.5*(jb100-ja100+
     1                   (ja100+jb100-2.*jz100)*xpiv)
      ylinef(i)=piv/100.
c     ---------          pour le canal final (sur 32) ncm
c                                                     ---
      enddo

      write(3,*)'cmd: xic,ja100,jb100,jz100,trj ',
     1            xic,ja100,jb100,jz100,trj

48    format(' xic,ylinef(iic) ',2f8.1)
      write(3,48)xic,ylinef(iic)
     
c-------------                           vmoy
              vmoy=0.
        write(3,*)' corv ',corv
      if(corv.ne.0)then
c        ---------
c                lecture du fichier j pour correction des lambdas
c          ce fichier, cree par dme.f, contient aussi les resultats
c          de correlations, que l'on doit sauter apres chaque lecture 
c          en cas de boucle Y

      
52    format(' Lecture fichier j pour correc lambdas: vvmoy(lpo,ny)')
      write(3,52)
      rewind(iuj)
        read(iuj,1)jjname
        type *,jjname
        write(3,1)jjname 
      read(iuj,1)comm

108       format(f12.1)
        read(iuj,108)vmoy
c          vmoy(lpo)=vmoy(lpo)+vvmoy(lpo,ly)/float(lymax)
  
110   format(' vmoy='f12.1)
      write(3,110)vmoy
      endif
c-------------------------------------
c                                                 ylinef
c                                            (calcul necessitant ltrj)

c           decaly=-(float(ltrj)/float(ltrjconver))*(nc-ncc)/100. !!!!

c           ylinepo(i,npo)=ylinepo(i,npo)+(ylinef(i,nw)+decaly)/pivd
c
c         Avant de faire la moyenne des ylinef des diveses fenetres nw,
c         on ramene ylinef a sa position dans le canal central final.
c         L'ecart a rattraper en unite de canal final est
c             nc-ncc   ou ncc est le rang n du couple (nwl,ncan)
c         ex: deux fenetres entrelacees:
c         nwl     121212121212121212
c         ncan    112233445566778899
c                          |          ncc=10 pour le canal 5 de la fenetre 2
c                
c         ex: une fenetre:
c              111111111
c              123456789
c                  |               ncc=5 pour le canal 5 de la fenetre 1
c
c         Comme le dlbd initial ne figure pas dans les headers des fichiers
c         b ni f, on calcule le facteur de conversion des ltrj par nwl:
c                121212121212121212
c                | |                  ltrjconver=2
c
c                111111111
c                ||                   ltrjconver=1

c           enddo
c         enddo
c59       format(' npo,vmoy,nw1,nw2,ltrjpo,',
c     1       'ylinepo(1,npo),ylinepo(imf,npo),ltrjconver,decaly'/
c     2            i4,f6.2,2i4,i6,2f8.2,i4,f8.2)
c         write(3,59)npo,vmoy(npo),nw1,nw2,ltrjpo(npo),
c     1              ylinepo(1,npo),ylinepo(imf,npo),
c     2              ltrjconver,decaly


c--------------------                   chargement tab
c  saut des headers

      nw1p=nw2
      nw2p=nw1


      do nw=nw1,nw2
        write(3,*)' cmd: nw ',nw
       
      iucw=30+nw
        rewind(iucw)
      read(iucw)kz
      enddo      
        write(3,*)' cmd: nw,kz ',nw,(kz(i),i=1,20)
c--------------  boucle canaux/fenetres
        ksig=1
        if(inverl.ne.0)ksig=-1

        nt=0       !!  canal nt

      do200 n=1,nmp
c                         entrelac
      do 200 nw=nw1p,nw2p,-1                 !!!!!!!!!!!!   invern.....
        nt=nt+1
            ymoy=-(vrot+vmoy)*ksig*float(lbda)*nwinp/(7500.*dlbd)
c                                       unite= 1/4 de canal final
c                                            (dlbd=resolution finale)
      write(3,*)' cmd: vmoy,ksig,lbda,ymoy ',vmoy,ksig,lbda,ymoy

      iucw=30+nw

      do j=1,jm
c         type*,'                                cmd    do 200 nw,j ',nw,j    
      read(iucw)(tab(i,j,nt),i=1,im)
c               -----------
      enddo
D 2
c                  type*,' cmd: do 200: nw, tab(50,50) ',
c     1                                 nw, tab(50,50,nw)
E 2
I 2
                  type*,' cmd: do 200: nw, tab(50,10) ',
     1                                 nw, tab(50,10,nw)
E 2
200   continue

      call ijlis(tab,im,jm,nmp,ilis,jlis)
c          -----      

c        write(3,*)'  fin boucle fenetres'

c-------------------  calcul des profils ou bissecteurs  pour npo

c                                          boucles j,i
      type *,'boucle 370 j1 :',j1,'j2 :',j2
      type *,'boucle 360 i1 :',i1,'i2 :',i2

        write(titre(1:20),'(a)')dname
c           write(3,*)' titre ',titre
      do370 j=j1,j2
c                                      boucle j1,j2
      do360 i=i1,i2
c                               boucle i1,i2

        itest=0
        if((i.eq.i1.or.i.eq.(i1+i2)/2.or.i.eq.i2).and.
     1     (j.eq.j1.or.j.eq.(j1+j2)/2.or.j.eq.j2))itest=1
        if(ides.eq.0)itest=0

        if(itest.eq.1)then
          write(titre(21:24),'(i4)')i
          write(titre(25:29),'(i4)')j
        endif
        call pgsch(2.)

        l1=1

        do l=1,lm
        v(l)=32000.
        enddo

        do n=1,nm
        p(n)=tab(i,j,n)
c       ---
        if(p(n).eq.32000)then
           if(itest.eq.1)call pgadvance 
           goto344
        endif
        enddo

        call spline4(p,1,nm,xp,piv1,piv2,nlis)
c            -------
c          write(3,*)'apres spline4' !!!!!!!!!!!!!!!!
         nnm=4*nm-3
      if(itest.eq.1)then
                 call protest(xp,nnm)
c                     -------
      call pglabel(' ',' ',titre)
        write(3,*),'                                       ',dname
        write(3,*)' protest: ncm,ylinef(i),trj ',
     1                       ncm,ylinef(i),trj
        write(3,*)'          profil p ',(p(nn),nn=1,nm)
      endif

        ksig=1
        if(inverl.ne.0)ksig=-1

c                                centre raie pour j = yyc (en unite canaux 125)
      yyc=4.*(ncm+float(kminm-1000)*nwinp/100.

     1            +nwinp*(j-0.5*(1.+jm)-ylinef(i))/trj)  -3.       

             do l=1,lm
             v(l)=32000.
             enddo
c                                             center,sum
      if(center.eq.1)then
      piv=pinter(xp,nnm,yyc,nlis,ymoy)
c                       ---      ----
c     out               in        in

      if(itest.eq.1)then
      write(3,*)' xp ',(xp(nn),nn=1,nnm)
      write(3,*)' nnm,yyc,nlis,ymoy,piv ',nnm,yyc,nlis,ymoy,piv 
      call pgsch(2.)
       pivx=yyc+ymoy
      call pgpoint(1,pivx,piv,2)
c          -------
      endif

         if(piv.ne.32000.)v(l1)=piv*2.+0.5
      endif

      if(sum.eq.1)then
      kpiv=0
        do n=1,nm
        kpiv=kpiv+p(n)+0.5
        enddo
      l=center+l1
      v(l)=0.5+float(kpiv)*2./float(nm)
      endif
c                                              profiles
c      type *,'etape profil' !!!!!!!!!!!
      ksig=1
      if(inverl.ne.0)ksig=-1

      if(lmp(1).ne.0)then !!!!!!lmo -> lmp
        do lb=1,lmp(1)
        lbd=lbd1(1)+(lb-1)*lbdstep(1)
        l=l1+center+sum+2*(lb-1)
        yy=yyc+lbd*ksig
c          ---
             xdes(2)=yy+ymoy
           piv=pinter(xp,nnm,yy,nlis,ymoy)
             ydes(2)=piv
           if(piv.ne.32000.)v(l)=piv*2.
        yy=yyc-lbd*ksig
c          ---
             xdes(1)=yy+ymoy
           piv=pinter(xp,nnm,yy,nlis,ymoy)
             ydes(1)=piv
           if(piv.ne.32000.)v(l+1)=piv*2.

          if(itest.eq.1)then
          call pgline(2,xdes,ydes,1)
c              ------
          endif

        enddo
      endif
c                                    sum and diff +/- lbd
      if(lmp(2).ne.0)then
        do lb=1,lmp(2)
        lbd=(lbd1(2)+(lb-1)*lbdstep(2))
        l=l1+center+sum+2*(lb-1+lmp(1))
        yy=yyc+lbd*ksig
c          ---
            xdes(2)=yy+ymoy
          piv1=pinter(xp,nnm,yy,nlis,ymoy)
            ydes(2)=piv1
        yy=yyc-lbd*ksig
c          ---
             xdes(1)=yy+ymoy
          piv2=pinter(xp,nnm,yy,nlis,ymoy)
             ydes(1)=piv2
        if(piv1.ne.32000..and.piv2.ne.32000.)then
          v(l)=piv1+piv2
          v(l+1)=piv1-piv2
        endif

          if(itest.eq.1)then
          call pgline(2,xdes,ydes,1)
c              ------
          endif

        enddo
      endif
c                                     bissector
c      type *,'etape bissector' !!!!!!!!!!!!!
      if(lmp(3).ne.0)then
        do lb=1,lmp(3)
        lbd=lbd1(3)+(lb-1)*lbdstep(3)
        lbd=2*lbd
        l=l1+center+sum+2*(lb-1+lmp(1)+lmp(2))
        lp=l+1
c              if(i.eq.i1)type*,' l1,center,sum,lb,lmp(1),lmp(2),lp ',
c     1                           l1,center,sum,lb,lmp(1),lmp(2),lp 
          call iv(xp,nnm,lbd,v(l),v(lp),1)
c              --

        if(itest.eq.1)then
         write(3,*)' protest iv: lbd,i,v ',lbd,v(l),v(lp)
          ydes(1)=v(l)
          ydes(2)=v(l)
          xdes(1)=v(lp)-float(lbd)/2.
          xdes(2)=v(lp)+float(lbd)/2.                
        call pgline(2,xdes,ydes,1)
        endif

          if(v(l).ne.32000..and.v(lp).ne.32000.)then
          v(l)=2.*v(l)
             if(v(l).ge.32000.)v(l)=32000.
          v(lp)=v(lp)-yyc
c                    ----
c          type *,'lbda :',lbda !!!!!!!!!!!!
            piv=lbda*nwinp
          v(lp)=-v(lp)*ksig*7500.*dlbd/piv
          v(lp)=v(lp)-vmoy
c                    -----
         if(itest.eq.1)then
         write(3,*)' protest iv: lbd,i,v(dam/s) ',lbd,v(l),v(lp)
         endif
             if(abs(v(lp)).gt.32000.)v(lp)=32000.
          endif
        enddo
      endif

      If (lmp(3).gt.1) then
c                            On rajoute la moyenne I,V
         vi=0.
         vv=0.
         poids=0.
         do lb=1,lmp(3)
            l= l1 + center + sum +2* (lb-1 + lmp(1) + lmp(2))
            If (v(l).ne.32000 . and . v(l+1).ne.32000) then
               poids = poids +1
               vi = vi + v(l)
               vv = vv + v(l+1)
            endif
         enddo
         v(lm-1) = 32000
         v(lm) = 32000
         If (poids.ne.32000) then
            v(lm) = vv/poids
            v(lm-1)= vi/poids
          endif
       endif

c-----
344   do l=1,lm
        kpiv=v(l)
        if(iabs(kpiv).gt.32000)kpiv=32000
      tab(i,j,l)=kpiv
c     ----------
      enddo

360   continue
I 2
       type *, ' j ',j

E 2
c                    fin boucle i1,i2
370   continue
c                        fin boucle j1,j2

c------------------------------ chargement d 
      type *,'chargement d et l'!!!!!!!!!
      type *,'lm:',lm,'inverj:',inverj,'i1:',i1,'i2:',i2
      type *,'j1: ',j1,'j2 : ',j2!!!!!!!!!!!!!!!

      do l=1,lm
381   format(' l='i2)
      write(3,381)l
        do j=j1,j2
        jp=j
        if(inverj.ne.0)jp=jm-j+1
          if(inveri.eq.0)then
          write(iud)(tab(i,jp,l),i=i1,i2)
c                    ----------
          else
          write(iud)(tab(i,jp,l),i=i2,i1,-1)
          endif
        enddo

      enddo
c-----------------------------  close
     
      do nw=nw1,nw2
        iuc=30+nw
        close(unit=iuc)
           if(idelc.eq.1)then
           call system('rm '//cna(nw))
           endif
        enddo

      return
      end
c********************************************************
 
      subroutine protest(y,nm)
c                -------
      dimension x(125),y(125)
      character*23 psname

1     format(' Profil centre du champ premiere image')
      write(3,1)
2     format(10f6.0)
      write(3,2)(y(n),n=1,nm)

      call pgadvance
      call pgvport(0.1,0.9,0.1,0.8)
          ymax=-100000.
        do n=1,nm
        x(n)=n
        ymax=amax1(ymax,y(n))
        enddo

      xnm=nm
      call pgsch(1.)
      call pgwindow(1.,xnm,0.,ymax)
      call pgbox('bcnts',0.,0,'bcnts',0.,0)
      call pglabel('lambda','intensity',' ')
      call pgline(nm,x,y)

      return
      end
c---------------------------------------------
      function PINTER(xp,nnm,yy,liscan,   yrot)
c              ------           inutilise
      dimension xp(1)
c
      nlis=0
c     inhibation de nlis

      pinter=0.
        nlis2=2*nlis+1
        yyy=yy+yrot-nlis-1
      do10 nl=1,nlis2
      y=yyy+nl
      k=y
        if(k.lt.1.or.y.gt.float(nnm))then
        pinter=32000.
        return
        endif
          if(y.eq.float(nnm))then
          pinter=xp(nnm)
          goto10
          endif
      alp=y-k
        if(k.eq.1.or.k.eq.nnm-1)then
        pinter=pinter+xp(k)*(1.-alp)+xp(k+1)*alp
        else
c                    cubique par 4 points 
          a=xp(k-1)
          b=xp(k)
          c=xp(k+1)
          d=xp(k+2)
        x=2.*alp-1.
        pinter=pinter+(-a+9.*b+9.*c-d)/16.+
     1        x*((a-27.*b+27.*c-d)/48.+
     2        x*((a-b-c+d)/16.+
     3        x*((-a+3.*b-3.*c+d)/48.)))
        endif
10      enddo
      pinter=pinter/float(nlis2)          
c
      return
      end
c----------------------
      subroutine spline4(p,n1,n2,xp,a,b,nlis)
c                -------
      dimension p(32),q(32),xp(125),c1(5),c2(5)
      km=4*(n2-n1+1)-3
c                           write(3,*)'nlis=',nlis

      if(nlis.eq.0)then
      do n=n1,n2
      q(n)=p(n)
      enddo
      endif

      if(nlis.eq.1)then
        s1=0.
        s2=0.
      do n=n1,n2,2
      s1=s1+p(n)
      s2=s2+p(n+1)
      enddo
        s=(s1+s2)/2.
        cor1=s/s1
        cor2=s/s2
      do n=n1,n2,2
      q(n)=p(n)*cor1
      q(n+1)=p(n+1)*cor2
      enddo
      endif

      if(nlis.eq.2)then
      q(n1)=.75*p(n1)+.25*p(n1+1)
      q(n2)=.75*p(n2)+.25*p(n2-1)
        do n=n1+1,n2-1
        q(n)=.25*p(n-1)+.5*p(n)+0.25*p(n+1)
        enddo
      endif

      if(nlis.eq.3)then
      q(n1)=.75*p(n1)+.25*p(n1+1)
      q(n2)=.75*p(n2)+.25*p(n2-1)
        do n=n1+1,n2-1
        q(n)=.33333*p(n-1)+.33333*p(n)+.33333*p(n+1)
        enddo
      endif

      call cpoly4(q,n1,c2)
c          ------
      do k=1,9
      x=k-9
      xp(k)=c2(1)+x*(c2(2)+x*(c2(3)+x*(c2(4)+x*c2(5))))
      enddo
c
      do n=n1,n2-5
          do i=1,5
          c1(i)=c2(i)
          enddo
          call cpoly4(q,n+1,c2)
c              ------
        do kp=1,3
        x1=kp
        x2=kp-4
        y1=c1(1)+x1*(c1(2)+x1*(c1(3)+x1*(c1(4)+x1*c1(5))))         
        y2=c2(1)+x2*(c2(2)+x2*(c2(3)+x2*(c2(4)+x2*c2(5))))
        k=kp+4*(n+2)-3
        xp(k)=(x1*y2-x2*y1)/4.
c           if(kp.eq.1)xp(k)=0.75*y1+0.25*y2
c           if(kp.eq.2)xp(k)=0.5*y1+0.5*y2
c           if(kp.eq.3)xp(k)=0.25*y1+0.75*y1
        enddo
        k=4*n+9
        xp(k)=q(n+3)
      enddo
c
c      call cpoly4(q,n2-4,c2)
c          inutile, car deja charge en c2
      do kp=1,8
      k=km+kp-8
      x=kp
      xp(k)=c2(1)+x*(c2(2)+x*(c2(3)+x*(c2(4)+x*c2(5))))
      enddo
c
      a=100000.
      b=-100000.
        do k=1,km
        a=amin1(a,xp(k))
        b=amax1(b,xp(k))
        enddo
      return
      end
c---------------------------------------
      subroutine cpoly4(p,n,cc)
      dimension p(32),cc(5)
      a=p(n)
      b=p(n+1)
      c=p(n+2)
      d=p(n+3)
      e=p(n+4)
c
      cc(1)=c
      cc(2)=(a-8.*b+8.*d-e)/48.
      cc(3)=(-a+16.*b-30.*c+16.*d-e)/384.
      cc(4)=(-a+2.*b-2.*d+e)/768.
      cc(5)=(a-4.*b+6.*c-4.*d+e)/6144.
      return
      end
c------------------------------------------------
      subroutine ijlis(tab,im,jm,nm,ilis,jlis)
c                -----
      integer*2 tab(400,60,32)
      integer*4 mem(400)
      if(ilis.gt.1)then
        do n=1,nm
        do j=1,jm
          do i=1,im
          mem(i)=tab(i,j,n)
          enddo
        call klis(mem,im,ilis)
c            ----
          do i=1,im
          tab(i,j,n)=mem(i)
          enddo
        enddo
        enddo
      endif

      if(jlis.gt.1)then
        do n=1,nm
        do i=1,im
          do j=1,jm
          mem(j)=tab(i,j,n)
          enddo
        call klis(mem,jm,jlis)
c            ----
          do j=1,jm
          tab(i,j,n)=mem(j)
          enddo
        enddo
        enddo
      endif

      return
      end
c--------------------------------------
      subroutine klis(mem,km,lis)
c                ----
      dimension mem(400),memp(500)

      do kd=1,50
      memp(kd)=mem(1)
      memp(km+kd+50)=mem(km)
      enddo

      do k=51,km+50
      memp(k)=mem(k-50)
      enddo

      if(lis.eq.2)then
        do k=1,km
        mem(k)=.25*memp(k+49)+.5*memp(k+50)+.25*memp(k+51)+.5
        enddo
      endif

      if(lis.eq.3)then
        do k=1,km
        mem(k)=.3333*memp(k+49)+.3333*memp(k+50)+.3333*memp(k+51)+.5
        enddo
      endif        

      if(lis.eq.4)then
        do k=1,km
        mem(k)=.125*memp(k+48)+.250*memp(k+49)+
     1         .250*memp(k+50)+
     2         .250*memp(k+51)+.125*memp(k+52)+.5
        enddo
      endif

      return
      end

         
E 1
