    	subroutine  rm_err_pick
	include 'common_arriv.inc'
	k11=neq_branch+1
        keq=meqa(k11)
	nobsp_org=nobsp
	nobss_org=nobss
cc	write(*,*) 'after dkread_pick_parm',iy,im,id,ih,imin,t,phi,ram,dep
	nobsallk=nobsall
		do 10 kk=1,nobsall
	     kst=kstn_plot_ord(kk)
	     if(kst.le.0) go to 10
	     nst=nstn(kst)
c       write(*,110) stnw(kst),nst,deltak(kst),po(kst),posakai(kst),
c     1   wp(kst),so(kst),sosakai(kst),ws(kst),resoc(kst,1),
c     2 res_sakai(kst,1),resoc(kst,2),res_sakai(kst,2),nc_removea(kst,1),
c     2 nc_removea(kst,2),iflag_repa(kst,1),iflag_repa(kst,2)
c      write(74,110) stnw(kst),ns,deltak(kst),po(kst),wp(kst),so(kst),
c     1 ws(kst),posakai(kst),res_sakai(kst,1),sosakai(kst),
c     2 res_sakai(kst,2),nc_removea(kst,1),nc_removea(kst,2),
c     3 iflag_repa(kst,1),iflag_repa(kst,2)
 110	   format(a7,i4,f6.1,' Ps,Theo,w=',6f7.2,' res=',4f8.3,
     3 ' ncer=',2i3,' flag=',2i3)

	     
	   do 20 kps=1,2
	      err_min2=err_min
	  if(kps.eq.2)err_min2= err_min2*1.5 
	      w=wp(kst)
	      if(kps.eq.2) w=ws(kst)
	      if(w.le.0.0001) go to 20
	      nst=nstn(kst)
	resvar=res_sakai(kst,kps)     !Varriance of theoretical arrival
	if(resvar.lt.err_min2) resvar=err_min2
	
	resmax=resvar*rat_stn_sd2
	res=abs(resoc(kst,kps))
	if(res.gt.resmax) then
c	   write(*,123) stnw(kst),kps,res_sakai(kst,kps),resoc(kst,kps),
c     1  resmax
 123	 format('removed stn ',a7,' kps=',i2,3f7.2)
	   if(kps.eq.1) wp(kst)=0.0
	   if(kps.eq.2) ws(kst)=0.0
	   wei_crossa(k11,nst,kps)=0.0
	   endif
 20	   continue
 10	   continue
cc	   read(*,*) iok
	   return
	  end



	subroutine  det_arriv_error
ccc     po(k)=po(2,nst)-po(1,nst)=999.9 if one of P are not measured.
ccc     travp(k)=trvp_cor(2,nst)-trvp_cor(1,nst)+t  Defined for all stn.
ccc     wp(k)=-abs(wp(k)) if large residual stn
ccc     resoc(k,1)={po(2,nst)-po(1,nst)}-travp(k)
ccc     if wp(1,nst)>0, wp(2,nst)<0, ptheo=po(1,nst)+travp(k)
ccc     res_dd_crossa(mmm,nst,1)=resoc(kst,1)
ccc     travp_dd_crossa(mmm,nst)=trvp_cor(2,nst)-trvp_cor(1,nst)+t
ccc	wei_dd_crossa(mmm,nst)=wp(kst)
ccc   posakai,sosakai;expected arrival time by hypo dd.
	include 'common_arriv.inc'
	include 'common_stncod.inc'
	dimension resa(na200)
	ratsd2=rat_stn_sd2
cc	write(*,*) 'det_arriv_error'
	nc_lim=neq_branch/10
	if(nc_lim.lt.3) nc_lim=3
	k11=neq_branch+1
              keq=meqa(k11)
ccc	       call read_ok_arrival(keq)       
	   np=0
	   ns=0
	   do k=1,nobsall
	      if(wp(k).gt.0.01) np=np+1
	      if(ws(k).gt.0.01) ns=ns+1
	      nc_removea(k,1)=0  !removed number in hypo_dd
	      nc_removea(k,2)=0
