 	subroutine cal_hypo_dd_cross
ccc     po=po_cor(2,nst)-po_cor(1,nst) : po=po(keq)-po(meq)
ccc     po(kst) for mmm event=po_crossa(mmm,nst) 
ccc     po_cor(2,nst)-po_cor(1,nst)=travp_crossa(mmm,nst) +res_crossa(mmm,nst,1)
ccc	p_theo=travp_crossa(mmm,nst) +res_crossa(mmm,nst,1)
ccc     If wp(kst)=0  p_theo=po_cor(1,nst)+travp_crossa(mmm,nst)
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    phi1_crossa(mmm); Hypo of the first event.
	include 'common_stncod.inc'
	include 'common_arriv.inc'
	dimension wswa(11,300)
	k11=neq_branch+1
ccc	write(*,*) 'sub cal_hypo_dd_cross','k11=',k11
cc	read(*,*) iok
c        call clear_branch_arriv
	   nc=0	  
	 sd_check=0.0
	    keq=meqa(k11)
	    call read_ok_arrival(keq) !meq ;master and large mag
c	 write(*,122) keq,iy,im,id,ih,imin,t,phi,ram,dep,
c     1 nobsp,nobss,nobsall,fmag,sd_all_cross

	    do k=1,nobsall
c	       if(wp(k).gt.0.01.or.ws(k).gt.0.01)  write(*,'(a7,6f7.2)')
c     1 stnw(k),deltak(k),po(k),wp(k),so(k),ws(k)
	       enddo

cc	       read(*,*) iok
	    
cc  Clear S arrival for large event with M>5.0	    
	    if(fmag.ge.5.0) then
	       do k=1,nobsall
		  ws(k)=0.0
	       enddo
	       endif
	       


	    
         call store_branch_arriv(k11)
c	write(*,122) keq,iy,im,id,ih,imin,t,phi,ram,dep,
c     1 nobsp,nobss,nobsall,fmag,sd
	nobspk=nobsp
	       do 20 mmm=1,neq_branch
cc	       do 20 mmm=1,1
c		  write(*,*) 'read ok arrival meq=',meq
		  meq=meqa(mmm)
		  call read_ok_arrival(meq)  !meq ;master and large mag
         call store_branch_arriv(mmm)
cc		  call set_arv_diff
c		  write(*,*) 'after set_arv_diff'
c	write(*,122) keq,iy,im,id,ih,imin,t,phi,ram,dep,
c     1 nobsp,nobss,nobsall,fmag,sd
	do k=1,nobsall
	   nst=nstn(k)
	   if(nst.le.0) nst=4900
c        write(*,147) stnw(k),nst,trvp_cor(1,nst),trvp_cor(2,nst),
c     1 po(k),wp(k),so(k),ws(k),resoc(k,1),resoc(k,2)  
	enddo

         call set_arv_diff_mem(mmm)

c	 write(*,*) 'after set_arv_diff_mem'
c	write(*,122) keq,iy,im,id,ih,imin,t,phi,ram,dep,
c     1 nobsp,nobss,nobsall,fmag,sd
 122	format(i8,2x,5i3,f6.1,3f7.2,3i5,2f7.2)
c	do k=1,nobsall
c	   nst=nstn(k)
c       write(*,147) stnw(k),nst,trvp_cor(1,nst),trvp_cor(2,nst),
c     1 po(k),wp(k),so(k),ws(k),resoc(k,1),resoc(k,2)  
c	enddo



cc	 read(*,*) iok

c	write(*,163) kkk_ps_err,mmm,nobsp,nobss
 163	format('kkps,mmm,nobsp,s=',10i5)
cc	read(*,*) iok
	nps=nobsp+nobss
cc	write(*,*) 'mmm=',mmm,' nobsp,nobss=',nobsp,nobss,nps
	if(nps.le.4) then
	   	 phi_crossa(mmm)=-99.9
		 sd_crossa(mmm)=9.9
		 go to 20
		 endif
		  call change_weight_near_stn
                 call  cal_hypo_sd_dif(phi,ram,dep,sd)
cc	 write(*,*) 'mmm=',mmm,' nobsp,nobss=',nobsp,nobss,' sd=',sd

		  do k=1,nobsall
		     nst=nstn(k)
	if(wp(k).gt.0.001.or.ws(k).gt.0.01) then
c	   write(*,147) stnw(k),nst,trvp_cor(1,nst),trvp_cor(2,nst),
c     1 po(k),wp(k),so(k),ws(k),resoc(k,1),resoc(k,2)   
 147	format(a7,i4,10f7.2)
	endif
	enddo


cc	read(*,*) iok

cc	WRITE(*,1361) phi,ram,dep,nobsp,nobss,sd
 1361	format(2f8.3,f6.1,2i4,f7.2)

		 if(phi.lt.-90.0)  sd=9.9
		  call dif_hypo_grid_liner
		  sd_org=sd
		  nobsp_org=nobsp
		  nobss_org=nobss
		  nobsp_cross0a(mmm)=nobsp
		  nobss_cross0a(mmm)=nobss
		  sd_cross0a(mmm)=sd
		  call  remove_2sd_res   !Put wp(k)=-abs(wp) if large oc
		nps=nobsp+nobss
cc	write(*,*) 'after 2sd ',nobsp,nobss,nps

cc	WRITE(*,137) phi,ram,dep,nobsp,nobss,nps,sd
 137	format(' after grid liner ',2f8.3,f6.1,3i4,f7.2)

	         if(nps.le.4) then
	   	 phi_crossa(mmm)=-99.9
		 sd_crossa(mmm)=9.9
		 go to 20
		 endif
	call dif_hypo_grid_liner	  
