PROGRAM NMRTAB
      IMPLICIT none
      integer*4 iargc,nat,iz,iop
      integer*4,allocatable:: qz(:)
      real*8,allocatable::nmr(:)
      real*8 sref
      character*80 s
      if(iargc().lt.2)then
       write(6,60)
60     format(' Usage: nmrtab iz  s_ref (iop)',/,
     1 ' iop   : 0(nothing) ... take all shifts as they are',/,
     2 '         1          ... delete OH/NH',/,
     3 '        10          ... sum up CH3 and CH2')
       stop
      endif
      call getarg(1,s)
      read(s,*)iz
      call getarg(2,s)
      read(s,*)sref
      if(iargc().gt.2)then
       call getarg(3,s)
       read(s,*)iop
      else
       iop=0
      endif
      call rdnmr(nat,nmr,qz,0)
      allocate(nmr(nat),qz(nat))
      call rdnmr(nat,nmr,qz,1)
      call fixh(iop,nat,nmr,qz,iz)
      call wnmrtab(iz,sref,nat,nmr,qz)
      end
c     ============================================================
      subroutine wnmrtab(ir,sref,nat,nmr,iz)
      implicit none
      integer*4 nat,i,iz(*),ir
      real*8 nmr(*),sref
      open(9,file='NMR.TAB')
      write(9,900)
900   format('Seleted NMR shifts',/)
      write(9,901)
901   format(60(1h-))
      do 1 i=1,nat
1     if(iz(i).eq.ir)write(9,902)i,sref-nmr(i),ir
902   format(i5,f12.4,i5)
      write(9,901)
      close(9)
      write(6,*)'NMR.TAB writen'
      return
      end
c     ============================================================
      subroutine rdnmr(nat,nmr,iz,ic)
      implicit none
      integer*4 nat,i,iz(*),ic
      real*8 nmr(*)
      open(9,file='FILE.NMR')
      read(9,*)nat
      if(ic.eq.1)then
       read(9,*)
       read(9,*)
       do 1 i=1,nat
1      read(9,*)iz(i),iz(i),nmr(i)
       read(9,*)
      endif
      close(9)
      write(6,*)'FILE.NMR read'
      return
      end
c     ============================================================
      subroutine fixh(iop,nat,nmr,qz,iz)
      implicit none
      integer*4 iop,nat,i,j,iz,ix,nhmax,nh,qz(*)
      parameter (nhmax=4)
      integer*4 hlist(nhmax)
      real*8 nmr(*),xi,yi,zi,xj,yj,zj,d,ave
      real*8,allocatable::r(:)
      if(iop.eq.0)return
      if(iz.ne.1)call report('iop possible only for iz = 1')
      open(9,file='FILE.X',status='old')
      read(9,*)
      read(9,*)
      allocate(r(3*nat))
      do 1 i=1,nat
1     read(9,*)r(1+3*(i-1)),(r(ix+3*(i-1)),ix=1,3)
      close(9)

c     delete OH,NH:      
      if(mod(iop,10).eq.1)then
       write(6,*)'deleting OH ON'
       do 2 i=1,nat
       if(qz(i).eq.1)then
        xi=r(1+3*(i-1))
        yi=r(2+3*(i-1))
        zi=r(3+3*(i-1))
        do 3 j=1,nat
        if(qz(j).eq.7.or.qz(j).eq.8)then
         xj=r(1+3*(j-1))
         yj=r(2+3*(j-1))
         zj=r(3+3*(j-1))
         d=dsqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2)
         if(d.lt.1.2d0)then
          nmr(i)=9999.9d0
          goto 2
         endif
        endif
3       continue
       endif
2      continue
      endif

c     sum up CH3,CH2,NH2,etc.
      if(mod(iop,100).gt.10)then
       write(6,*)'summing up XH_N'
       do 4 i=1,nat
       if(qz(i).ne.1)then
        xi=r(1+3*(i-1))
        yi=r(2+3*(i-1))
        zi=r(3+3*(i-1))
        nh=0
        do 5 j=1,nat
        if(qz(j).eq.1)then
         xj=r(1+3*(j-1))
         yj=r(2+3*(j-1))
         zj=r(3+3*(j-1))
         d=dsqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2)
         if(d.lt.1.2d0)then
          nh=nh+1
          if(nh.gt.nhmax)call report('nh > nhmax')
          hlist(nh)=j
         endif
        endif
5       continue
        if(nh.gt.2)then
         ave=0.0d0
         do 6 j=1,nh
6        ave=ave+nmr(hlist(j))
         ave=ave/dble(nh)
         do 7 j=1,nh
7        nmr(hlist(j))=ave
        endif
       endif
4      continue
      endif

      return
      end
c     ============================================================
      subroutine report(o)
      character*(*) o
      write(6,*)o
      stop
      end
c     ============================================================