c	 write(*,'(a7,2i4,4f7.2)') stnw(k),np,ns,po(k),
c     1 wp(k),so(k),ws(k)
	      enddo
	      nobsp=np
	      nobss=ns


c       write(*,195) keq,iy,im,id,ih,imin,t,phi,ram,dep,nobsp,nobss,
c     1 fmag,sd
 195	      format('# ',i5,1x,5i2,f7.1,2f8.2,f7.1,2i4,f6.1,2f7.2)

	      

		  do 1 kst=1,nobsall
	 posakai(kst)=999.9   !Expected arrivals from other events 
	 sosakai(kst)=999.9
	 res_sakai(kst,1)=99.9
	 res_sakai(kst,2)=99.9
	  nst=nstn(kst)
	 if(nst.le.0) go to 1	  
	       do 2 kps=1,2
	if(wei_crossa(k11,nst,kps).le.0.01) go to 26
ccc    Stn with P,S arrival times
		  av=0.0
		  av2=0.0
		  nc=0
		  n_all=0
		  do mmm=1,neq_branch
		     if(phi_crossa(mmm).gt.-90.0) then
		 if(wei_crossa(mmm,nst,kps).gt.0.0001) n_all=n_all+1 
		  if(wei_dd_crossa(mmm,nst,kps).ge.0.0001) then
		nc=nc+1
		resa(nc)=res_dd_crossa(mmm,nst,kps)
		  endif
		  endif
		enddo

ccc  Number of event with arrivals but removed because of large residual 
		ner=n_all-nc
		rat=float(ner)/(float(n_all)+0.1)
		if(rat.gt.0.4) nc_removea(kst,kps)=ner
cc      write(*,167) stnw(kst),(wei_crossa(m,nst,kps),m=1,neq_branch)
c		if(kps.eq.1.and.wp(kst).gt.0.01) then
c      write(*,167) stnw(kst),po(kst),wp(kst),(wei_dd_crossa(m,nst,kps),
c     1 m=1,neq_branch)
 167	 format(a7,20f7.2/7x,20f7.2/7x,20f7.2)
cc	 endif
         if(nc.le.nc_lim) go to 2  !No refference data. res_sakai=9.0
	call cal_av_rem_2sd(nc,resa,ratsd2,av,av2)
		     if(kps.eq.1) posakai(kst)=po(kst)+av
		     if(kps.eq.2) sosakai(kst)=so(kst)+av
		    res_sakai(kst,kps)=av2

c       if(kps.eq.1) write(*,'(a7,f6.1,2f7.2)') stnw(kst),deltak(kst),
c     1 po(kst),posakai(kst)
		     go to 2
ccc   Stn without P arrival ***************
 26		     continue
ccc  define posakai(theoretical arrival time)  for no P data.
ccc  res_dd_crossa(mmm,nst,kps)=
ccc  po_crossa(k11,nst)-po_crossa(mmm,nst)-trav_dd_crossa(mmm,nst,1)
ccc  posakai(kst)=po_crossa(mmm,nst)+trav_dd_crossa(mmm,nst,1)
	 call condel(ram,phi,ramobs(nst),phiobs(nst),delram)		
         call  delhh(heitk,dep,delram,th,tpram,dtdr,dtddel,kps)
		
		  do mmm=1,neq_branch
		       if(phi_crossa(mmm).gt.-90.0) then
		   if(wei_crossa(mmm,nst,kps).gt.0.01) then
		      pw=po_crossa(mmm,nst)
		      if(kps.eq.2) pw=so_crossa(mmm,nst)
	  ph=phi1_crossa(mmm)
	  ra=ram1_crossa(mmm)
	  de=dep1_crossa(mmm)