c       write(27,107) kkk_ps_err,kkk,mmm,phi,ram,dep,t,fmag,nobsp_org,
c     1 nobss_org,nobsp,nobss,sd_org,sd,sd2_cross
c       write(*,107) kkk_ps_err,kkk,mmm,phi,ram,dep,t,fmag,nobsp_org,
c     1 nobss_org,nobsp,nobss,sd_org,sd,sd2_cross
 107    format(3i3,2f8.4,3f6.1,' PS org=',2i5,' ps=',2i4,' sd=',4f8.2)
	sd_check=sd_check+sd
		  do k=1,nobsall
c	if(wp(k).gt.0.001.or.ws(k).gt.0.01) write(*,147) stnw(k),
c     1 po(k),wp(k),so(k),ws(k)   
	enddo
        

 21	continue
        call set_theo_dif(mmm)   ! resoc(kst,1)=po(kst)-travp(kst) if abs(wp)>0.001
c	write(*,*) 'after set',nobsp,nobss,sd
c	if(mmm.eq.1) then
	
c            call  cal_hypo_sd_dif_test(phi,ram,dep,sd)
c	    write(*,*) 'call set_arv_diff_mem'
c           call set_arv_diff_mem(mmm)
c            call  cal_hypo_sd_dif_test(phi,ram,dep,sd)
cc	    read(*,*) iok
c	 write(*,136) mmm,phi,ram,dep,t,nobsp,nobss,sd,sd2_cross
c	 endif
 136	 format(i3,2f9.4,2f6.1,2i4,2f8.2)
cc	if(kkk_ps_err.eq.2) read(*,*) iok
		
cc		   read(*,*) iok
ccc     store_cross_res  **************
cc    store res to res_crossa(keq,meq,nst,kps) even large res data 
cc    sd_crossa(keq,meq)=sd	     
cc           if(sd.gt.9.0) return
c	   do 1 nst=1,nstobs
c	if(wei_crossa(k11,nst,1).le.0.01.and.wei_crossa(mmm,nst,1).
c     1 le.0.01) go to 1
c	write(*,113) stnobs(nst),po_crossa(k11,nst),wei_crossa(k11,nst,1),
c     1 po_crossa(mmm,nst),wei_crossa(mmm,nst,1),
c     2 trav_dd_crossa(mmm,nst,2),res_dd_crossa(mmm,nst,1)
c 113	format(a7,8f7.2)
c 1	continue
cc	read(*,*) iok
 20	continue
c	do mmm=1,neq_branch
c	 phi=phi_crossa(mmm)
c	 ram=ram_crossa(mmm)
c	 dep=dep_crossa(mmm)
c         call set_arv_diff_mem(mmm)
c            call  cal_hypo_sd_dif(phi,ram,dep,sd)
c	 write(27,136) mmm,phi,ram,dep,t,nobsp,nobss,sd,sd2_cross
c 136	 format(i3,2f9.4,2f6.1,2i4,2f8.2)
c	 enddo
cc	 close(27)
cc	 stop
cc	 read(*,*) iok
cc	close(27)
cc	write(*,'(10f8.3)') (sd_crossa(m),m=1,neq_branch)
	return
	end
	subroutine set_theo_dif(mmm)                   
	include 'common_stncod.inc'
	include 'common_arriv.inc'
	dimension nst_to_kst(na200)
cc	write(*,*) 'set theo diff ',phi,ram,dep,' t=',t
	nobsp=0
	nobss=0
	do nst=1,nstobs
	   nst_to_kst(nst)=0
	   enddo

        k11=neq_branch+1          
	do kst=1,nobsall
	if(wp(kst).gt.0.0001) nobsp=nobsp+1
	if(ws(kst).gt.0.0001) nobss=nobss+1
	nst=nstn(kst)
	if(nst.ge.1) nst_to_kst(nst)=kst
	enddo

        do 1 nst=1,nstobs
	wei_dd_crossa(mmm,nst,1)=0.0
	wei_dd_crossa(mmm,nst,2)=0.0
	kst=nst_to_kst(nst)
	if(kst.le.0) go to 1
	   if(wp(kst).le.0.0001.and.ws(kst).le.0.0001) go to 1
	call condel(ram,phi,ramobs(nst),phiobs(nst),del)
cc	write(*,*) 'ram,phi=',ram,phi,ramobs(nst),phiobs(nst),'del=',del
CC	WRITE(*,*) 
	 heitk=altobs(nst)
cc	 deltak(kst)=del-deltak_cor(1,nst)
        call  delhh(heitk,dep,del,th,tpq,dtdr,dtddel,1)
cc	write(*,*) 'delhh=',heitk,depq,del,th,tpq,dtdr,dtddel
	travpq=tpq-trvp_cor(1,nst)+t
        call  delhh(heitk,dep,del,th,tsq,dtdr,dtddel,2)
	travsq=tsq-trvs_cor(1,nst)+t
 100	format(a7,6f7.2)
	pq=po_crossa(k11,nst)-po_crossa(mmm,nst)
	sq=so_crossa(k11,nst)-so_crossa(mmm,nst)
	res_dd_crossa(mmm,nst,1)=po(kst)-travpq
        res_dd_crossa(mmm,nst,2)=so(kst)-travsq
	wei_dd_crossa(mmm,nst,1)=wp(kst)
	wei_dd_crossa(mmm,nst,2)=ws(kst)	
	trav_dd_crossa(mmm,nst,1)=travpq
	trav_dd_crossa(mmm,nst,2)=travsq
c	if(wei_crossa(k11,nst,1).ge.0.01.or.wei_crossa(mmm,nst,1).
c     1 ge.0.01) then
c       write(*,101) stnobs(nst),nst,po_crossa(mmm,nst),
c     1 po_crossa(k11,nst),
c    1 travpq,res_dd_crossa(mmm,nst,1),tpq,trvp_cor(1,nst)
 101	format(a7,i4,10f7.2)
c     1 wp(kst),so(kst),ws(kst),resoc(kst,1),resoc(kst,2)

cc	endif


 1	continue
	n=nobsp+nobss
