      CHARACTER*(*) SOURCE
      PARAMETER(SOURCE='#'//CHAR(13)//CHAR(10)//'# Time domain 15 degree
     & (t'',x'',tau'') migration of zero-offset section '//CHAR(13)//CHA
     &R(10)//'#'//CHAR(13)//CHAR(10)//'#Usage: <Headin T15mig n1=nt= d1=
     &dt=  n2=nx=  d2=dx= vc=/vel='//CHAR(13)//CHAR(10)//'#        trick
     &=.14 taper=no stepdown=1 npad=1 >Headout'//CHAR(13)//CHAR(10)//'#'
     &//CHAR(13)//CHAR(10)//' ')
c Time domain 15 degree (t',x',tau') migration of zero-offset section 
c
cUsage: <Headin T15mig n1=nt= d1=dt=  n2=nx=  d2=dx= vc=/vel=
c        trick=.14 taper=no stepdown=1 npad=1 >Headout
c%
      integer hetch,fetch,getch,tetch,infd,outfd,headfd
      integer fsbrk,loc,locore,input,output,head,idum
      integer lcrpos
      double precision basecore(1)
CMF$LAYOUT BASECORE(:SERIAL)
      common/falloc/basecore
c  the retarded coordinate transform:
c       t'      =       t + tau
c       x'      =       x
c       tau'    =       tau       ( tau = 2*z/v -- psudo-depth )
c
cn1=nt,n2=nx    trace length, number of traces
cd1=dt,d2=dx    sampling intervals
cvc=            constant velocity(v=vc); must specify either vc or vel
cvel=           headfile specifying velocity profile(s) of length nt
c                 when v is v(t) or v(t,x)
ctrick=.14      1/6. trick in more accurate finite difference approximation 
c                 of second x derivative;  trick<0 when dipfilter needed
ctaper=yes      taper side traces to reduce boundary reflections
cnpad=1		minimum number of zero traces to pad on each side to
c		reduce wraparound.
cstepdown=1     increment of downward step
c  
c author:       Zhiming Li, 8/17/84
c
c       Modified:       
c               stew        --- transfered to convex and took out ap code 
c               li 12/6/85  --- transpose data and velocity array       
c               stew 4/18/86 -- vectorized rewrite based on SEP-38 article
c               stew 4/25/86 -- convert interval vels to tau step internally
c               stew 5/6/86  -- hand roll rhs and triply outer loops for speed
c               stew 5/26/86 -- somewhat faster recoding
c               marta 6/17/86---changed rite to npad+1 instead of npad
c               stew&marta 9/30/86---fixed alp(ntau,nx) alloc and declaration
c                                 ---changed min0(ittop,(it-1)/(stepdown-1))
c                                    to min0(ittop,(it-1)/(stepdown-1)-1)
c               marta 10/1/86--- fixed bug that left beta undefined
c               stew 7/27/89 --- For CM2 Fortran 8x
c               stew 8/6/89 --- Pull tridiagonal code inline so I can use
c                               CMF$LAYOUT directives.
c		stew 8/11/89 --- Riemann sheeted input/output section
c                                to greatly reduce communication cost.
c		stew 9/19/89 --- added modpar() and reinitialization
c		stew 9/20/89 --- added frame buffer display of migration
c		stew 9/21/89 --- moved display into migration loop
c		stew 9/26/89 --- reduce stepdown for accuracy when it
c				 won't affect speed.
c
c ----------
c  Keyword:  15-degree time-domain migration
c ----------
      integer nt,nx,stepdown,npad
      integer nxpow2,nxlog2,ntaupow2,nustep
      real eps,dt,dx,beta,vc,clip
      character *10 taper
      call initpar
      call doc( source )
      idum=hetch('Gobbledegook','s','Lets get started.')
      lcrpos=1
      if(.not.(fetch('n1','i',nt).eq.0))goto 23000
         call erexit('need n1:nt')
23000 continue
      if(.not.(fetch('n2','i',nx).eq.0))goto 23002
         call erexit('need n2:nx')
23002 continue
      if(.not.(fetch('stepdown','i',stepdown).eq.0))goto 23004
         stepdown=1
23004 continue
      if(.not.(fetch('d1','r',dt).eq.0))goto 23006
         call erexit('need d1:dt')
23006 continue
      if(.not.(fetch('d2','r',dx).eq.0))goto 23008
         call erexit('need d2:dx')
23008 continue
      if(.not.(getch('beta','r',beta).eq.0))goto 23010
         beta=.14
23010 continue
      call putch('From par: beta = beta','r',beta)
      if(.not.(getch('vc','r',vc).eq.0))goto 23012
         vc=0.0
23012 continue
      call putch('From par: vc = vc','r',vc)
      if(.not.(getch('eps','r',eps).eq.0))goto 23014
         eps=.0001
23014 continue
      call putch('From par: eps = eps','r',eps)
      if(.not.(getch('clip','r',clip).eq.0))goto 23016
         clip=5.0
23016 continue
      call putch('From par: clip = clip','r',clip)
      if(.not.(getch('taper','s',taper).eq.0))goto 23018
         taper='no'
23018 continue
      call putch('From par: taper = taper','s',taper)
      if(.not.(getch('npad','i',npad).eq.0))goto 23020
         npad=1
23020 continue
      call putch('From par: npad = npad','i',npad)
      npad=1
      nx=nx+npad*2
      ntau=(nt-1)/stepdown
      nxpow2=1
      nxlog2=0
c     while
23022 if(.not.(nxpow2.lt. nx))goto 23023
         nxlog2=nxlog2+1
         nxpow2=nxpow2*2
         goto 23022
c     endwhile
23023 continue
      ntaupow2=1
c     while
23024 if(.not.(ntaupow2.lt.ntau))goto 23025
         ntaupow2=ntaupow2*2
         goto 23024
c     endwhile
23025 continue
      nustep = (nt-1)/ntaupow2
      nustep=max(1,nustep)
c     while
23026 if(.not.(ntaupow2 .lt. ((nt-1)/nustep)))goto 23027
         nustep = nustep+1
         goto 23026
c     endwhile
23027 continue
      if(.not.(nustep.gt.0 .and. nustep.lt.stepdown))goto 23028
         stepdown=nustep
         ntau=(nt-1)/stepdown
         call putch('stepdown','i',stepdown)
23028 continue
      lvfe = lcrpos
      lcrpos = lcrpos + (((ntau*nx+nt-ntau))*4+7)/8
      lpfe = lcrpos
      lcrpos = lcrpos + (((nt)*(nx))*4+7)/8
      lfactor1 = lcrpos
      lcrpos = lcrpos + (((nxlog2+1))*4+7)/8
      lfactor2 = lcrpos
      lcrpos = lcrpos + (((nxlog2+1))*4+7)/8
      call cm_init
      call cm_set_safety_mode(0)
      locore=fsbrk((lcrpos-1)*8)
      locore=(loc(basecore)-locore)/8
      call parmig(eps,dt,dx,beta,nt,nx,vc,taper,stepdown,ntau,npad,
     & nxpow2,nxlog2,ntaupow2,basecore(lpfe-locore),basecore(lvfe-
     & locore),basecore(lfactor1-locore),basecore(lfactor2-locore),clip)
      stop
      end
      subroutine parmig(eps,dt,dx,beta,nt,nx,vc,taper,stepdown,ntau,
     & npad,nxpow2,nxlog2,ntaupow2,pfe,vfe,factor1,factor2,clip)
      implicit character (a-z)
      integer fetch, hetch, tetch, getch, auxpar
      integer nt,nx,ix,it,itau,infd,outfd,velfd
      integer nxpow2,nxlog2,ntaupow2
      integer stepdown,ntau,ittop,min0,npad,i,j
      real eps,clip,ovc
c       real efac(nxpow2)
      real sc,beta,dt,dx,vc,factor,exp,float
      integer factor1(*),factor2(*)
      integer dispid
      real vfe(ntau,nx),pfe(nt,nx)
      integer itemp,input,output
      integer jt,irem
      integer iada,iadb,iadc
      integer iadislot,iadoslot,iadxedges
      integer iadp1,iadp2,iadp3,iadp4
      integer iadr1,iadr2,iadalp
      integer iadn1,iadn2
      integer iadpsec,iadpin,iadpout
      integer iadpsave, iadtemp
      integer iadnonzero
      integer iadnonzero1
      integer iadnonzero2
      integer iadnonzero3
      integer iadplot
      integer loopcnt,nsheet
      integer lbound,p1addr,vp,geom
      integer modpar
      integer n3
      real tepsilon
      logical rpt,dodata
      character*(*) taper
      integer dscarr(7,2)
       include '/usr/include/cm/paris-configuration-fort.h'
      infd=input()
      outfd=output()
       nsheet=1+(nt-1)/ntaupow2
       if(stepdown.eq.1) nsheet=1
      dispid=CMFB_make_display_id()
      itemp=CMFB_attach_display('D312',dispid)
      call CMFB_initialize_display(dispid,8,1)
      call CMFB_overlay_disable(dispid)
      call CMFB_set_colors_linear(dispid,0,127,
     1          128,128,128,255,255,255)
      call CMFB_set_colors_linear(dispid,128,255,
     1          0,0,0,127,127,127)
      i=max(1,1024/nxpow2)
      j=max(1,1024/(nsheet*ntaupow2))
      call CMFB_set_pan_and_zoom(dispid,0,0,i-1,j-1,0)
      n3=0
      dscarr(1,1)=nxpow2
      dscarr(2,1)=4
      i=getch('xweight','i',dscarr(2,1))
      dscarr(3,1)=cm_news_order
      dscarr(4,1)=0
      dscarr(6,1)=0
      dscarr(1,2)=ntaupow2
      dscarr(2,2)=2
      i=getch('tweight','i',dscarr(2,2))
      dscarr(3,2)=cm_news_order
      dscarr(4,2)=0
      dscarr(6,2)=0
       geom=CM_create_detailed_geometry(dscarr,2)
       vp=CM_allocate_vp_set(geom)
       call CM_set_vp_set(vp)
       iada=CM_allocate_stack_field(32)
       iadb=CM_allocate_stack_field(32)
       iadc=CM_allocate_stack_field(32)
       iadp1=CM_allocate_stack_field(32)
       iadp2=CM_allocate_stack_field(32)
       iadp3=CM_allocate_stack_field(32)
       iadp4=CM_allocate_stack_field(32)
       iadr1=CM_allocate_stack_field(32)
       iadn1=CM_allocate_stack_field(32)
       iadr2=CM_allocate_stack_field(32)
       iadn2=CM_allocate_stack_field(32)
       iadalp=CM_allocate_stack_field(32)
       iadpsec=CM_allocate_stack_field(32*nsheet)
       iadpsave=CM_allocate_stack_field(32*nsheet)
       iadplot=CM_allocate_stack_field(8)
       iadislot=CM_allocate_stack_field(1)
       iadoslot=CM_allocate_stack_field(1)
       iadxedges=CM_allocate_stack_field(1)
       iadnonzero=CM_allocate_stack_field(1)
       iadnonzero1=CM_allocate_stack_field(1)
       iadnonzero2=CM_allocate_stack_field(1)
       iadnonzero3=CM_allocate_stack_field(1)
c
c     xedges=.false.
c     xedges(:,1)=.true.
c     xedges(:,nx)=.true.
c
       call CM_set_context
       call CM_my_news_coordinate_1L(iadn2,0,32)
       call CM_u_eq_zero_1L(iadn2,32)
       call CM_store_test_always(iadxedges)
       call CM_u_eq_constant_1L(iadn2,nx-1,32)
       call CM_logior_test_always(iadxedges)
       call CM_store_test_always(iadxedges)
      if(.not.(nx .lt. nxpow2))goto 23030
       call CM_u_ge_constant_1L(iadn2,nx,32)
       call CM_store_test_always(iadnonzero1)
23030 continue
       call CM_u_eq_constant_1L(iadn2,nx-1,32)
       call CM_store_test_always(iadnonzero2)
       call CM_u_eq_zero_1L(iadn2,32)
       call CM_store_test_always(iadnonzero3)
       call CM_my_news_coordinate_1L(iadn1,1,32)
cccc   ix=CM_set_system_leds_mode(CM_leds_throb)
c
      print *,'reading input section'
c read input section
      iadpin=iadpsave
      dodata = .true.
      call CM_set_context
      call getsec(infd,pfe,npad,nt,nx,iadpin,ntaupow2,dodata)
      dodata = .false.
c interactive loop
c     while
23032 if(.not.(0.eq.modpar()))goto 23033
         n3=n3+1
         call putch('n3','i',n3)
      call CMFB_clear(dispid,128,CMFB_green)
c%      call CM_u_move_const_always_1L(iadplot,128,8)
c       do i=1,2048,nxpow2 {
c       do j=1,1024,ntaupow2 {
c%      call CMFB_write_preshuffled_always(dispid,CMFB_green,iadplot,
c%     1     i-1,j-1)
c       } }
c
c set input and output slot pointers
c
c      islot=.false.
c      islot(1+mod(nt,ntaupow2),:)=.true.
c      if(stepdown.gt.1){
c         oslot=islot
c             } else {
c          oslot=.true.
c          oslot(1,:)=.false.
c     }
       call CM_u_eq_constant_1L(iadn1,mod(nt,ntaupow2),32)
       call CM_store_test_always(iadislot)
         if(.not.(stepdown.gt.1))goto 23034
       call CM_u_move_always_1L(iadoslot,iadislot,1)
            goto 23035
c        else
23034       continue
       call CM_u_ne_zero_1L(iadn1,32)
       call CM_store_test_always(iadoslot)
23035    continue
c copy input section to computational input/output grid
      call CM_set_context
         iadtemp=iadpsave
         iadpout=iadpsec
         do 23036 i=1,nsheet 
          call cm_f_move_always_1L(iadpout,iadtemp,23,8)
          iadtemp=cm_add_offset_to_field_id(iadtemp,32)
          iadpout=cm_add_offset_to_field_id(iadpout,32)
23036       continue
         iadpin=iadpsec
         dodata=.false.
         call getsec(infd,pfe,npad,nt,nx,iadpin,ntaupow2,dodata)
         iadpout=iadpin
c
c STILL NEED TO CONVERT THIS
c
c taper sides to reduce boundary artifacts
c
c       if(taper(:1) == 'y') {
c            efac=0.0
c            do ix=1+npad,nx/32 {
c                factor=float(ix-1-nx/32)/float(nx/32)
c                efac(ix)=factor
c                efac(nx+1-ix)=factor
c                }
c            efac=exp(-efac**2)
c	    r1=spread(efac,1,ntaupow2)
c           iadpout=iadpsec
c	    do it=1,nt,ntaupow2 {
c%      if(it.gt.1)iadpout=CM_add_offset_to_field_id(iadpout,32)
c%      call CM_f_multiply_always_2_1L(iadpout,iadr1,23,8)
c	    }
c        }
c read migration velocities
         print *,'setting up velocity field'
         if(.not.(getch('vc','r',vc).eq.0))goto 23038
            vc=vc
23038    continue
         call putch('From par: vc = vc','r',vc)
         if(.not.(getch('beta','r',beta).eq.0))goto 23040
            beta=beta
23040    continue
         call putch('From par: beta = beta','r',beta)
         if(.not.(getch('eps','r',eps).eq.0))goto 23042
            eps=eps
23042    continue
         call putch('From par: eps = eps','r',eps)
         if(.not.(getch('clip','r',clip).eq.0))goto 23044
            clip=clip
23044    continue
         call putch('From par: clip = clip','r',clip)
         ovc=127.0/clip
         sc = - (1./32.)*(dt/dx)**2
c       alp = 0.0
      call CM_f_move_const_always_1L(iadalp,0.0,23,8)
         if(.not.(vc .eq. 0.0))goto 23046
            call getvel(nt,nx,npad,ntau,vfe,stepdown,sc,iadalp)
            goto 23047
c        else
23046       continue
c constant velocity
            sc = sc * stepdown * vc**2
c           alp=sc
      call CM_f_move_const_always_1L(iadalp,sc,23,8)
23047    continue
c
c set up and prefactor tridiagonal systems
c
         print *,'factoring tridiagonal system'
c       a=beta+alp
c       a(:,nx+1:nxpow2)=0.0
c       b=1.0-2.0*a
c       c=a
c       c(:,1)=-1.0
c       a(:,nx)=-1.0
c
       call CM_f_add_const_always_3_1L(iada,iadalp,beta,23,8)
         if(.not.(nx .lt. nxpow2))goto 23048
       call CM_load_context(iadnonzero1)
       call CM_f_move_zero_1L(iada,23,8)
23048    continue
       call CM_set_context
       call CM_f_mult_const_add_const_1L(iadb,iada,-2.0,1.0,23,8)
       call CM_f_move_always_1L(iadc,iada,23,8)
       call CM_load_context(iadnonzero2)
       call CM_f_move_constant_1L(iada,-1.0,23,8)
       call CM_load_context(iadnonzero3)
       call CM_f_move_constant_1L(iadc,-1.0,23,8)
c
c       where(xedges)b=1.0
       call CM_load_context(iadxedges)
       call CM_f_move_constant_1L(iadb,1.0,23,8)
       call CM_set_context
c
c debug:
c      call putcomp(nxpow2,ntaupow2,iada,vfe)
c      call putcomp(nxpow2,ntaupow2,iadb,vfe)
c      call putcomp(nxpow2,ntaupow2,iadc,vfe)
c
         call trifact(iada,iadb,iadc,ntaupow2,nxpow2,factor1,factor2,
     &    nxlog2,eps,loopcnt)
         print *,nxpow2,ntaupow2,nxlog2,loopcnt
       call cm_reset_timer(0)
       call cm_start_timer(0)
c
c simplify rhs tridiagonal coefficients
c
c       alp=beta-alp
       call CM_set_context
       call CM_f_subfrom_const_always_2_1L(iadalp,beta,23,8)
c initialize parallel subdiagonals
         print *,'starting migration nt=',nt,' nx=',nx,' ntau=',ntau,
     &    ' ntaupow2=',ntaupow2
c       p1=0
c       p2=0
c       p3=0
c       p4=0
      call CM_f_move_zero_always_1L(iadp1,23,8)
      call CM_f_move_zero_always_1L(iadp2,23,8)
      call CM_f_move_zero_always_1L(iadp3,23,8)
      call CM_f_move_zero_always_1L(iadp4,23,8)
c
c now for the "parallel" migration
c
c      P1 P4
c      P2 P3
c
c      (1+(beta+alpha)T) (P4+P2) = (1+(beta-alpha)T) (P1+P3)
c
c At top of loop, P1, P2, and P3 are in correct position
c    P1 runs from it to it+ntaupow2-1
c    P2 runs from it+1 to it+ntaupow2
c    P3 runs from it+1 to it+ntaupow2
c Logic:
c    shift input pointer, slot psec(it) into p1(it)
c    compute P4
c    P4 runs from it to it+ntaupow2-1
c    shift output pointer, slot p4(jt) into psec(jt) as needed
c    shift P1 to P2
c    copy  P4 to P1
c    shift P1 to P3
c
         do 23050 it=nt,1+stepdown,-1 
c work from bottom up
c
c
c         islot=cshift(islot,1,1)
       call CM_get_from_news_always_1L(iadislot,iadislot,1,
     1      CM_upward,1)
            if(.not.(mod(it,ntaupow2).eq.0))goto 23052
       iadpin=CM_add_offset_to_field_id(iadpin,-32)
23052       continue
       call CM_load_context(iadislot)
       call CM_f_move_1L(iadp1,iadpin,23,8)
       call CM_set_context
            if(.not.(.true.))goto 23054
c         R1=P1+P3
       call CM_f_add_always_3_1L(iadr1,iadp1,iadp3,23,8)
c         r2=r1-CSHIFT(r1,2,1)
           call CM_f_news_sub_always_3_1L(iadr2,iadr1,iadr1,0,
     1          CM_upward,23,8)
c         r2=r2-CSHIFT(r2,2,-1)
           call cm_f_news_sub_always_2_1L(iadr2,iadr2,0,
     1          cm_downward,23,8)
c         r2=r1-alp*r2
           call cm_f_mult_subf_1L(iadr2,iadalp,iadr2,iadr1,23,8)
c         where(xedges) r2=0.0
       call CM_load_context(iadxedges)
       call CM_f_move_zero_1L(iadr2,23,8)
       call CM_set_context
c
c####    call factapp(r2,ntaupow2,nxpow2,factor1,factor2,nxlog2,loopcnt)
               lbound=1
               do 23056 i=1,loopcnt 
c         p3=r2
       call CM_f_move_always_1L(iadp3,iadr2,23,8)
c          r1=cshift(p3,dim=2,shift=-lbound)
       call getfrompowertwo(iadr1,iadp3,0,i-1,
     1      CM_downward,32)
          call CM_f_mult_add_1L(iadr2,factor1(i),iadr1,iadr2,23,8)
c           r1=cshift(p3,dim=2,shift=lbound)
       call getfrompowertwo(iadr1,iadp3,0,i-1,
     1      CM_upward,32)
          call CM_f_mult_add_1L(iadr2,factor2(i),iadr1,iadr2,23,8)
                  lbound=lbound*2
c%       call CM_f_ne_zero_1L(iadr2,23,8)
c%       call CM_store_test(iadnonzero)
c%       call CM_latch_leds_always(iadnonzero)
23056             continue
c          p4=r2-p2
       call CM_f_mult_sub_1L(iadp4,iadr2,factor1(loopcnt+1),
     1      iadp2,23,8)
c           if(it>325 && it<365) {
c%         call CM_f_move_always_1L(iada,iadp4,23,8)
c	       call putcomp(nxpow2,ntaupow2,iada,pfe)
c	       }
23054       continue
            if(.not.(stepdown .eq. 1))goto 23058
               if(.not.(it .eq. 2))goto 23060
       call CM_load_context(iadoslot)
       call CM_f_move_1L(iadpsec,iadp4,23,8)
23060          continue
               goto 23059
c           else
23058          continue
               itau=(it-2)/(stepdown-1)
               irem=it-2-itau*(stepdown-1)
               jt=it+itau-1
               if(.not.(itau.gt.0 .and. itau .le. ntau .and. jt .le. nt)
     &          )goto 23062
c      print *,' it=',it,' itau=',itau,' irem=',irem,' jt=',jt
c		    oslot=cshift(oslot,1,1)
       call CM_get_from_news_always_1L(iadoslot,iadoslot,1,
     1      CM_upward,1)
                  if(.not.(mod(jt,ntaupow2).eq.0))goto 23064
      call CM_set_context
      call cm_f_max_constant_3_1L(iadr1,iadpout,-clip,23,8)
      call cm_f_min_constant_2_1L(iadr1,clip,23,8)
      call cm_f_multiply_const_always_2_1L(iadr1,ovc,23,8)
      call cm_s_f_floor_2_2L(iadplot,iadr1,8,23,8)
      call CMFB_write_always(dispid,CMFB_green,iadplot,
     1     0,jt)
       iadpout=CM_add_offset_to_field_id(iadpout,-32)
23064             continue
       call CM_load_context(iadoslot)
       call CM_f_move_1L(iadpout,iadp4,23,8)
23062          continue
               if(.not.(irem.eq.0))goto 23066
                  itau=itau-1
                  jt=jt-1
                  if(.not.(itau.gt.0 .and. itau .le. ntau .and. jt .le. 
     &             nt))goto 23068
c      print *,' it=',it,' itau=',itau,' irem=',irem,' jt=',jt
c		    oslot=cshift(oslot,1,1)
       call CM_get_from_news_always_1L(iadoslot,iadoslot,1,
     1      CM_upward,1)
                     if(.not.(mod(jt,ntaupow2).eq.0))goto 23070
      call CM_set_context
      call cm_f_max_constant_3_1L(iadr1,iadpout,-clip,23,8)
      call cm_f_min_constant_2_1L(iadr1,clip,23,8)
      call cm_f_multiply_const_always_2_1L(iadr1,ovc,23,8)
      call cm_s_f_floor_2_2L(iadplot,iadr1,8,23,8)
      call CMFB_write_always(dispid,CMFB_green,iadplot,
     1     0,jt)
       iadpout=CM_add_offset_to_field_id(iadpout,-32)
23070                continue
       call CM_load_context(iadoslot)
       call CM_f_move_1L(iadpout,iadp4,23,8)
23068             continue
23066          continue
23059       continue
       call CM_set_context
c          p2=cshift(p1,1,1)
       call CM_get_from_news_always_1L(iadp2,iadp1,1,
     1      CM_upward,32)
c	   p1=p4
       call CM_f_move_always_1L(iadp1,iadp4,23,8)
c	   p3=cshift(p1,1,1)
       call CM_get_from_news_always_1L(iadp3,iadp1,1,
     1      CM_upward,32)
c
23050       continue
c all done migrating; write output section
      call cm_stop_timer(1)
         print *,'writing migrated section'
         iadpout=iadpsec
         call putsec(outfd,nt,nx,npad,ntaupow2,iadpout,pfe,.true.)
c deallocate tridiagonal factors
         do 23072 i=0,loopcnt 
       call cm_deallocate_heap_field(factor1(i+1))
23072       continue
c put top of migrated image on display
      call CM_set_context
      call cm_f_max_constant_3_1L(iadr1,iadpsec,-clip,23,8)
      call cm_f_min_constant_2_1L(iadr1,clip,23,8)
      call cm_f_multiply_const_always_2_1L(iadr1,ovc,23,8)
      call cm_s_f_floor_2_2L(iadplot,iadr1,8,23,8)
      call CMFB_write_always(dispid,CMFB_green,iadplot,
     1     0,0)
         goto 23032
c     endwhile
23033 continue
c
      print *,'done'
      return
      end
      subroutine getfrompowertwo ( d , s , a , l2d , dir , len )
      implicit none
      integer fetch, hetch, tetch, getch, auxpar,auxputch
      integer d,s,a,l2d,dir,len
      integer n,i
      if(.not.(l2d.lt.1))goto 23074
      call CM_get_from_news_always_1L(d,s,a,dir,len)
c     if(l2d.ge.0) return
         n=2**l2d
         do 23076 i=2,n 
      call CM_get_from_news_always_1L(d,d,a,dir,len)
23076       continue
         goto 23075
c     else
23074    continue
      call CM_get_from_power_two_always_1L(d,s,a,l2d,dir,len)
23075 continue
      return
      end