cc	  write(*,*) 'EQ=',ph,ra,de
	call conazd(ra,ph,ramobs(nst),phiobs(nst),del,az)
	heitk=altobs(nst)
        call  delhh(heitk,de,del,th,tpra,dtdr,dtddel,kps)
        if(kps.eq.1) trvp_cor(1,nst)=tpra
	if(kps.eq.2) trvs_cor(1,nst)=tpra
ccc  ***********************
	travpq=tpram-tpra+t_crossa(mmm)
	trav_dd_crossa(mmm,nst,kps)=travpq
	endif
	endif
	 enddo


		    nc=0
		    pq=0.0
		  do mmm=1,neq_branch
		       if(phi_crossa(mmm).gt.-90.0) then
		   if(wei_crossa(mmm,nst,kps).gt.0.0001) then
		      q=po_crossa(mmm,nst)
		      if(kps.eq.2) q=so_crossa(mmm,nst)
		  pq=pq+q+trav_dd_crossa(mmm,nst,kps)
		   nc=nc+1
		   resa(nc)=q+trav_dd_crossa(mmm,nst,kps)
		   endif
		   endif
		   enddo
		   if(nc.le.nc_lim) go to 2
	call cal_av_rem_2sd(nc,resa,ratsd2,av,av2)
cc		   write(*,*) stnw(kst),po(kst),nc,nc_lim
cc		   write(*,'(20f7.2)') (resa(j),j=1,nc)
		    if(kps.eq.1) posakai(kst)=pq/float(nc)
		    if(kps.eq.2) sosakai(kst)=pq/float(nc)
		     res_sakai(kst,kps)=av2 
 2		    continue
 1		    continue
cc		    write(*,'(10f8.1)') (t_crossa(m),m=1,neq_branch)

cc		    read(*,*) iok

ccc     defile iflag *********************************
		    rat_flag=2.0
		    
	do nst=1,nstobs	    
       iflag_repa(nst,1)=0
       iflag_repa(nst,2)=0
	enddo
cc	write(*,*) 'call  get_plot_ord'
	call get_plot_ord	
cc         write(*,'(20i4)') (kstn_plot_ord(j),j=1,nobsall)
cc	 write(*,*) iok
ccc    determine maximum distance for P S reading ****************
	np80=nobsp*9/10
	if(np80.le.10) np80=10
	ns80=nobss*9/10
	if(ns80.le.10) ns80=10
	   do 5 kps=1,2
	      kmax=np80
	      if(kps.eq.2) kmax=ns80
cc	      write(*,*) 'kps,kmax=',kps,kmax
	    	   nps=0  
	do 6 kk=1,nobsall
	     kst=kstn_plot_ord(kk)
	     if(kst.le.0) go to 6
	   nst=nstn(kst)
	   if(nst.le.0) go to 6
	   w=wp(kst)
	      if(kps.eq.2) w=ws(kst)
	      if(w.gt.0.01) nps=nps+1
	if(kps.eq.1) then
         sa=abs(po(kst)-posakai(kst))
	 resoc(kst,1)=po(kst)-posakai(kst)
	 if(posakai(kst).gt.900.0) resoc(kst,1)=9.9
	 if(posakai(kst).gt.900.0) go to 6
	 else
	 sa=abs(so(kst)-sosakai(kst))
	 resoc(kst,2)=so(kst)-sosakai(kst)
	 if(sosakai(kst).gt.900.0) resoc(kst,2)=9.9 
	 if(sosakai(kst).gt.900.0) go to 6
	 endif
         resmax=res_sakai(kst,kps)*rat_flag   !rat_flag=2.0
	 if(resmax.lt.err_min) resmax=err_min

ccc     Far station ************************
	 if(nps.gt.kmax.and.w.lt.0.01) go to 6   !Far stn ,No P or S , not repick.
ccc      repick stn  1) large residual,2) No picked data
	 if(w.gt.0.001) then
		if(sa.gt.resmax) iflag_repa(kst,kps)=iflag_repa(kst,kps)+1
		go to 6
		endif