cc	if(nobsp.le.2) sd=9.9
	if(n.le.5) sd=9.9 
	if(sd.gt.9.0) then
	   do kst=1,nobsall
	       resoc(kst,1)=99.9
	       resoc(kst,2)=99.9    
	enddo
	endif
		
		 phi_crossa(mmm)=phi   !phi(kkk)-phi(mmm)
		 ram_crossa(mmm)=ram
		 dep_crossa(mmm)=dep
		 sd_crossa(mmm)=sd
		  nobsp_crossa(mmm)=nobsp
		  nobss_crossa(mmm)=nobss
		  t_crossa(mmm)=t
	   return
	   end



   	subroutine dif_hypo_grid_liner
ccc	subroutine to calculate next hypocenter location and sd.
cc	input *****************************************************
ccc	phiq,ramq,depq; initial hypocenter location.
cc	dl; distance to calculate derivative in km.
ccc
ccc	output	***************************************************
ccc	sdq; standard deviation at initial hypocenter location.
ccc	phi1,ram1,dep1; least squares solution

	include 'common_arriv.inc'
	include 'common_stncod.inc'
	include 'common_channel.inc'
	common /cliner/cmat(4,9),cmatk(4,9)
	dimension xyzt(5),abf(na200,2),bbf(na200,2),cbf(na200,2)
	dimension trav_phi(na200,2),trav_ram(na200,2),trav_dep(na200,2)
	pideg=3.14159265/180.0
        sdk=99.9
	nobsp=0
	nobss=0
	do l=1,nobsall
	do kps=1,2
	resoc(l,kps)=9.9
	enddo
	if(wp(l).gt.0.0) nobsp=nobsp+1
	if(ws(l).gt.0.0) nobss=nobss+1
	enddo
        nq=nobsp+nobss
cc       if(nobsp_min_cal_hypo.le.1) write(*,*) nobsp,nobss,phiq,ramq,depq
	if(nq.le.4) go to 999
	if(phi.lt.-90.0) go to 999
	phiq=phi
	ramq=ram
	depq=dep
	phik=phiq
	ramk=ramq
	depk=depq
	dl=0.3
	call cal_hypo_sd_dif(phiq,ramq,depq,sdmin)
	sclx=cos(phiq*pideg)*111.18
	scly=111.18
	dp_grid=dl/scly
	dr_grid=dl/sclx
	dz_grid=dl
	do 888 kkmmk=1,2
	do 666 kkk=1,100
	   phiq=phik
	   ramq=ramk
	   depq=depk
	   phib=phik
	   ramb=ramk
	   depb=depk

	dp_grid=dl/scly
	dr_grid=dl/sclx
	dz_grid=dl

c	call cal_hypo_sd(phiq,ramq,depq,sdq)
c	write(*,*) 'org',phiq,ramq,depq,sdq

	phid=phiq+dp_grid
	ramd=ramq+dr_grid
	depd=depq+dz_grid
	call cal_hypo_sd_dif(phid,ramq,depq,sdq)
c	write(*,*) 'phi',phid,ramq,depq,sdq
	if(sdq.lt.sdmin) then
	   phik=phid
	   sdmin=sdq
	   endif
	do k=1,nobsall
	   trav_phi(k,1)=travp(k)
	   trav_phi(k,2)=travs(k)	   
	enddo

	call cal_hypo_sd_dif(phiq,ramd,depq,sdq)
c	write(*,*) 'ram',phiq,ramd,depq,sdq
	do k=1,nobsall
	   trav_ram(k,1)=travp(k)
	   trav_ram(k,2)=travs(k)	   
	enddo
	if(sdq.lt.sdmin) then
	   ramk=ramd
	   sdmin=sdq
	   endif
	call cal_hypo_sd_dif(phiq,ramq,depd,sdq)
c	write(*,*) 'dep',phiq,ramq,depd,sdq
	do k=1,nobsall
	   trav_dep(k,1)=travp(k)
	   trav_dep(k,2)=travs(k)	   
	enddo
	if(sdq.lt.sdmin) then
	   depk=depd
	   sdmin=sdq
	   endif
	call cal_hypo_sd_dif(phiq,ramq,depq,sdq)
	if(sdq.lt.sdmin) then
	   phik=phiq
	   ramk=ramq
	   depk=depq
	   sdmin=sdq
	   endif

	do k=1,nobsall
	abf(k,1)=(trav_ram(k,1)-travp(k))/dr_grid  !dtdram-unit  deg
	abf(k,2)=(trav_ram(k,2)-travs(k))/dr_grid   !dtdram-unit
	bbf(k,1)=(trav_phi(k,1)-travp(k))/dp_grid		!dtdphi-unit
	bbf(k,2)=(trav_phi(k,2)-travs(k))/dp_grid		!dtdphi-unit
	cbf(k,1)=(trav_dep(k,1)-travp(k))/dz_grid		!dtddep-unit
	cbf(k,2)=(trav_dep(k,2)-travs(k))/dz_grid		!dtddep-unit
	enddo
c	n4=3
c	n9=7
	n4=4
	n9=9

      do 33 k=1,n4
      do 34 m=1,n9
      cmat(k,m)=0.0
   34 continue
      mm=k+n4+1
      cmat(k,mm)=1.0
   33 continue

	wtot=0.0
      do  kps=1,2
	nst=0
	do 48 l=1,nobsall
	if(kps.eq.1) w=wp(l)
	if(kps.eq.2) w=ws(l)
	if(w.le.0.0) go to 48
	nst=nst+1
	oc=resoc(l,kps)
	if(oc.gt.100.0.or.oc.lt.-100.0) go to 48
      xyzt(1)=abf(l,kps)	!dt/dram(deg)
      xyzt(2)=bbf(l,kps)		!dt/dphi(deg)
      xyzt(3)=cbf(l,kps)		!dt/dz
c      xyzt(4)=oc
      xyzt(4)=1.0
      xyzt(5)=oc


	sdq=sdq+oc*oc*w
	wtot=wtot+w

