diff --git a/tecs/tecs_plot.f b/tecs/tecs_plot.f
new file mode 100644
index 00000000..f9386a8f
--- /dev/null
+++ b/tecs/tecs_plot.f
@@ -0,0 +1,711 @@
+subroutine tecs_plot(auxpar)
+
+  character(len=*) auxpar
+
+  integer dmax, nmax, tmax, amax, nmenu, chartperiod, naux
+  parameter (dmax=400, nmax=9, tmax=8, amax=3, nmenu=13, chartperiod=5, naux=1)
+  integer minRange, maxRange, oneDay
+  parameter (minRange=60, maxRange=7*24*3600, oneDay=24*3600)
+  integer zoom, right, live
+  parameter (zoom=1, right=2, live=3)
+  real winmin, undef
+  parameter (winmin=60., undef=-1.125/1024./1024./1024.)
+
+  real*4 x1,x2,xmin,xmax,ymin(2),ymax(2),window
+  real*4 xd(dmax, nmax),yd(dmax,nmax), yy0(nmax), yy1(nmax)
+  real*4 ylast1,ylast2,y1,y2
+  real*4 ex,ey,fx,fy,row,ticks,tim0,tim1,menuwid
+  integer l,j,i,n,t,leng,i1,i2,rl,startday,thisday
+  integer ncol, nset, mode
+  integer first,last,step,tbase,lastj
+  integer colorList(nmax)/5,3,2,4,6,8,14,15,8/
+  integer color(nmax)
+  integer retLen(nmax)
+  integer sel/0/, sel1, sel2, auxsel/1/
+  character key*1
+  character text(2,nmenu)*16/ &
+     'L'    ,'live off'  &
+    ,'z'    ,'sel. zoom' &
+    ,'+'    ,'zoom in'   &
+    ,'-'    ,'zoom out'  &
+    ,'x'    ,'show all'  &
+    ,'2d'   ,'2 days'   &
+    ,'1h'   ,'1 hour'   &
+    ,'15m'  ,'15 min'  &
+    ,'31.7t','goto date' &
+    ,'c'    ,'show P/He/Aux' &
+    ,'s'    ,'select T' &
+    ,'f'    ,'write file' &
+    ,'q'    ,'quit'/
+  character weekdays(7)*4/'Mon','Tue','Wed','Thu','Fri','Sat','Sun'/
+  character buf*8, device*8, name*40, filnam*128, numb*16, title*64, pars*64
+  character(len=4) tpar(tmax)/'Te', 'Tr', 'Tm', 'Ts', 'T1', 'T2', 'T3', 'T4'/
+  character(len=4) apar(amax)/'P', 'He', 'Aux'/
+  character(len=16) parnam(nmax)
+  logical gap, done
+  logical saveit
+  integer iret, lund, numl, mon, day
+  integer pars_len, title_len, text_len, name_len
+
+! functions
+  integer sys_gmt_off, myc_now, myc_time, myc_date, get_data, tecs_get_mult, tecs_get_par
+
+  data window/0./
+
+  if (window==0) window=1800.
+  saveit=.false.
+  mode=live
+  call pgopen(' ')
+
+  call pgqinf('TYPE', device, l)
+  if (device=='NULL') then
+    call pgclos
+    call pgopen('?')
+    call pgqinf('TYPE', device, l)
+    if (device=='NULL') then
+      print *,'No PGPLOT-Device defined'
+      goto 9
+    endif
+  endif
+  
+  call pgqcol(i,ncol)
+!  print *,i,ncol,' colors ',device
+  if (ncol>=8 .and. device /= 'VT125') then
+    call pgscr(0, 1.0, 1.0, 1.0)
+    call pgscr(1, 0.0, 0.0, 0.0)
+    call pgscr(2, 1.0, 0.0, 0.0)
+    call pgscr(3, 0.0, 1.0, 0.0)
+    call pgscr(4, 0.0, 0.0, 1.0)
+    call pgscr(5, 0.0, 1.0, 1.0)
+    call pgscr(6, 1.0, 0.0, 1.0)
+    call pgscr(7, 1.0, 1.0, 0.0)
+    call pgscr(8, 1.0, 0.5, 0.0)
+  endif
+
+  call pgask(.false.)
+  call pgupdt
+  l=0
+  x1=0
+  step=0
+  do i=1,amax
+    if (auxpar == apar(i)) then
+      auxsel=i
+    endif
+  enddo
+
+1 pars=' '
+  nset=0
+  do i=1,tmax
+    nset=nset+1
+    color(nset)=colorList(i)
+    parnam(nset)=tpar(i)
+    call str_trim(pars, pars, pars_len)
+    pars=pars(1:pars_len)//' '//parnam(nset)
+  enddo
+  nset=nset+1
+  color(nset)=colorList(3)
+  parnam(nset)=apar(auxsel)
+  call str_trim(pars, pars, pars_len)
+  pars=pars(1:pars_len)//' '//parnam(nset)
+
+  iret=tecs_get_mult(pars, t, nset, yy0)
+  if (iret < 0) goto 99
+  if (window == 0) then
+    last=t
+    step=maxRange/dmax
+    window=maxRange
+    first=t-min(dmax*step-1,maxRange-step)
+  else if (mode >= right) then
+    step=window/(dmax-2)+0.99
+    last=t
+    first=t-min(dmax*step-1,nint(window))
+  else
+    if (mode==zoom) then
+      x2=(x1+x2+window)/2
+      x1=x2-window
+    endif
+    if (x1 .gt. x2-minRange) x1=x2-minRange
+    step=(x2-x1)/(dmax-2)+0.99
+    last=x2+tbase
+    first=x1-step+tbase
+  endif 
+  if (step == 0) step=1
+
+  if (step>60) then   ! normalize step
+    step=(step+59)/60*60
+  else if (step>30) then
+    step=60
+  elseif (step>20) then
+    step=30
+  else
+    step=(step+4)/5*5
+  endif
+  first=last-(last-first+step-1)/step*step ! round first
+!  print *,step,last-first
+  tbase=first-mod(first,7*oneDay)
+  iret=get_data(pars, first, last, step, tbase, xd, yd, dmax, nmax, retLen)
+  if (iret < 0) goto 99
+
+  x2 = last - tbase
+  if (mode >= right) then
+    x1=x2-window
+  else
+    x1 = first - tbase
+  endif
+  tim0=t-tbase
+  do i=1,nset
+    leng=retLen(i)
+    do while (leng > 1 .and. yd(leng,i) == undef)
+      leng=leng-1
+    enddo
+    if (leng == 1) leng=0
+    retLen(i)=leng
+    if (mode==live .and. leng>0) then
+      xd(leng,i)=tim0
+      yd(leng,i)=yy0(i)
+    endif
+  enddo
+  if (sel /= 0) then
+    do while (sel < nset .and. retLen(sel) == 0)
+      sel=sel+1
+    enddo
+    if (sel >= nset) sel = 0
+  endif
+  if (saveit) goto 9
+  if (mode==live) then
+    x2=max(tim0,x2)+min(1800., window*0.5)
+  endif
+
+  if (window>50*3600) then
+    ticks=8*3600
+  elseif (window>25*3600) then
+    ticks=4*3600
+  else
+    ticks=0.0 ! automatic
+  endif
+
+  i1=1
+  i2=nset-naux
+  if (sel==0) then
+    sel1=i1
+    sel2=i2
+  else
+    sel1=sel
+    sel2=sel
+  endif
+  do rl=1,2
+    if (mode==zoom .and. rl==1) then
+      ymin(1)=y1
+      ymax(1)=y2
+    else
+      ymin(rl)=1e30
+      ymax(rl)=-1e30
+      ylast1=ymin(rl)
+      ylast2=ymax(rl)
+      do i=sel1,sel2
+        do j=1,retLen(i)
+          if (yd(j,i)/=undef) then
+            ymin(rl)=min(ymin(rl),yd(j,i))
+            ymax(rl)=max(ymax(rl),yd(j,i))
+          endif
+        enddo
+        do j=max(1,retLen(i)-4),retLen(i)
+          if (yd(j, i)/=undef) then
+            ylast1=min(ylast1,yd(j, i))
+            ylast2=max(ylast2,yd(j, i))
+          endif
+        enddo
+      enddo
+
+      ey=(ymax(rl)-ymin(rl))
+      fy=abs(ymax(rl))
+      ymax(rl)=ymax(rl)+max(fy*0.0075,ey*0.01)
+      ymin(rl)=ymin(rl)-max(fy*0.005,ey*0.01)
+      if (mode==live) then
+        ymin(rl)=min(ymin(rl),max(0.0,ylast1-ey*0.4))
+        ymax(rl)=max(ymax(rl),ylast2+ey*0.4)
+      endif
+    endif
+
+    if (ymax(rl) < ymin(rl)) then
+      ymax(rl)=1.0
+      ymin(rl)=0
+    elseif (ymax(rl) == ymin(rl)) then
+      ymax(rl)=ymin(rl)*1.00001+1.0
+      ymin(rl)=-1.0e-3
+    endif
+
+    call set_win(rl,x1,x2,ymin(rl),ymax(rl))
+
+    do i=i1,i2
+      call pgsci(color(i))
+      l=0
+      lastj=1
+      do j=1,retLen(i)
+        if (yd(j,i)==undef) then
+          if (j>lastj) then
+            call pgline(j-lastj, xd(lastj,i), yd(lastj,i))
+          endif
+          lastj=j+1
+        else
+          l=j
+        endif
+      enddo
+      if (retLen(i) > lastj) call pgline(retLen(i)+1-lastj, xd(lastj,i), yd(lastj,i))
+      retLen(i)=l
+    enddo
+    call pgsci(1)
+    if (rl == 1) then
+      call pgsch(1.0)
+      call pgtbox('ZHXYBINST', ticks, 0, 'BCINMST', 0.0, 0)
+      call pgtbox('C', 0.0, 0, ' ', 0.0, 0)
+      ey=0.0
+      do i=i1,i2
+        if (retLen(i) > 0) then
+          name=parnam(i)
+          if (name=='Tm') then
+            name='Main Sensor'
+          elseif (name=='Ts') then
+            name='Sample Sensor'
+          elseif (name=='Tr') then
+            name='SetPoint'
+          endif
+          call str_trim(name, name, name_len)
+          if (sel == i) then
+            name=name(1:name_len)//'*'
+            call str_trim(name, name, name_len)
+          endif
+          call pglen(5, name(1:name_len), fx, fy)
+          call pgsci(color(i))
+          call pgmtxt ('L', 2.5, ey, 0.0, name(1:name_len))
+          ey=ey+fy+0.04
+        endif
+      enddo
+      call pgsci(1)
+      call pgmtxt ('L', 2.5, ey, 0.0, 'T [K]')
+    else
+      call pgsch(0.7)
+      call pgtbox('ZCIST', ticks, 0, 'BCVINMST', 0.0, 0)
+      call pgtbox('B', 0.0, 0, ' ', 0.0, 0)
+      call pgsci(color(nset))
+      call pgsch(1.0)
+      if (parnam(nset) == 'P' .or. parnam(nset) == 'p') then
+        call pgmtxt ('L', 2.5, 0.5, 0.5, 'Power [W]')
+      elseif (parnam(nset) == 'He') then
+        title='%'
+        iret=tecs_get_par('heUnits', title, 0)
+        call str_trim(title, title, title_len)
+        call pgmtxt ('L', 2.5, 0.5, 0.5, 'Helium ['//title(1:title_len)//']')
+      else
+        call pgmtxt ('L', 2.5, 0.5, 0.5, parnam(nset))
+      endif
+    endif
+    i1=nset-naux+1
+    i2=nset
+    sel1=i1
+    sel2=i2
+  enddo
+
+  call pgsch(0.7)
+  rl=1
+  call set_win(rl,x1,x2,ymin(rl),ymax(rl))
+
+  call pgsci(1)
+  call pgsclp(0)
+  if (mode==live) then
+    text(2,1)='live off'
+  else
+    text(2,1)='live on'    
+  endif
+  text(2,10)='show '//apar(mod(auxsel,3)+1)
+  menuwid=0.0
+  do i=1,nmenu
+    call str_trim(text(2,i), text(2,i), text_len)
+    call pglen(5, text(2,i)(1:text_len), fx, fy)
+    call pgmtxt('T', 3.0, menuwid, 0.0, '|'//text(1,i))
+    call pgmtxt('T', 2.5, menuwid, 0.0, '|')
+    call pgmtxt('T', 2.0, menuwid, 0.0, '|'//text(2,i))
+    menuwid=menuwid+fx+0.01
+  enddo
+  call pgmtxt('T', 3.0, menuwid, 0.0, '|')
+  call pgmtxt('T', 2.5, menuwid, 0.0, '|')
+  call pgmtxt('T', 2.0, menuwid, 0.0, '|')
+  menuwid=menuwid+0.01
+
+  call pgsch(0.8)
+  done=.false.
+  row=(ymax(rl)-ymin(rl))/30.
+  ey=ymin(rl)-row*3.5
+  i=max(0,int((x1+oneDay/2)/oneDay))
+
+  do
+    ex=(i+0.5)*oneDay
+    if (ex > x2) EXIT
+    done=.true.
+    write(buf,'(i8.8)') myc_date(nint(ex)+tbase)
+    call pgptxt(ex, ey, 0.0, 0.5, weekdays(mod(i,7)+1)//buf(7:8)//'.'//buf(5:6))
+    ex=ex-12*3600
+    if (ex > x1) then
+      call pgmove(ex, ey)
+      call pgdraw(ex, ey+row)
+    endif
+    ex=ex+oneDay
+    if (ex < x2) then
+      call pgmove(ex, ey)
+      call pgdraw(ex, ey+row)
+    endif
+    i=i+1
+  enddo
+  if (.not. done) then
+    n=nint(x2)/oneDay*oneDay
+    i=nint(x1)-n
+    j=nint(x2)-n
+    if (i < 0) then
+      if (-i > j) then
+        ex=0.0
+        i=nint(x1)+tbase
+      else
+        ex=1.0
+        i=nint(x2)+tbase
+      endif
+    else
+      ex=0.5
+      i=nint(x2)+tbase
+    endif
+    thisday=mod(i/oneDay,7)+1
+    write(buf,'(i8.8)') myc_date(i)
+!    call pgptxt((x1+x2)/2, ey, 0.0, ex, weekdays(mod(thisday,7)+1)//buf(7:8)//'.'//buf(5:6))
+    call pgmtxt('B', 3.5, ex, ex, weekdays(thisday)//buf(7:8)//'.'//buf(5:6))
+  endif
+
+  iret=tecs_get_par('device', title, 0)
+  if (iret < 0) goto 99
+  i=index(title, '(')
+  if (i > 2 ) then
+    title=title(1:i-1)
+  else
+    title='test - no device'
+  endif
+
+  call pgmtxt('T', -1.5, 0.02, 0.0, title)
+  call pgsclp(1)
+
+  call purge_keys ! purge buffer
+
+  numl=0
+  numb=' '
+7 ex=undef
+  if (mode==live) then
+!    if (device(1:1)=='X') then
+!      call pgmtxt('T', 0.5, 0.0, 0.0, 'LIVE MODE (click on text window before pressing any further key)')
+!    endif
+!    call get_key(key, 0, chartperiod)
+    i=chartperiod-mod(myc_now(), chartperiod)
+    call get_cursor(ex, ey, key, -i)
+    
+    do while (key == char(0) .or. key == ' ') ! no key or space key pressed
+      iret=tecs_get_mult(pars, t, nset, yy1)
+      if (iret<0) goto 99
+      tim1=t-tbase
+      if (tim1 > x2) then
+        call pgpage
+	window=x2-x1
+        goto 1
+      endif
+      if (tim1 > tim0) then
+        i1=1
+        i2=nset-naux
+        do rl=1,2
+          call set_win(rl,x1,x2,ymin(rl),ymax(rl))
+          do i=i1,i2
+            if (yy0(i) /= undef .and. yy1(i) /= undef) then
+              if ((sel==0 .or. sel==i) .and. (yy1(i) < ymin(rl) .or. yy1(i) > ymax(rl))) then
+                call pgpage
+	        window=x2-x1
+                goto 1
+              endif
+              call pgsci(color(i))
+              call pgmove(tim0, yy0(i))
+	      call pgdraw(tim1, yy1(i))
+            endif
+            yy0(i)=yy1(i)
+          enddo
+          i1=nset-naux+1
+          i2=nset
+        enddo
+        tim0=tim1
+
+      endif
+      i=chartperiod-mod(myc_now(), chartperiod)
+      call get_cursor(ex, ey, key, -i)
+    enddo
+  else
+    call get_cursor(ex, ey, key, 0)
+  endif
+  rl=1
+  call set_win(rl,x1,x2,ymin(rl),ymax(rl))
+
+8 if (key>='a') key=char(ichar(key)-32)
+  if (key=='-') then
+    window=min(window*2, 8.0*oneDay)
+    if (mode==zoom) then
+      x1=x1-(x2-x1)/2
+      x2=x2+(x2-x1)/3
+      y1=y1-(y2-y1)/2
+      y2=y2+(y2-y1)/3
+    endif
+  elseif (key=='X') then
+    window=0
+    mode=0
+  elseif (key=='+' .or. key==',') then
+    window=max(winmin,window/2)
+    if (mode==zoom) then
+      if (ex==undef) then
+        ex=(x1+x2)/2
+        ey=(min(y2,ymax(1))+max(y1,ymin(1)))/2
+      end if
+      fx=max(winmin,x2-x1)
+      fy=max(y2-y1,1e-3,y2*1e-5)
+      x1=ex-fx/4
+      x2=ex+fx/4
+      y1=ey-fy/4
+      y2=ey+fy/4
+    endif
+  elseif (key=='Z') then
+    call pgsci(1)
+    if (ex==undef) then
+      call pgmtxt('T', 0.5, 0.0, 0.0, 'click on two opposite corners of a selection rectangle')
+      call get_cursor(ex, ey, key, 0)
+    else
+      call pgmtxt('T', 0.5, 0.0, 0.0, 'click on second corner of selection rectangle')
+    endif
+    call pgsci(6)
+    xmin=x1
+    xmax=x2
+      call pgmove(xmin, ey)
+      call pgdraw(xmax, ey)
+      call pgmove(ex, ymin(rl))
+      call pgdraw(ex, ymax(rl))
+    if (device(1:1) .eq. 'X') then
+      fx=ex
+      fy=ey
+      call get_cursor(fx, fy, key, 2)
+    else
+      call get_cursor(fx, fy, key, 0)
+    endif
+    x1=max(xmin,min(ex,fx))
+    x2=min(xmax,max(ex,fx))
+    if (x1>=x2) then
+      x1=xmin
+      x2=xmax
+    endif
+    window=x2-x1
+    y1=max(ymin(1),min(ey,fy))
+    y2=min(ymax(1),max(ey,fy))
+    if (y1>=y2) then
+      y1=ymin(1)
+      y2=ymax(1)
+    endif
+    mode=zoom
+  elseif (key >= '0' .and. key <= '9' .or. key == '.') then ! number
+    if (numl .lt. len(numb)) then
+      numl=numl+1
+      numb(numl:numl)=key
+      call pgsci(1)
+      call pgsch(0.8)
+      call pgmtxt('T', 2.0, menuwid, 0.0, numb(1:numl))
+    endif
+!    call get_cursor(ex, ey, key, -chartperiod)
+!    if (key/=char(0)) goto 8
+    if (mode==zoom) mode=0
+    goto 7
+  elseif (key == 'D') then
+    ex=1
+    read(numb, *, iostat=i) ex
+    window=min(maxRange,max(minRange, nint(oneDay*ex)))
+    if (mode < right) mode=right
+    x1=0
+  elseif (key == 'H') then
+    ex=1
+    read(numb, *, iostat=i) ex
+    window=min(maxRange,max(minRange, nint(3600*ex)))
+    if (mode < right) mode=right
+    x1=0
+  elseif (key == 'M') then
+    ex=1
+    read(numb, *, iostat=i) ex
+    window=min(maxRange,max(minRange, nint(60*ex)))
+    if (mode < right) mode=right
+    x1=0
+  elseif (key == 'T' .or. numl>0 .and. (key==char(13) .or. key==char(10))) then
+    j=index(numb,'.')
+    day=0
+    if (j > 1 .and. j < numl) then
+      read(numb(1:j-1), *, iostat=i) day
+      mon=0
+      read(numb(j+1:numl), *, iostat=i) mon
+      tbase = myc_time(day+mon*100)
+    else
+      read(numb, *, iostat=i) day
+      tbase = myc_time(day)
+    endif
+    x1=0
+    x2=oneDay
+    window=x2
+    mode=0
+  elseif (key == 'L') then
+    if (mode == live) then
+      mode=right
+    else
+      mode=live
+    endif
+  elseif (key == 'F') then
+    saveit=.true.
+  elseif (key=='Q' .or. key==char(13) .or. key==char(10)) then
+    goto 9
+  elseif (key == 'S') then
+    sel=sel+1
+    if (sel > tmax) sel=0
+    if (mode==zoom) mode=0
+  elseif (key == 'C') then
+    auxsel=auxsel+1
+    if (auxsel > amax) auxsel=1
+  elseif (mode==live) then
+    goto 7
+  endif
+  numl=0
+  numb=' '
+  call pgpage
+  goto 1
+99 call tecs_write_msg(6)
+9 continue
+  call pgclos
+  call purge_keys
+  print *
+  if (saveit) then
+    lund=41
+    print '(x,a,$)', 'Filename: '
+    read(*,'(a)') filnam
+    open(lund, file=filnam, status='unknown') ! , carriagecontrol='list')
+
+    l=0
+    i2=0
+    do i1=1,nset
+      if (i1 > nset-naux) then
+        j=i1
+      else
+        j=nset-i1+(1-naux)
+      endif
+      
+      gap=.false.
+      do i=1,retLen(j)
+        if (yd(i,j)==undef) then
+          if (gap) then
+            write(lund, '(f9.4,2a)') xd(i,j)/3600., char(9), ' '
+            gap=.false.
+          endif
+        else
+          if (i2/=i1) then
+            if (l > 0) write(lund, *)
+            call str_trim(parnam(j), parnam(j), text_len)
+            write(lund, *) ' time [h]',char(9), '   ',parnam(j)(1:text_len)
+            l=l+1
+            i2=i1
+          endif
+          write(lund, '(f9.4,a,f9.4)') xd(i,j)/3600., char(9), max(-999.,min(9999.,yd(i,j)))
+          l=l+1
+          gap=.true.
+        endif
+      enddo
+    enddo
+    close(lund)
+    print *, l, ' lines written to ',filnam(1:48)
+  endif
+
+end subroutine
+
+subroutine get_cursor(x, y, key, mode)
+  character*1 key
+  real*4 x, y
+  integer mode
+  integer with_timeout/-1/
+  integer l
+  character res*32
+
+  if (with_timeout<0) then
+    with_timeout=0
+    call pgqinf('VERSION', res, l)
+    if (res(l:l)=='+') then
+      call pgqinf('TYPE', res, l)
+      if (res(1:1)=='X') then
+        with_timeout=1
+      end if
+    end if
+  end if
+  if (with_timeout>0 .or. mode>=0) then
+    call pgband(mode, 0, x, y, x, y, key)
+  else
+    call sys_get_key(key, -mode)
+  endif
+end subroutine
+
+subroutine purge_keys
+  character key*1
+  key=' '
+  do while (key/=char(0))
+    call sys_get_key(key, 0)
+  end do
+end subroutine
+
+subroutine set_win(rl, x1, x2, y1, y2)
+  integer rl
+  real x1, x2, y1, y2
+
+  if (rl == 1) then
+    call pgsvp(0.07,0.93,0.3,0.9)
+  else
+    call pgsvp(0.07,0.93,0.01,0.20)
+  endif
+
+  call pgswin(x1,x2,y1,y2)
+end subroutine
+
+integer function get_data(pars, first, last, step, tbase, xd, yd, dmax, nmax, retlen)
+  character pars*(*)
+  integer first, last, step, tbase, dmax, nmax, retlen(nmax)
+  real*4 xd(dmax,nmax), yd(dmax,nmax)
+
+  integer oneDay, maxn
+  parameter (oneDay = 24*3600, maxn=9)
+  integer tecs_get_data
+
+  integer i,j,rl(maxn),m,k,n,mm
+
+  if (nmax > maxn) stop 'get_data: nmax>maxn'
+  if (last-first <= oneDay) then
+    get_data=tecs_get_data(pars, first, last, step, tbase, xd, yd, dmax, nmax, retLen)
+  else
+    do j=1,nmax
+      retlen(j)=0
+    enddo
+    m=0
+    do i=first/oneDay,last/oneDay
+      get_data=tecs_get_data(pars, max(first,i*oneDay), min(last,(i+1)*oneDay-step), step, tbase &
+                        , xd(m+1, 1), yd(m+1, 1), dmax, nmax, rl)
+      if (get_data<0) return
+      mm=0
+      do j=1,nmax
+        n=retlen(j)
+        do k=m+1,m+rl(j)
+          n=n+1
+	  xd(n,j)=xd(k,j)
+	  yd(n,j)=yd(k,j)
+        enddo
+        retlen(j)=n
+        mm=max(mm,n)
+      enddo
+!      print *,mm-m,' points read'
+      m=mm
+    enddo
+  endif
+end function