ccc    No arrival time station
       if(kps.eq.1) then
       if(wp(kst).lt.0.01) iflag_repa(kst,1)=99
	resoc(kst,1)=po(kst)-posakai(kst)
	if(wp(kst).lt.0.01) resoc(kst,1)=9.9
	endif
       if(kps.eq.2) then
	if(ws(kst).lt.0.01) iflag_repa(kst,2)=99
	resoc(kst,2)=so(kst)-sosakai(kst)
	if(ws(kst).lt.0.01) resoc(kst,2)=9.9
	endif
 6	continue
 5	continue
	ns=0
cc	write(*,*) 'nobsall=',nobsall
	do 10 kk=1,nobsall
	     kst=kstn_plot_ord(kk)
	     if(kst.le.0) go to 10
	     if(ws(kst).gt.0.01) ns=ns+1
	     travp(kst)=posakai(kst)
	     travs(kst)=sosakai(kst)
	     if(posakai(kst).gt.900.0) then
		wp(kst)=0.0
		 iflag_repa(kst,1)=0
		endif
	     if(sosakai(kst).gt.900.0) then
		ws(kst)=0.0
		 iflag_repa(kst,2)=0
		endif


c      write(*,110) stnw(kst),ns,deltak(kst),po(kst),posakai(kst),
c     1   wp(kst),so(kst),sosakai(kst),ws(kst),resoc(kst,1),
c     2 res_sakai(kst,1),resoc(kst,2),res_sakai(kst,2),nc_removea(kst,1),
c     2 nc_removea(kst,2),iflag_repa(kst,1),iflag_repa(kst,2)
c      write(74,110) stnw(kst),ns,deltak(kst),po(kst),wp(kst),so(kst),
c     1 ws(kst),posakai(kst),res_sakai(kst,1),sosakai(kst),
c     2 res_sakai(kst,2),nc_removea(kst,1),nc_removea(kst,2),
c     3 iflag_repa(kst,1),iflag_repa(kst,2)
 110	   format(a7,i4,f6.1,' Ps,Theo,w=',6f7.2,' res=',4f8.3,
     3 ' ncer=',2i3,' flag=',2i3)
 10	   continue
cc	   write(*,*) 'before det_arriv_error end'
ccc	read(*,*) iok
       return
	end








	

	subroutine get_plot_ord
	include 'common_arriv.inc'
	dimension pwa(na200),nwa(na200),nstnwa(na200)
	
	   do k=1,nobsall
            pwa(k)=deltak(k)
            if(deltak(k).lt.0.01) pwa(k)=999.9
cc            if(jrec_start(k,2).le.0) pwa(k)=999.9
cc	if(phi.lt.-90.0) pwa(k)=abs(po(k))/10000.0-float(ntrg_flaga(k))
cc	    if(amp_bidou(k,1).le.0.001) pwa(k)=9999.9
	    nwa(k)=k
c         if(ipola(k).eq.'P '.or.ipola(k).eq.' P') pwa(k)=travp(k)+1000.0
c         if(ipola(k).eq.'  ') pwa(k)=travp(k)+1000.0
	    enddo

            if(phi.lt.0.0) then
	   do kst=1,nobsall
            pwa(kst)=900.0+float(nstn(kst))
            ntrgk=nphase_ps(kst)    
           if(ntrgk.gt.0) pwa(kst)=p_psa(kst,1)
           enddo
        endif



	    call plot_order(nobsall,nwa,pwa,kstn_plot_ord)
	    if(dep.gt.50.0.or.phi.le.-99.0) return
	    np=0
	    nc=0
	    ncmax=2
	    do kkk=1,nobsall
	   k=kstn_plot_ord(kkk)
	   if(wp(k).gt.0.01) np=np+1
	   if(np.ge.ncmax) go to 1
	   if(wp(k).lt.0.01) nc=nc+1
	   enddo
 1	   continue
	return
	end
  