c	write(*,'(a7,i3,5f8.2)') stnw(l),kps,(xyzt(m),m=1,5)

      do 39 k=1,n4
      do  m=1,n4+1
      qx=xyzt(k)*xyzt(m)*w
      cmat(k,m)=cmat(k,m)+qx
	enddo
   39 continue
   48   continue
	enddo      

      call liner(n4,n9)

      dr=cmat(1,n4+1)
      dp=cmat(2,n4+1)
      dz=cmat(3,n4+1) 



      ram1=ramq+dr
      phi1=phiq+dp
      dep1=depq+dz    
	if(dep1.lt.-2.0) dep1=-2.0

	call cal_hypo_sd_dif(phi1,ram1,dep1,sdq)
c	write(*,100) kkk,dp,dr,dz,phi1,ram1,dep1,sdq,sdmin
 100	format('ans=',i3,10f9.3)
	if(sdq.lt.sdmin) then
	   phik=phi1
	   ramk=ram1
	   depk=dep1
	   sdmin=sdq
	   endif
           nst4=nobsp+nobss-4
           if(nst4.lt.1) nst4=1
           sd2=sqrt(sdq/float(nst4))
      dphi=sdq*sqrt(abs(cmat(1,n4+2)))
      dram=sdq*sqrt(abs(cmat(2,n4+3)))
      ddep=sdq*sqrt(abs(cmat(3,n4+4)))
      dt=sdq*sqrt(abs(cmat(4,n4+5)))
	if(n4.eq.3) dt=0.0
      if(dt.gt.9.9) dt=9.9     
	if(dphi.gt.9.9) dphi=9.9
	if(dram.gt.9.9) dram=9.9
	if(ddep.gt.99.9) ddep=99.9   
	 dx=(ram1-ramq)*sclx
	 dy=(phi1-phiq)*scly
	 dz=dep1-depq
	 dlk=dl
	 dl=sqrt(dx*dx+dy*dy+dz*dz)
cc	 write(*,*) kkk,phik,ramk,depk,sdmin
	 if(dl.lt.0.2) go to 667
	 if(dl.gt.dlk) dl=dlk
	 if(dl.lt.1.0) dl=1.0
	  sa=abs(phib-phik)+abs(ramb-ramk)+(depb-depk)
	  if(sa.lt.0.0001) go to 667
 666	 continue
 667	 continue
	 phi=phik
	 ram=ramk
	 dep=depk
	call cal_hypo_sd_dif(phi,ram,dep,sd)
	 do k=1,nobsall
	   q=abs(resoc(k,1))
	   if(q.gt.3.5) wp(k)=-abs(wp(k))
	   q=abs(resoc(k,2))
	   if(q.gt.3.5) ws(k)=-abs(ws(k))   
	   enddo



 888	 continue

cc	read(*,*) iok
     
      return
  999 continue
      sdk=99.9
      dt=0.0
      dphi=0.0
      dram=0.0
      ddep=0.0
	phi=-99.9
	ram=0.0
	dep=0.0
      return
      end
	 subroutine  ini_multi_eq
	include 'common_stncod.inc'
	include 'common_arriv.inc'
	call initial_multi_hypo  !phieri,rameri,deperi
cc	write(*,*) 'auto hypo=',phi,ram,dep
cc	write(*,*) 'initial hypo=',phieri,rameri,deperi
c	write(*,'(10f8.3)') (phi_crossa(m),m=1,neq_branch)
c	write(*,'(10f8.3)') (ram_crossa(m),m=1,neq_branch)
cc	write(*,'(10f8.3)') (dep_crossa(m),m=1,neq_branch)
         call cal_all_cross_sd(phieri,rameri,deperi,sd_ini)
c	 write(*,100) phieri,rameri,deperi,sd_ini
 100	 format('initial hypo=',2f9.4,f6.1,' sd2=',f8.1)
ccc     test
cc	 call cal_all_cross_sd2(sd2_tot)
c	 write(*,*) 'sd2_tot=',sd2_tot
c	 write(*,'(10f8.3)') (sd_crossa(m),m=1,neq_branch)
c	 read(*,*) iok
c	 write(*,*) 'av2phi=',av2phi,av2ram,av2dep
c	 phi1=phieri-av2phi
c	 phi2=phieri+av2phi
c	 ram1=rameri-av2ram
c	 ram2=rameri+av2ram
c	 dep1=deperi-av2dep
c	 dep2=deperi+av2dep
	 sdmin=sd_ini
	 dp=av2phi/3.0
	 dr=av2ram/3.0
	 dz=av2dep/3.0
	 do k=1,10
         call det_xyz_min(phieri,rameri,deperi,dp,dr,dz,
     1 phiq,ramq,depq,sd2)
c	 write(*,*) 'det_xyz_min ',phieri,rameri,deperi,dp,dr,dz,
c     1 phiq,ramq,depq,sd2
	 iw=0
	 if(sd2.gt.sdmin) iw=1
c	 if(phiq.lt.phi1.or.phiq.gt.phi2) iw=1
c	 if(ramq.lt.ram1.or.ramq.gt.ram2) iw=1
c	 if(depq.lt.dep1.or.depq.gt.dep2) iw=1
	 if(iw.le.0) then
	 phieri=phiq
	 rameri=ramq
	 deperi=depq
	 sdmin=sd2
	 else
	    call grid_xyz_min(phieri,rameri,deperi,dp,dr,dz,
     1 phiq,ramq,depq,sd2) 	
c	    write(*,*) 'grid_xyz_',phieri,rameri,deperi,dp,dr,dz,
c     1 phiq,ramq,depq,sd2
	 phieri=phiq
	 rameri=ramq
	 deperi=depq
	 sdmin=sd2
	    dp=dp/2.0
	    dr=dr/2.0
	    dz=dz/2.0
	    if(dr.lt.0.0005) go to 15
	    endif
 101	 format('located hypo=',2f9.4,f6.1,' dp=',f7.4,' sd2=',2f8.3)
	 enddo
 15	 continue
cc	 write(*,101) phieri,rameri,deperi,dp,sd2,sdmin
	 phi=phieri
	 ram=rameri
	 dep=deperi
	 sd_all_cross=sdmin
	 sd=sdmin
c	 sd_hyp=0.0
c	 do mmm=1,neq_branch
c	  call set_arv_diff_mem(mmm)
c        call cal_hypo_sd_dif(phi_crossa(mmm),ram_crossa(mmm),
c     1 dep_crossa(mmm),sd)
c	sd_hyp=sd_hyp+sd2_cross
c	write(*,167) mmm,phi_crossa(mmm),ram_crossa(mmm),
c     1 dep_crossa(mmm),nobsp,nobss,sd,sd2_cross
c	enddo
c 167	format('check ',i3,2f9.4,f6.1,2i5,' sd=',2f8.2)
c	write(*,*) 'sd_hyp=',sd_hyp
	 return
	 end

	subroutine initial_multi_hypo
ccc     po=po_cor(2,nst)-po_cor(1,nst) : po=po(keq)-po(meq)
ccc     res_crossa(keq,meq,nst,1)=(po(k)-travp(k))
ccc     res_crossa(keq,meq,nst,1)=po(keq)-po(meq)-travp
	include 'common_stncod.inc'
	include 'common_arriv.inc'
	dimension phia(1000),rama(1000),depa(1000)
	nc=0
	phieri=-99.9
	rameri=0.0
	deperi=0.0
	do 1 mmm=1,neq_branch
	   q=phi_crossa(mmm)
	   if(q.lt.-90.0) go to 1
	      nc=nc+1
	    phia(nc)= phi_crossa(mmm)
	    rama(nc)= ram_crossa(mmm)
	    depa(nc)= dep_crossa(mmm)
 1	    continue
	      if(nc.le.0) return
	      ratsd2=2.0
	call cal_av_rem_2sd(nc,phia,ratsd2,phieri,av2phi)
	call cal_av_rem_2sd(nc,rama,ratsd2,rameri,av2ram)
	call cal_av_rem_2sd(nc,depa,ratsd2,deperi,av2dep)

	if(av2phi.le.0.001)  av2phi=0.001
	if(av2ram.le.0.001)  av2ram=0.001
	if(av2dep.le.0.3) av2dep=0.3
c	write(*,*) 'phi=',nc
c	write(*,'(20f8.4)') (phia(m),m=1,neq_branch),phieri
c        write(*,*) 'av2phi=', av2phi,av2ram,av2dep
	
cc	read(*,*) iok

	return
	end
   	subroutine cal_all_cross_sd(phik,ramk,depk,sd2_tot)
ccc     po=po_cor(2,nst)-po_cor(1,nst) : po=po(keq)-po(meq)
	include 'common_arriv.inc'
	include 'common_stncod.inc'
	 sdt2=0.0
	 wei=0.0
	       do 20 mmm=1,neq_branch
		  if(phi_crossa(mmm).le.-90.0) go to 20
		  call set_arv_diff_mem(mmm)
cc	  write(*,*) 'after set arv diff',nobsall,nobsp,nobss
	   call change_weight_near_stn	  
cc	     write(*,*) 'cal_hypo_sd_dif',phik,ramk,depk
                 call  cal_hypo_sd_dif(phik,ramk,depk,sd)
		        call  remove_2sd_res  
          call  cal_hypo_sd_dif(phik,ramk,depk,sd)
cc	 write(*,178) mmm,phik,ramk,depk,nobsp,nobss,sd,sd2_cross
 178	 format(i3,2f9.4,f6.1,2i4,' sd=',2f8.2)
ccc     sd_net:sum of res*res
		sdt2=sdt2+sd2_cross
		wei=wei+wei_cross
 20	continue
	 sd2_tot=sqrt(sdt2/(wei+0.00001))
cc	 write(*,*) ' sd2_tot=', sd2_tot
	return
	end
	subroutine cal_all_cross_sd2(sd2_tot)
ccc     po=po_cor(2,nst)-po_cor(1,nst) : po=po(keq)-po(meq)
	include 'common_arriv.inc'
	include 'common_stncod.inc'
	 sdt2=0.0
	 wei=0.0
	       do 20 mmm=1,neq_branch
		  if(phi_crossa(mmm).le.-90.0) go to 20
		  call set_arv_diff_mem(mmm)
                  call change_weight_near_stn
       call  cal_hypo_sd_dif(phi_crossa(mmm),ram_crossa(mmm),
     1 dep_crossa(mmm),sd_crossa(mmm))
cc           call  remove_2sd_res 
      call  cal_hypo_sd_dif(phi_crossa(mmm),ram_crossa(mmm),
     1 dep_crossa(mmm),sd_crossa(mmm))
  
cc	  write(*,*) 'after set arv diff',nobsall,nobsp,nobss
		  do kst=1,nobsall
cc	     write(*,125) stnw(kst),po(kst),wp(kst),so(kst),ws(kst)
 125	     format(a7,10f7.2)
	     enddo
cc	     write(*,*) 'cal_hypo_sd_dif',phik,ramk,depk
                 call  cal_hypo_sd_dif(phi_crossa(mmm),ram_crossa(mmm),
     1 dep_crossa(mmm),sd_crossa(mmm))
cc	 write(*,178) mmm,phik,ramk,depk,nobsp,nobss,sd,sd2_cross
 178	 format(i3,2f9.4,f6.1,2i4,' sd=',2f8.2)
ccc     sd_net:sum of res*res
		sdt2=sdt2+sd_crossa(mmm)
		wei=wei+wei_cross
 20	continue
	 sd2_tot=sqrt(sdt2/(wei+0.00001))
cc	 write(*,*) ' sd2_tot=', sd2_tot
	return
	end
	subroutine cal_av_rem_2sd(nn,resa,ratsd2,av,av2)
	dimension resa(nn),resa2(100000)
	av=0.0
	av2=0.0
	if(nn.le.0) return
	do j=1,nn
	   av=av+resa(j)
	   enddo
	   av=av/float(nn)
	   do j=1,nn
	      q=resa(j)-av
	      av2=av2+q*q
	      enddo
	   av2=sqrt(abs(av2/float(nn)))
	   amp=av2*ratsd2
cc	   write(*,*) 'amp=',amp,av2
	   z=0.0
	   zz=0.0
	   nc=0
	do j=1,nn
	   sa=abs(resa(j)-av)
	   if(sa.lt.amp) then
	      nc=nc+1
	   z=z+resa(j)
	   zz=zz+resa(j)*resa(j)
	   resa2(nc)=resa(j)
	   endif
	   enddo
	   if(nc.ge.1) then
	      av=z/float(nc)
	      z=0.0
	      do j=1,nc
		 q=resa2(j)-av
		 z=z+q*q
		 enddo
	      av2=sqrt(abs(z/float(nc)))
	      endif
c	      write(*,*) (resa(j),j=1,nn)
c	      write(*,*) 'nc=',nc,nn
c	      write(*,*) 'av,av2=',av,av2
c	      read(*,*) iok
	      return
	      end
         subroutine change_weight_near_stn
	 include 'common_arriv.inc'
	 include 'common_stncod.inc'
        dimension dela(na200),ich_ord(na200),nwa(na200),ndelord(na200)
	 if(nobsp.le.10) return
	 
	 do 1 k=1,nobsall
	    deltak(k)=9999.9
	    nst=nstn(k)
	    if(nst.le.0) go to 1
	  call conazd(rap,phi,ramobs(nst),phiobs(nst),deltak(k),az)
 1	  continue
	     call del_order(nobsall,deltak,ich_ord)
	     depk=dep
	     if(depk.lt.10.0) depk=10.0
	     delk=depk
	     nc=0
          do 5 kkk=1,nobsall
	     k=ich_ord(kkk)
	    if(k.le.0) go to 5
	  if(deltak(k).lt.delk) go to 5
	       if(wp(k).gt.0.1) then
		  nc=nc+1
		  w=1.0
	       if(nc.gt.10) w=0.3
	       if(nc.gt.20) w=0.1
	       if(nc.gt.30) w=0.03
		  wp(k)=wp(k)*w
		  ws(k)=ws(k)*w
		  endif
 5		  continue
		  return
		  end

    

	subroutine cal_hypo_sd_dif(phiq,ramq,depq,sdq)
	include 'common_stncod.inc'
	include 'common_arriv.inc'
	tt=0.0
	w=0.0
cc	write(*,*) 'cal_hypo_sd_dif ',phiq,ramq,depq
	do 1 kst=1,nobsall
	nst=nstn(kst)
	if(nst.le.0) go to 1
	travp(kst)=0.0
	travs(kst)=0.0
	resoc(kst,1)=0.0
	resoc(kst,2)=0.0
	call condel(ramq,phiq,ramobs(nst),phiobs(nst),del)
cc	write(*,*) 'ram,phi=',ramq,phiq,ramobs(nst),phiobs(nst),'del=',del
CC	WRITE(*,*) 
	 heitk=altobs(nst)
	 deltak(kst)=del-deltak_cor(1,nst)
	 if(wp(kst).gt.0.000001) then
        call  delhh(heitk,depq,del,th,tpq,dtdr,dtddel,1)
cc	write(*,*) 'delhh=',heitk,depq,del,th,tpq,dtdr,dtddel
	travp(kst)=tpq-trvp_cor(1,nst)
	tt=tt+(po(kst)-travp(kst))*wp(kst)
	w=w+wp(kst)
	sa=po(kst)-travp(kst)
	tsa=trvp_cor(2,nst)-trvp_cor(1,nst)
cc	write(*,*) heitk,depq,' del=',del,'tpq=',tpq
	res=po(kst)-tsa
c	write(*,146) stnw(kst),del,po_cor(2,nst),po_cor(1,nst),
c     1  trvp_cor(2,nst),trvp_cor(1,nst),tpq,po(kst),tsa,res
 146	format(a7,f6.1,8f7.2)
	endif
	 if(ws(kst).gt.0.000001) then
        call  delhh(heitk,depq,del,th,tsq,dtdr,dtddel,2)
	travs(kst)=tsq-trvs_cor(1,nst)
	tt=tt+(so(kst)-travs(kst))*ws(kst)
	w=w+ws(kst)
	endif


 1	continue
c	write(*,*) 'w=',w
	if(w.le.0.01) return
	t=tt/w
	sdq=0.0
c	write(*,*) 't=',t
	w=0.0
	do 2 kst=1,nobsall
	nst=nstn(kst)
	if(nst.le.0) go to 2
	if(wp(kst).gt.0.000001) then
	   travp(kst)=travp(kst)+t
	   res=po(kst)-travp(kst)
	   resoc(kst,1)=res
	   sdq=sdq+res*res*wp(kst)
	   w=w+wp(kst)
cc	   write(*,'(a7,3hres,3f7.2)') stnw(kst),po(kst),travp(kst),res
	   endif
	if(ws(kst).gt.0.000001) then
	   travs(kst)=travs(kst)+t
	   res=so(kst)-travs(kst)
	   sdq=sdq+res*res*ws(kst)
	   w=w+ws(kst)
	   resoc(kst,2)=res
	   endif

 2	   continue
	   sd_net=sdq
	   sdq=sqrt(sdq/w)
	   sd2_cross=sd_net
	   wei_cross=w
cc	   read(*,*) iok

	   return
	   end
      subroutine remove_2sd_res
	include 'common_arriv.inc'
	include 'common_stncod.inc'
	sdw=sd
	if(sdw.le.0.15) sdw=0.15
	if(sdw.gt.1.5) sdw=1.5
	respmax=sdw*2.0
	ressmax=sdw*2.5


	do k=1,nobsall
	   if(wp(k).gt.0.001) then
	      resp=abs(resoc(k,1))
	      if(resp.gt.respmax) wp(k)=-abs(wp(k))
	      endif
	   if(ws(k).gt.0.001) then
	      ress=abs(resoc(k,2))
	      if(ress.gt.ressmax) ws(k)=-abs(ws(k))
	      endif
	      enddo
	      return
	      end
	subroutine det_xyz_min(x0,y0,z0,dx,dy,dz,x,y,z,sdmin)
ccc     x0,y0,z0m  initial hypocenter
ccc     x,y,z :solution
ccc     y=ax*x+bx+c =a(x-b/a/2)**2+q  x=b/a/2
	dimension funxx(3),funyy(3),funxyz(3,3)
	funxyz(1,1)=x0
	funxyz(2,1)=y0
	funxyz(3,1)=z0
	funxyz(1,2)=x0+dx
	funxyz(2,2)=y0+dy
	funxyz(3,2)=z0+dz
	funxyz(1,3)=x0+dx*2.0
	funxyz(2,3)=y0+dy*2.0
	funxyz(3,3)=z0+dz*2.0
cc	write(*,*) 'det_xyz_min '
	do  k=1,3
cc	   write(*,'(i4,3f12.7)') k,(funxyz(m,k),m=1,3)
	   enddo

	   do m=1,3
	     funxx(m)=funxyz(1,m)
	     call cal_xyz_val(funxx(m),y0,z0,funyy(m))
	     enddo
	      call fun_err_min(funxx,funyy,x)
cc	      write(*,*) 'phi=',x
	   do m=1,3
	     funxx(m)=funxyz(2,m)
	     call cal_xyz_val(x0,funxx(m),z0,funyy(m))
ccc	     write(*,*) 'cal_xyz_val ',x0,funxx(m),z0,funyy(m)
	     enddo
	      call fun_err_min(funxx,funyy,y)
cc	      write(*,*) 'ram=',(funxx(m),m=1,3),y
	   do m=1,3
	     funxx(m)=funxyz(3,m)
	     call cal_xyz_val(x0,y0,funxx(m),funyy(m))
	     enddo
	     call fun_err_min(funxx,funyy,z)
	     if(z.gt.700.0) z=700.0
	     if(z.lt.-1.0) z=-1.0
cc	      write(*,*) 'z=',z
	       call cal_xyz_val(x,y,z,sdmin)

	      return

	      end
	subroutine cal_xyz_val(phik,ramk,depk,sd2_tot)
	 call cal_all_cross_sd(phik,ramk,depk,sd2_tot)
        return
	end
	  subroutine  grid_xyz_min(phi,ram,dep,dp,dr,dz,
     1 phimin,rammin,depmin,sdmin) 
	  sdmin=999999.9
       ngrid=5
	nm=ngrid/2+1
	grid=float(ngrid)
	dpq=dp/grid*2.0
	drq=dr/grid*2.0
	dzq=dz/grid*2.0
	do kp=1,ngrid
	   do kr=1,ngrid
	      do kz=1,ngrid
		 phiq=phi-(kp-nm)*dpq
		 ramq=ram-(kr-nm)*drq
		 depq=dep-(kz-nm)*dzq
		 if(depq.lt.-0.1) depq=-0.1
         call cal_xyz_val(phiq,ramq,depq,sdq)
	 if(sdq.lt.sdmin) then
	    sdmin=sdq
	    phimin=phiq
	    rammin=ramq
	    depmin=depq
	    endif
	    enddo
	    enddo
	    enddo
	    return
	    end
	subroutine fun_err_min(funxx,funyy,b)
ccc	function to get position of error minimum.
	dimension funxx(3),funyy(3)
ccc	determine x giving minimum of second-degree function.
ccc	y=a(x-b)*(x-b)+c	
	x1=funxx(1)
	x2=funxx(2)
	x3=funxx(3)
	y1=funyy(1)
	y2=funyy(2)
	y3=funyy(3)
cc	if(x1.eq.x2.or.x2.eq.x3)  then
cc	write(*,*) 'argument error in sub. fun_err_min',x1,x2,x3
cc	read(*,*) iok
cc	endif
	dy1=(y2-y1)/(x2-x1)
	dy2=(y3-y2)/(x3-x2)
	q=x1-x3
	c1=(dy1-dy2)/q
	c2=((x1+x2)*dy2-(x2+x3)*dy1)/q
	if(c1.le.0.0) go to 1
	b=-c2/c1/2.0
	go to 2
    1   b=x1+x1-x2
	if(y3.lt.y1) b=x3+x3-x2
    2   continue
	return
	end

 

	 subroutine  check_hypo_td
	include 'common_stncod.inc'
	include 'common_arriv.inc'
	dimension phiqa(200),depqa(200),ramqa(200),sdqa(200)
ccc     determine initial hypo by averaging dd hypo
cc	write(*,*) 'check_hypo_td ',nobsall
	if(nobsall.le.0) go to 999
c		 write(*,122) keq,iy,im,id,ih,imin,t,phi,ram,dep,
c     1 nobsp,nobss,nobsall,fmag,sd_all_cross
 122		 format(i6,5i3,f6.1,2f9.4,f6.1,3i5,f6.1,f7.2)
		 del=0.0
	    do k=1,nobsall
c	       write(*,'(a7,6f7.2)') stnw(k),po(k),so(k),deltak(k)
	       del=del+deltak(k)
	    enddo
	    delav=del/float(nobsall)
c	    write(*,*) 'delav=',delav,del_lim_sd,neq_branch
c	    write(*,*) 'P'
c	    write(*,'(10i6)') (nobsp_crossa(mmm),mmm=1,neq_branch)
c	    write(*,*) 'S'
c	    write(*,'(10i6)') (nobss_crossa(mmm),mmm=1,neq_branch)
c	    write(*,*) 'SD'
cc	    write(*,'(10f7.2)') (sd_crossa(mmm),mmm=1,neq_branch)
c	    write(*,'(10f9.4)') (phi_crossa(mmm),mmm=1,neq_branch)
c	    write(*,'(10f9.4)') (ram_crossa(mmm),mmm=1,neq_branch)
c	    write(*,'(10f9.4)') (dep_crossa(mmm),mmm=1,neq_branch)
ccc     Use only template eq with low SD
	    nc=0
	    del3=del_lim_sd*3.0*delav/del_lim_sd
	
	    
	    
	    do mmm=1,neq_branch
	       call conazd(ram,phi,ram_crossa(mmm),phi_crossa(mmm),del,az)
	      if(del.gt.del3) sd_crossa(mmm)=9.9
	       if(sd_crossa(mmm).le.9.0) then
		  nc=nc+1
		  sdqa(nc)=sd_crossa(mmm)
	       endif
	    enddo
	    rat=float(nc)/float(mmm)
	  if(nc.le.1) go to 999  
	  call cal_av_rem_2sd(nc,sdqa,ratsd2,sdav,sdav2)
	  sdavk=sdav
	  if(sdav.lt. err_min) sdav=err_min
	  sdk=sdav*rat_stn_sd2 
	    do mmm=1,neq_branch
	       if(sd_crossa(mmm).gt.sdk) then
c		  write(*,*) 'remove set=',mmm,sd_crossa(mmm),sdk,err_min
		  sd_crossa(mmm)=9.9
		  phi_crossa(mmm)=-99.9
	       endif
	       enddo
c	       write(*,*) 'sdk=',sdk,sdavk,sdav
c	    write(*,'(10f7.2)') (sd_crossa(mmm),mmm=1,neq_branch)

cc	    read(*,*) iok

	    nc=0
	    nok=0
	    sdq=0.0
	    sdk=sd_remove
	    dlk=dl_remove
	    delk=del_lim_sd
	    delk2=delk
	    if(delav.gt.delk2)  sdk=sd_remove*delav/del_lim_sd
	    if(delav.gt.delk2)  dlk=dl_remove*delav/del_lim_sd
	    if(fmag.gt.5.0) then
	       sdk=sdk*2.0
	       dlk=dlk*5.0
	       endif

	    
	    do k=1,neq_branch
	       if(sd_crossa(k).gt.sdk) nc=nc+1
	    if(sd_crossa(k).lt.9.0) then
		  sdq=sdq+sd_crossa(k)
		  nok=nok+1
		  phiqa(nok)=phi_crossa(k)
		  ramqa(nok)=ram_crossa(k)
		  depqa(nok)=dep_crossa(k)
		  endif
	    enddo
	    if(nok.le.1) go to 999
	    sq=sqrt(float(nok-1))
	     rat=float(nok)/float(neq_branch)
	    sdq=sdq/float(nok)
	    call cal_av_rem_2sd(nok,phiqa,ratsd2,ph,dph)
	    call cal_av_rem_2sd(nok,ramqa,ratsd2,rh,dra)
	    call cal_av_rem_2sd(nok,depqa,ratsd2,dh,dz)
c	    write(*,'(10f8.3)') (phiqa(j),j=1,nok)
c	    write(*,'(10f8.3)') (ramqa(j),j=1,nok)
c	    write(*,'(10f8.3)') (depqa(j),j=1,nok)
c           write(*,*) 'ph=',ph,dph,rh,dra,dh,dz
	    dphi=dph/sq
	    dram=dra/sq
	    ddep=dz/sq
cc	    write(*,*) dph,dra,dz,sq,' dphi,dram,ddep=',dphi,dram,ddep
c	    read(*,*) iok
	    pideg=3.14159/180.0
	    dy=dph*111.18
	    dx=dra*111.18*cos(ph*pideg)
cc	    write(*,*) 'dx,dy,dz=',dx,dy,dz,' sd,rat=',sdq,rat
	    dl=sqrt(dx*dx+dy*dy+dz*dz/4.0)
cc	    write(*,*) 'sdq,dl,rat=',sdq,dl,rat,sdk,dlk,rat_err_eq
	    if(rat.gt.0.7) then
	       sdq=sdq/1.5
	       dl=dl/1.5
	       endif
	       if(dep.gt.100.0) then
	       sdq=sdq/1.5
	       dl=dl/1.5
	       rat=rat*1.5
	    endif
	    
	    if(sdq.gt.sdk) go to 999
	    if(dl.gt.dlk) go to 999
	    if(rat.lt.rat_err_eq) go to 999
	    sd_all_cross=sdq
	    phieri=phi
	    rameri=ram
	    deperi=dep
	    go to 888
 999	    continue
	      phieri=-99.9
	      rameri=0.0
	      deperi=0.0
	      write(*,168)  iy,im,id,ih,imin,t,phi,ram,dep,
	1 nobsp,nobss,nobsall,fmag,delav,dl,rat,sdq,sd_all_cross
	      sd_all_cross=9.9
 168	      format('removed event ',5i2,f6.1,2f10.4,f6.1,3i5,6f7.2)

 888	      continue
c	      	open(56,file='check_hypo.dat',access='append')
c	      write(56,167)  iy,im,id,ih,imin,t,phi,ram,dep,
c	1 nobsp,nobss,nobsall,fmag,delav,dl,dlk,rat,rat_err_eq,sdq,sdk,
c	1	     sd_all_cross
c	    	      write(*,167)  iy,im,id,ih,imin,t,phi,ram,dep,
c	1 nobsp,nobss,nobsall,fmag,delav,dl,dlk,rat,rat_err_eq,sdq,sdk,
c	1	     sd_all_cross  
 167	    format(i5,4i2,f6.1,3f7.2,3i4,2f6.1,3(2f8.2,2x),4f8.2)
c	    close(56)
cc	    if(sd_all_cross.gt.9.0) read(*,*) iok
	   return
	   end
