



	subroutine inv_matrix_data(nn,mm,ndata,qmat,xa,errora,sd)
	dimension qmat(nn,ndata),xa(nn),errora(nn)
	dimension cmat(100,200)
ccc     initialize
	      do k=1,nn-1
		 do m=1,mm
		  cmat(k,m)=0.0
		  enddo
		  cmat(k,nn+k)=1.0
		  enddo
ccc     calculate matrix.
		  do 10 kkk=1,ndata
		  do k=1,nn-1
		     do m=1,nn
			cmat(k,m)=cmat(k,m)+qmat(k,kkk)*qmat(m,kkk)
			enddo
			enddo
 10			continue
			call liner4(nn-1,mm,cmat)
			sd=0.0
			do k=1,nn-1
			   xa(k)=cmat(k,nn)
			   errora(k)=abs(cmat(k,nn+k))
			   enddo
			   xa(nn)=-1.0
			   do j=1,ndata
			      z=0.0
			      do m=1,nn
				 z=z+xa(m)*qmat(m,j)
				 enddo
				 sd=sd+z*z
ccc      if(j.le.10) write(*,*) j,qmat(2,j),qmat(3,j),z
				 enddo
				 q=float(ndata-nn-1)
				 if(q.le.0.5) q=1.0
				 sd=sqrt(sd/q)
				 do k=1,nn
				 errora(k)=sqrt(errora(k))*sd
				 enddo
	      return
	      end

      subroutine liner4(m,nn,cmat)
      dimension cmat(100,200)
      do 10 k=1,m
      w=cmat(k,k)
      if(w.eq.0.0) w=0.001
      kk=k+1
      do 20 j=kk,nn
      cmat(k,j)=cmat(k,j)/w
   20 continue
      do 30 i=1,m
      if(i.eq.k) go to 30
      w=cmat(i,k)
      do 40 j=kk,nn
      cmat(i,j)=cmat(i,j)-w*cmat(k,j)
   40 continue
   30 continue
   10 continue
      return
      end

      subroutine delth2(r1,del,th,t,dtdr,dtddel,kps)
      common /cvel/al(2,2),be(2,2),r,rm,ntype_velmodel
      va1=al(1,kps)
      va2=al(2,kps)
      be1=be(1,kps)
      be2=be(2,kps)
      v1=va1*(r1/r)**be1
      if(r1.lt.rm) v1=va2*(r1/rm)**be2
      vd=va1*(rm/r)**be1
      q1=r1/v1
      qm=rm/va2
      qd=rm/vd
      qs=r/va1
      a1=1.0-be1
      a2=1.0-be2
      delt=del/r
      call hadel(qm,a1,qs,qd,crdel)
      if(r1.lt.rm) go to 1
      qs1=q1-qs
      qsd=qs-qd
      s=sin(a1*delt/2.0)
      tg=sqrt2(qs1*qs1+4.0*qs*q1*s*s)/a1
      call hadel(qm,a1,q1,qd,d1)
      d2=delt-d1-crdel
      if(d2.lt.0.0) go to 2
ccc	travel time of Pn	***********************************
      do 3 k=1,30
      dth=d2/2.0*a2
      hp=cos(dth)*qm
      call hadel(hp,a1,q1,qd,d11)
      call hadel(hp,a1,qs,qd,d12)
      dsa=abs(delt-d2-d11-d12)*r
      d2=delt-d11-d12
      if(dsa.lt.0.1) go to 4
    3 continue
      go to 999
    4 continue
      s1=sin(d11*a1/2.0)
      s2=sin(d12*a1/2.0)
      tq=(sqrt2((q1-qd)*(q1-qd)+4.0*q1*qd*s1*s1)+sqrt2(qsd*qsd+
     1 4.0*qs*qd*s2*s2))/a1+2.0*qm*sin(d2*a2/2.0)/a2
      if(tq.gt.tg) go to 2
    7 t=tq
      s=hp/q1
      c=sqrt2(1.0-s*s)
      th=atan3(s,c)
      go to 998
    2 th=atan3(sin(a1*delt),q1/qs-cos(a1*delt))
      sg=del/(abs(r-r1)+0.01)
      if(sg.lt.0.01.and.th.lt.1.5) th=3.14159265
      sg=sin(th)-qd/q1
      if(th.lt.1.57.and.sg.lt.0.0) then
	if(d2.gt.0.0) go to 7		!r1>radius
	endif
      delmax=3.14/a1
      if(delt.gt.delmax) go to 7
      t=tg
      go to 998
    1 continue
      qds2=(qd-qs)*(qd-qs)
      qm2=(q1-qm)*(q1-qm)
      d1=delt-crdel
      if(d1.gt.0.0) go to 11
	if(r.eq.r1) go to 11
      d1=delt*(rm-r1)/(r-r1)
      dh=rm-r1
      if(dh.gt.0.2) go to 11
      d2=delt
      thd=atan3(sin(a1*delt),qd/qs-cos(a1*delt))
      s=sin(d2*a1/2.0)
      t=sqrt2((qs-qd)*(qs-qd)+4.0*qs*qd*s*s)/a1
      s=v1/vd*sin(thd)*rm/r1
      c=-sqrt2(abs(1.0-s*s))
      th=atan3(s,c)
      t=t-dh*c/v1
      go to 998
   11 continue
      do 5 k=1,30
      th=atan3(sin(a2*d1),q1/qm-cos(a2*d1))
           hp=sin(th)*q1
      hp2=hp*hp
      call hadel(hp,a1,qs,qd,d2)
      dsa=abs(delt-d1-d2)*r
      sc=sin(a1*d2/2.0)
      sm=sin(a2*d1/2.0)
      tc=sqrt2(qds2+4.0*qd*qs*sc*sc)/a1
      tm=sqrt2(qm2+4.0*q1*qm*sm*sm)/a2
      if(dsa.lt.0.1) go to 6
	qe=(qs*qs-hp2)*(qd*qd-hp2)
	if(qe.gt.-0.00001.and.qe.le.0.00001) qe=0.00001
	if(tm.gt.-0.00001.and.tm.le.0.00001) tm=0.00001	
      ddd1=1.0+sqrt2((q1*q1-hp2)*(qm*qm-hp2)/qe)*tc/tm
      db=d1
	if(ddd1.gt.-0.00001.and.ddd1.le.0.00001) ddd1=0.00001		
      d1=d1+(delt-d1-d2)/ddd1
      if(d1.lt.0.000001) d1=db/2.0
      if(d1.gt.delt) d1=delt
    5 continue
      go to 999
    6 t=tc+tm
      th=atan3(sin(a2*d1),q1/qm-cos(a2*d1))
  998 continue
      if(th.lt.0.01) th=3.14159265-th
      dtddel=sin(th)/v1
      dtdr=cos(th)/v1
      return
  999 t=9999.0
      return
      end
      subroutine delth_tc_tm(r1,del,th,tc,tm,dtdr,dtddel,kps)
      common /cvel/al(2,2),be(2,2),r,rm,ntype_velmodel
      va1=al(1,kps)
      va2=al(2,kps)
      be1=be(1,kps)
      be2=be(2,kps)
      v1=va1*(r1/r)**be1
      if(r1.lt.rm) v1=va2*(r1/rm)**be2
      vd=va1*(rm/r)**be1
      q1=r1/v1
      qm=rm/va2
      qd=rm/vd
      qs=r/va1
      a1=1.0-be1
      a2=1.0-be2
      delt=del/r
      call hadel(qm,a1,qs,qd,crdel)
      if(r1.lt.rm) go to 1
      qs1=q1-qs
      qsd=qs-qd
      s=sin(a1*delt/2.0)
      tg=sqrt2(qs1*qs1+4.0*qs*q1*s*s)/a1
      call hadel(qm,a1,q1,qd,d1)
      d2=delt-d1-crdel
      if(d2.lt.0.0) go to 2
ccc	travel time of Pn	***********************************
      do 3 k=1,30
      dth=d2/2.0*a2
      hp=cos(dth)*qm
      call hadel(hp,a1,q1,qd,d11)
      call hadel(hp,a1,qs,qd,d12)
      dsa=abs(delt-d2-d11-d12)*r
      d2=delt-d11-d12
      if(dsa.lt.0.1) go to 4
    3 continue
      go to 999
    4 continue
      s1=sin(d11*a1/2.0)
      s2=sin(d12*a1/2.0)
      tq=(sqrt2((q1-qd)*(q1-qd)+4.0*q1*qd*s1*s1)+sqrt2(qsd*qsd+
     1 4.0*qs*qd*s2*s2))/a1+2.0*qm*sin(d2*a2/2.0)/a2
      if(tq.gt.tg) go to 2
    7 t=tq
      tm=2.0*qm*sin(d2*a2/2.0)/a2
       tc=t-tm     
      s=hp/q1
      c=sqrt2(1.0-s*s)
      th=atan3(s,c)
      go to 998
    2 th=atan3(sin(a1*delt),q1/qs-cos(a1*delt))
      sg=del/(abs(r-r1)+0.01)
      if(sg.lt.0.01.and.th.lt.1.5) th=3.14159265
      sg=sin(th)-qd/q1
      if(th.lt.1.57.and.sg.lt.0.0) then
	if(d2.gt.0.0) go to 7		!r1>radius
	endif
      delmax=3.14/a1
      if(delt.gt.delmax) go to 7
      t=tg
      tc=t
      tm=0.0
      go to 998
    1 continue
      qds2=(qd-qs)*(qd-qs)
      qm2=(q1-qm)*(q1-qm)
      d1=delt-crdel
      if(d1.gt.0.0) go to 11
	if(r.eq.r1) go to 11
      d1=delt*(rm-r1)/(r-r1)
      dh=rm-r1
      if(dh.gt.0.2) go to 11
      d2=delt
      thd=atan3(sin(a1*delt),qd/qs-cos(a1*delt))
      s=sin(d2*a1/2.0)
      t=sqrt2((qs-qd)*(qs-qd)+4.0*qs*qd*s*s)/a1
      tc=t
      s=v1/vd*sin(thd)*rm/r1
      c=-sqrt2(abs(1.0-s*s))
      th=atan3(s,c)
      t=t-dh*c/v1
      tm=t-tc
      go to 998
   11 continue
      do 5 k=1,30
      th=atan3(sin(a2*d1),q1/qm-cos(a2*d1))
           hp=sin(th)*q1
      hp2=hp*hp
      call hadel(hp,a1,qs,qd,d2)
      dsa=abs(delt-d1-d2)*r
      sc=sin(a1*d2/2.0)
      sm=sin(a2*d1/2.0)
      tc=sqrt2(qds2+4.0*qd*qs*sc*sc)/a1
      tm=sqrt2(qm2+4.0*q1*qm*sm*sm)/a2
      if(dsa.lt.0.1) go to 6
	qe=(qs*qs-hp2)*(qd*qd-hp2)
	if(qe.ne.0.0.and.tm.ne.0.0) then
      ddd1=1.0+sqrt2((q1*q1-hp2)*(qm*qm-hp2)/qe)*tc/tm
	endif
      db=d1
      if(ddd1.ne.0.0) d1=d1+(delt-d1-d2)/ddd1
      if(d1.lt.0.000001) d1=db/2.0
      if(d1.gt.delt) d1=delt
    5 continue
      go to 999
    6 t=tc+tm
      th=atan3(sin(a2*d1),q1/qm-cos(a2*d1))
  998 continue
      if(th.lt.0.01) th=3.14159265-th
      dtddel=sin(th)/v1
      dtdr=cos(th)/v1
      return
  999 t=9999.0
      return
      end

	subroutine delth4(r1,del,t) !!!add 03/10/30 by kamimura
	dimension trav(71,201)
	save trav,ibgn
	data ibgn /999/
	if(ibgn.eq.999) then
	 ibgn=0
         open(19,file='./delth2.out',access='direct',recl=300)
	 do idel=1,201
         read(19,rec=idel,iostat=iostat) deldel,
     1   (trav(idep,idel),idep=1,71)
	 if(iostat.ne.0) go to 10
	 enddo
	 close(19)
	 go to 88
 10	 continue
	 write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
	 write(*,*) 'No delth2.out file'
	 stop
	endif
 88	continue
         radius=6371.05
	dep=radius-r1
	idep=ifix(dep/10.0)+1
        idel=ifix(del/10.0)+1
	if(idel.gt.200.or.idep.gt.70) then
	   call delth2(r1,del,th,t,dtdr,dtddel,1)
	   return
	   endif
        tra1=trav(idep,idel)
        tra2=trav(idep,idel+1)
        tra3=trav(idep+1,idel)
        tra4=trav(idep+1,idel+1)
        t1=tra1+((del/10.0-idel+1)*(tra2-tra1))
        t2=tra3+((del/10.0-idel+1)*(tra4-tra3))
        t=t1+((dep/10.0-idep+1)*(t2-t1))
	return
	end

      subroutine delth_routine(r1,del,th,t,dtdr,dtddel,kps)
      common /cvel/al(2,2),be(2,2),r,rm,ntype_velmodel
	call set_vel_routine
      call delth2(r1,del,th,t,dtdr,dtddel,kps)
 	return
	end
	subroutine set_vel_routine
	common /cvel/al(2,2),be(2,2),r,rm,ntype_velmodel
	dimension aq(2,2),bq(2,2)
	data aq(1,1),aq(2,1),aq(1,2),aq(2,2) /5.4,7.75,3.15,4.35/
	data bq(1,1),bq(2,1),bq(1,2),bq(2,2) /-59.0,-2.3,-59.0,-2.3/
	data rq,rmq /6371.05, 6340.05/
	do k=1,2
	do kps=1,2
	al(k,kps)=aq(k,kps)
	be(k,kps)=bq(k,kps)
	enddo
	enddo
	r=rq 
	rm=rmq
	return
	end
	subroutine liner(m,nn)
      common /cliner/cmat(4,9)
      do 10 k=1,m
      w=cmat(k,k)
      if(w.eq.0.0) w=0.001
      kk=k+1
      do 20 j=kk,nn
      cmat(k,j)=cmat(k,j)/w
   20 continue
      do 30 i=1,m
      if(i.eq.k) go to 30
      w=cmat(i,k)
      do 40 j=kk,nn
      cmat(i,j)=cmat(i,j)-w*cmat(k,j)
   40 continue
   30 continue
   10 continue
      return
      end

      subroutine conaz(re,pe,rs,ps,del,az)
c     calculation of delta, azimuth
      pideg=3.14159/180.0
      x1=sin(ps)*sin(re-rs)
      c=sin(pe)*sin(ps)*cos(re-rs)+cos(pe)*cos(ps)
      y1=cos(pe)*sin(ps)*cos(re-rs)-sin(pe)*cos(ps)
      s=sqrt2(x1*x1+y1*y1)
cc	if(s.lt.0.00001) go to 1
      del=atan3(s,c)*6371.05
      az=270.0-atan3(y1,x1)/pideg
      azb=az-360.0
      if(az.gt.360.0) az=azb
   
	if(del.lt.0.3) then
	   s=cos(ps)
	   dx=(rs-re)*6371.05
	   dy=(ps-pe)*6371.05
	   del=sqrt(dx*dx+dy*dy)
	   endif
	return
      end

      subroutine conazd(ram,phi,rams,phis,del,az)
c     calculation of delta, azimuth
ccc   phi,ram degree
      pideg=3.14159/180.0
	if(phi.lt.10.0.and.rams.lt.10.0) go to 5
	pe=(90.0-phi)*pideg
	re=ram*pideg
	ps=(90.0-phis)*pideg
	rs=rams*pideg
	go to 6
    5   continue
	pe=phi
	re=ram
	ps=phis
	rs=rams
    6   continue
      x1=sin(ps)*sin(re-rs)
      c=sin(pe)*sin(ps)*cos(re-rs)+cos(pe)*cos(ps)
      y1=cos(pe)*sin(ps)*cos(re-rs)-sin(pe)*cos(ps)
      s=sqrt2(x1*x1+y1*y1)
cc	if(s.lt.0.00001) go to 1
      del=atan3(s,c)*6371.05
      az=270.0-atan3(y1,x1)/pideg
      azb=az-360.0
      if(az.gt.360.0) az=azb
	if(del.lt.0.5) call condel(ram,phi,rams,phis,del)
	return
      end

	subroutine condel(ram,phi,rams,phis,del)
	parameter(scly=111.195)
	parameter(pid2=3.141592/180.0/2.0)
      dx=(rams-ram)*cos((phi+phis)*pid2)*scly
      dy=(phis-phi)*scly
       del=sqrt(dx*dx+dy*dy)
	return
	end




      subroutine hadel(hp,alf,q1,q2,del)
ccc	subroutine to calculate distance of two points.
ccc	alf=1.0-be(k,kps),  hp=(R-H0)/V0
      s1=hp/q1
      s2=hp/q2
      c1=sqrt2(1.0-s1*s1)
      c2=sqrt2(1.0-s2*s2)
      del=abs(atan3(s1,c1)-atan3(s2,c2))/alf
      return
      end
      subroutine azin(fx,fy,fz,faz,fin)
      fq=sqrt2(fx*fx+fy*fy)
      fin=atan3(fq,fz)
      if(fin.lt.0.0001) go to 1
      faz=atan3(fy,fx)
      if(faz.gt.0.0) return
      faz=faz+2.0*3.14159265
      return
    1 faz=0.0
      return
      end
      subroutine azind(fx,fy,fz,faz,fin)
      fq=sqrt2(fx*fx+fy*fy)
      fin=atan3(fq,fz)
      if(fin.lt.0.0001) go to 1
      faz=atan3(fy,fx)
      if(faz.gt.0.0) go to 2
      faz=faz+2.0*3.14159265
      go to 2
    1 faz=0.0
    2 continue
	faz=faz/3.14159265*180.0
	fin=fin/3.14159265*180.0
      return
      end
      subroutine prout(p,r,dth,daz,xg,yg,zg)
      azm=-3.14159265/2.0-daz
      x=sin(dth)*cos(azm)
      y=sin(dth)*sin(azm)
      z=cos(dth)
      sr=sin(r)
      cr=cos(r)
      sp=sin(p)
      cp=cos(p)
      xg=sr*x+cr*cp*y+cr*sp*z
      yg=-cr*x+sr*cp*y+sr*sp*z
      zg=-sp*y+cp*z
      return
      end
      subroutine azin_prout(p,r,dth,daz,faz,fin)
ccc	inputp,r,dth,daz; radian, output;faz,fin radian
c   routine to calculate the colatitude (fin) and longitude (faz) of
c   a point b, when point a with colatitude p and longitude r, and the
c   azimuth daz and anglar distance dth from a to b are given.
      eps  = 1.0e-6
      azm  = -1.570796-daz
      sdth = sin(dth)
      x    = sdth*cos(azm)
      y    = sdth*sin(azm)
      z    = cos(dth)
      sr   = sin(r)
      cr   = cos(r)
      sp   = sin(p)
      cp   = cos(p)
      xg   = sr*x+cr*cp*y+cr*sp*z
      yg   = -cr*x+sr*cp*y+sr*sp*z
      zg   = -sp*y+cp*z
      qg   = sqrt(xg*xg+yg*yg)
      fin  = atan2(qg,zg)
      faz  = 0.0
      if(qg.gt.eps)   faz = atan2(yg,xg)
      if(faz.lt.0.0)  faz = faz+6.283185
      return
      end
        subroutine chara_zero(text,texta)
        character*80 text
        character*1 texta(80)
        read(text,'(80a1)') texta
        k=0
        do l=1,80
        if(texta(l).ne.' ') k=l
        enddo
        if(k.eq.0) return
        do l=1,k
        if(texta(l).eq.' ') texta(l)='0'
        enddo
        return
        end
	subroutine set_file_name(dkfile,file)
	character*80 dkfile,file
	character*1 cha1(80),cha2(80)
	read(dkfile,'(80a1)') cha1
	n=0
	do k=1,80
	if(cha1(k).ne.' ') then
	n=n+1
	cha2(n)=cha1(k)
	endif
	enddo
	write(file,'(80a1)') (cha2(l),l=1,n)
	return
	end
	subroutine	inside (phi0, ram0, phcnr, rmcnr, n, index2)

!	**********************************************************
!	**	inside		88/10/31	ver.2.0		**
!	**			89/02/12	ver.2.1		**
!	**				debugged		**
!	**							**
!	**	 Subroutine for checking whether a earthquake 	**
!	**	is located in a region given			**
!	**							**
!	**	 This version can be used even if the region	**
!	**	has comlex shape so long as the boundary lines	**
!	**	do not cross each other.			**
!	**	 Algorithm of this routine is originally	**
!	**	given by Dr. Y. Kuwahara.			**
!	**							**
!	**		88/10/31	T. Matsuzawa		**
!	**********************************************************

!	------ Input --------
	REAL	PHI0,		! Latitude of the epicenter
     1 RAM0,		! Longitude of the epicenter
     2 PHCNR(n),	! Latitudes of the corners of the region
     3 RMCNR(n)	! Longitude of the corners of the region
	INTEGER	N		! Total number of corners

!	------ Output -------
	INTEGER	index2		! = 1 : epicenter is inside of the region
				! = 0 : epicenter is outside of the region

!	------ n > 2 ? ------
	if (n.lt.2)	return

	if (n.eq.2)	then	! iF n = 2, phcnr AND rmcnr ARE INTERPRETED
	    index2	= 1	! AS THE LIMITS OF LAT. AND LON. OF THE REGION
	    if (phcnr(1).gt.phcnr(2)) then
		phil 	= phcnr(2)
		phih	= phcnr(1)
	    else
		phil	= phcnr(1)
		phih	= phcnr(2)
	    endif
	    if (rmcnr(1).gt.rmcnr(2)) then
		raml	= rmcnr(2)
		ramh	= rmcnr(1)
	    else
		raml	= rmcnr(1)
		ramh	= rmcnr(2)
	    endif
	    if (phi0.gt.phih .or. phi0.lt.phil)	index2	= 0
	    if (ram0.gt.ramh .or. ram0.lt.raml)	index2	= 0 

	else
!	----- n >= 3 ---------------------
!	   A brief explanation of the algorithm used here is as follows:
!	If the epicenter is located inside of the region, longitudal line
!	of the epicenter must cross the boundary of the region odd times
!	north of the epicenter.  Thus all you have to do is count the crossing
!	point north of the epicenter.

	    PHI00	= PHI0
	    RAM00	= RAM0
	    ICROSS	= 0	! Counter 

!	--- If the longitude of the corner of the boundary coinsides with  
!	  that of the epicenter, shift the epicenter a little --------

	    do i = 1,n
		if (rmcnr(i).eq.ram00) ram00 = ram00 + 0.0001
	    enddo

	    do 10 i = 1, n
		j	= i + 1
		if (j.gt.n)	then
		    phi2	= phcnr(1)
		    ram2	= rmcnr(1)
		else
		    phi2	= phcnr(j)
		    ram2	= rmcnr(j)
		endif
		phi1	= phcnr(i)
		ram1	= rmcnr(i)
		if (phi1.eq.phi0 .and. ram1.eq.ram0)	go to 2000

		if (phi1.lt.phi0 .and. phi2.lt.phi0)	go to 10

		if (ram2.eq.ram0  .and. ram1.eq.ram0) then
		    if (phi1.lt.phi0 .or. phi2.lt.phi0)	go to 2000
		endif

		if (ram1.lt.ram00 .and. ram2.lt.ram00)	go to 10
		if (ram1.gt.ram00 .and. ram2.gt.ram00)	go to 10
		if (phi1.eq.phi00 .and. phi2.eq.phi00)	go to 2000
		if (phi1.ge.phi00 .and. phi2.ge.phi00)	go to 1500

		phi00x	= (ram00-ram1)*(phi2-phi1)/(ram2-ram1) + phi1
		if (phi00x.lt.phi00) go to 10
1500		continue
		icross	= icross + 1
10	    continue
	    index2	= mod(icross,2)
	endif
	return

2000	index2	= 1
	return
	end
	subroutine make_file_name(iy1,im1,id1,ih1,imin1,cha_im,cha_imin)
	character*8 cha_im
	character*14 cha_imin
	character*80 text
	character*1 texta(80)
ccc	open event file  *****************************************
	write(text,'(5i2)') iy1,im1,id1,ih1,imin1
	call chara_zero(text,texta)
  120   format(2x,'texta=',40a1)
	write(cha_im,100) (texta(l),l=1,4)
	write(cha_imin,110) (texta(l),l=1,10)
  100   format(4a1,'.dat')
  110   format(10a1,'.dat')
 	return
	end
       subroutine filter(filin,filout,ndata)
ccc	band pass filter by saito
ccc	*********************************************************
ccc	input							*
cc	filin; an array of input data				*
cc	ndata; number of input data				*
cc	filter_device; defined when you call sub. read_device 	*
cc		       which  read device.dat			*
ccc	*********************************************************
cc
cc	*********************************************************
cc	output							*
cc	filout; an array for out put of filtered data		*
cc	fl,fh,fs,samp are defined in this sub. at first calling	*
cc	*********************************************************
	common /c_filter/ flw,fhw,fsw
	common /cfc/nmem,samp,nstflm,nad,nmemt,nstobs,nch80
	DIMENSION HFIL(50)
	dimension FILOUT(ndata),FILIN(ndata)
	data istate/99/
	if(istate.ne.99) go to 1
	call read_filter_parm
	ap=0.5
	as=5.0
	istate=0
    3   write(*,*) ' fl,fh,fs,ap,as=',flw,fhw,fsw,ap,as
	samp1=samp
	write(*,*)  'samp=',samp
	if(samp.eq.0.0) then
	write(*,*)  'samp is not defined'
	stop
	endif
       dt=1./samp
       fl=flw*dt 
       fh=fhw*dt
       fs=fsw*dt
	if(flw.le.0.0) go to 1
       call butpas(hfil,mfil,gn,nfil,fl,fh,fs,ap,as)
C	write(*,*)  hfil
	WRITE(*,*) ' GN=',gn
	istate=0
    1  continue
	if(samp.ne.samp1) go to 3
	if(fl.le.0.0) then
		if(flw.ne.-50.0) go to 5
		nd=ifix(samp/49.0)
		do l=1,ndata
		l1=l-nd
		if(l1.le.1) l1=1
		filout(l)=filin(l)-filin(l1)
		enddo
		go to 6
    5   do l=1,ndata
	filout(l)=filin(l)
	enddo
    6   continue

	else
      call tandem(filin,filout,ndata,hfil,mfil,1)
	do l=1,ndata
	filout(l)=filout(l)*gn
	enddo
	endif

      return
      end
       subroutine filter_oppsite(filin,filout,ndata)
ccc	band pass filter by saito from opposite direction.
ccc	*********************************************************
ccc	input							*
cc	filin; an array of input data				*
cc	ndata; number of input data				*
cc	filter_device; defined when you call sub. read_device 	*
cc		       which  read device.dat			*
ccc	*********************************************************
cc
cc	*********************************************************
cc	output							*
cc	filout; an array for out put of filtered data		*
cc	fl,fh,fs,samp are defined in this sub. at first calling	*
cc	*********************************************************
	common /c_filter/ flw,fhw,fsw
	common /cfc/nmem,samp,nstflm,nad,nmemt,nstobs,nch80
	DIMENSION HFIL(50)
	dimension FILOUT(ndata),FILIN(ndata)
	data istate/99/
	if(istate.ne.99) go to 1
	call read_filter_parm
	ap=0.5
	as=5.0
	istate=0
    3   write(*,*) ' fl,fh,fs,ap,as=',flw,fhw,fsw,ap,as
	samp1=samp
	write(*,*)  'samp=',samp
	if(samp.eq.0.0) then
	write(*,*)  'samp is not defined'
	stop
	endif
       dt=1./samp
       fl=flw*dt 
       fh=fhw*dt
       fs=fsw*dt
	if(flw.le.0.0) go to 1
       call butpas(hfil,mfil,gn,nfil,fl,fh,fs,ap,as)
C	write(*,*)  hfil
	WRITE(*,*) ' GN=',gn
	istate=0
    1  continue
	if(samp.ne.samp1) go to 3
	if(fl.le.0.0) then
		if(flw.ne.-50.0) go to 5
		nd=ifix(samp/49.0)
		do l=1,ndata
		l1=l-nd
		if(l1.le.1) l1=1
		filout(l)=filin(l)-filin(l1)
		enddo
		go to 6
    5   do l=1,ndata
	filout(l)=filin(l)
	enddo
    6   continue

	else
      call tandem(filin,filout,ndata,hfil,mfil,-1)
	do l=1,ndata
	filout(l)=filout(l)*gn
	enddo
	endif

      return
      end

	subroutine read_filter_parm
	common /c_filter/ flw,fhw,fsw
	character*80 file
	file='filter.dat'
	open(7,file='filter.dat',status='old')
CCC    FILTER
        READ(7,100)  flw
	read(7,100) fhw
	read(7,100) fsw
	write(*,*)  'read filter parm',flw,fhw,fsw
  100   format(f10.1)  
	close(7)
	return
	end
        subroutine filter2(samp,filin,filout,ndata)
	common /c_filter/ flw,fhw,fsw
	DIMENSION HFIL(50)
	dimension FILOUT(ndata),FILIN(ndata)
	data istate/99/
	if(istate.ne.99) go to 1
	call read_filter_parm
	ap=0.5
	as=5.0
    3   write(*,*) ' fl,fh,fs,ap,as=',flw,fhw,fsw,ap,as
	samp1=samp
	write(*,*)  'samp=',samp
	if(samp.eq.0.0) then
	write(*,*)  'samp is not defined'
	stop
	endif
       dt=1./samp
       fl=flw*dt 
       fh=fhw*dt
       fs=fsw*dt
	if(flw.le.0.0) go to 1
       call butpas(hfil,mfil,gn,nfil,fl,fh,fs,ap,as)
C	write(*,*)  hfil
	WRITE(*,*) ' GN=',gn
	istate=0
    1  continue
	if(samp.ne.samp1) go to 3
	if(fl.le.0.0) then
		if(flw.ne.-50.0) go to 5
		nd=ifix(samp/49.0)
		do l=1,ndata
		l1=l-nd
		if(l1.le.1) l1=1
		filout(l)=filin(l)-filin(l1)
		enddo
		go to 6
    5   do l=1,ndata
	filout(l)=filin(l)
	enddo
    6   continue

	else
      call tandem(filin,filout,ndata,hfil,mfil,1)
	do l=1,ndata
	filout(l)=filout(l)*gn
	enddo
	endif

      return
      end
       subroutine filter_band(samp,filin,filout,ndata)
	common /c_filter/ flw,fhw,fsw
ccc	band pass filter by saito
ccc	*********************************************************
ccc	input							*
cc	filin; an array of input data				*
cc	ndata; number of input data				*
cc	flw,fhw,fsw; filtering range in Hz.
ccc	samp; sampling frequency in Hz.
cc	*********************************************************
cc	output							*
cc	filout; an array for out put of filtered data		*
	DIMENSION HFIL(50)
	dimension FILOUT(ndata),FILIN(ndata)
	ap=0.5
	as=5.0
       dt=1./samp
       fl=flw/samp
       fh=fhw/samp
       fs=fsw/samp
       call butpas(hfil,mfil,gn,nfil,fl,fh,fs,ap,as)
      call tandem(filin,filout,ndata,hfil,mfil,1)
	do l=1,ndata
	filout(l)=filout(l)*gn
	enddo
      return
      end


c	subroutine programme - tandem -
c	    recursive filtering in series 
	subroutine tandem(x,y,n,h,m,nml) 
	dimension h(50) 
	dimension x(n),y(n) 
c		arguments 
c	x	: input time series 
c	y	: output time series (may be equivalent to x) 
c	n	: length of x & y 
c	h	: coefficients of filter
c	m 	: order of filter
c	nml	: >0 ; for normal  direction filtering
c		  <0 ;     reverse direction filtering
c	subroutine required  :  recfil
	if (n .le. 0 .or. m .le. 0) goto 2
	if(m.gt.12) m=12
c***  	 1-stc all
	call recfil(x,y,n,h,nml)
	if(m .le. 1) return
c***	 2-nd and after
	do 1 i=2,m
	  call recfil (y,y,n,h(i*4-3),nml)
1	continue
	return
c*** 	error 
2	write(6,3) n,m 
3	format('(tandem) invalid input',3x,'n = ',2i5)
	return 
	end 
c 
c	subroutine programme - recfil - 
c	  recursive filtering  : f(z) = (1+a*z+aa*z**2)/(1*b*z+bb*z**2)
	subroutine recfil(x,y,n,h,nml) 
	dimension x(n),y(n) 
	dimension h(50) 
	K=14
	if(n .le. 0) goto 4 
	if(nml .ge. 0) goto 1 
c***	reverse filtering 
	j  =  n 
	jd = -1
	  goto 2 
c***   	normal filtering 
1	j  = 1
	jd = 1 
c 
2	a  = h(1) 
	aa = h(2) 
	b  = h(3) 
	bb = h(4) 
	u1 = 0. 
	u2 = 0. 
	v1 = 0.
	v2 = 0. 
c*** 	filtering 
	    do 3 i = 1,n 
		u3 = u2 
		u2 = u1 
		u1 = x(j) 
		v3 = v2 
		v2 = v1
		v1 = u1 + a*u2 + aa*u3 - b*v2 - bb*v3
		y(j)= v1
		j  = j + jd
3	    continue
c***	exit 
	return 
c***	error 
4	write(k,5) n 
5	format(//1x,5('?'),3x,'(recfil)',3x,'invalid input',3x,'n =',2i5) 
	return 
	end
c	 
c	subroutine programme - butpas - 
c	   butterworth band pass filter coefficients 
	subroutine butpas(hfil,mfil,gn,nfil,fl,fh,fs,ap,as) 
c 
	complex r(2),oj,cq 
	dimension hfil(50)
	data pi/3.141593/,hp/1.570796/ 
ccc
ccc	write(*,*)  fl,fh,fs,ap,as
ccc
c 
c	   arguments 
c	hfil	: filterc oefficients
c	mfil	: order of filter 
c	gn	: gain factor 
c	nfil	: order of butterworth function 
c	fl	: low frequency c ut-off  ( non-dimensional ) 
c	fh	: high frequencyc ut-off 
c	fs	: stop band frequency 
c	ap	: max. attenuation in pass band 
c	as	: min. attenuation in stop band 
c 
c	tan(x) = sin(x)/cos(x) 
	wl = amin1(abs(fl),abs(fh))*pi 
	wh = amax1(abs(fl),abs(fh))*pi 
	ws = abs(fs)*pi 
	if(wl .eq. 0.  .or. wl .eq. wh) go to 100
	if(wh .ge. hp .or. ws .eq. 0.0) go to 100 
        if(ws .ge. hp .or. (ws - wh)*(ws - wh) .le. 0.) goto 100 
c*** 	determine n &c  
	clh = 1./(cos(wl)*cos(wh))
	op  = sin(wh - wl)*clh 
	ww  = tan(wl)*tan(wh) 
	ts  = tan(ws)
	os  = abs(ts - ww/ts) 
	pa  = amin1(abs(ap),abs(as)) 
	sa  = amax1(abs(ap),abs(as)) 
	if (pa .eq. 0.) pa = 0.5 
	if (sa .eq. 0.) sa = 5. 
	nfil  = max0(2,ifix(abs(alog(pa/sa)/alog(op/os)) + 0.5)) 
	cc = exp(alog(pa*sa)/float(nfil))/(op*os) 
	c  = sqrt2(cc) 
	ww = ww*cc 
c 
		dp = hp/float(nfil) 
		k  = nfil/2 
		mfil  = k*2 
		l  = 0 
		g  = 1. 
		fj = 1. 
c
	do 2 j = 1,k
	oj = cmplx(cos(dp*fj),sin(dp*fj))*0.5
	fj = fj+2.
	cq =c sqrt(oj**2 + ww)
	r(1) = oj + cq
	r(2) = oj - cq
 	g = g *cc
		do 1 i = 1,2
		re =  real(r(i))**2 
		ri = aimag(r(i)) 
		a  = 1./((c+ri)**2 + re) 
		g  = g*a 
		hfil(l+1) = 0. 
		hfil(l+2) = -1.  
		hfil(l+3) = 2.*((ri - c)*(ri + c) + re)*a 
		hfil(l+4) = ((ri - c)**2 + re)*a 
		l = l + 4 
1		continue 
2	continue 
c***	exit
	gn = g 
	if(nfil .eq. mfil) return
c*** 	for odd nfil
	mfil = mfil+1
	wpc = cc*cos(wh - wl)*clh
	wmc =-cc*cos(wh + wl)*clh 
		a  = 1./(wpc + c) 
		gn = g*c*a 
		hfil(l+1) = 0. 
		hfil(l+2) = -1. 
		hfil(l+3) = 2.*wmc*a 
		hfil(l+4) = (wpc - c)*a 
	return 
c***	error 
100	write(*,101) fl,fh,fs 
101	format(/1x,'(butpas)   invalid input',3e14.5)
	return 
	end 

        subroutine funagc(n200,ndata,nwind,amax,ibuf,buf)
	dimension ibuf(8000),buf(8000),filout(8000),agca(8000)
cc	calculate average **************************************
	zab=0.0
	zav=0.0
	nn=ndata/5+1
	if(nn.gt.ndata) nn=ndata
	do l=1,nn
	zav=zav+float(ibuf(l))
	enddo
	ib=ifix(zav/float(nn))
	do l=1,ndata
	z=float(ibuf(l)-ib)
	buf(l)=z
	enddo
cc	call filtering subroutine *****************************
	call filter(buf,filout,ndata)
cc	calculate agc filter************************************	
	zab=0.0
	nn=nwind
	if(nn.le.0) then
cc	in the case of not using agc filter
	do l=1,ndata
	buf(l)=filout(l)
	enddo
cc	in the case of using agc filter
	else 
	do l=1,nn
	z=abs(filout(l))
	agca(l)=z
	zab=zab+z
	enddo
	if(zab.le.0.001) zab=0.001
	do l=1,nn
	buf(l)=filout(l)/zab
	enddo
	nn1=nn+1
	j=0
	do l=nn1,ndata
	buf(l)=filout(l)/zab
	j=j+1
	if(j.gt.nn) j=1
	aa=abs(filout(l))
	zab=zab+aa-agca(j)
	if(zab.le.0.001) zab=0.01
	agca(j)=aa
	enddo
	endif
cc	find amax  ***********************************
	amax=0.001
	n1=n200
	if(n1.le.0) n1=1
	do l=n1,ndata
	a=abs(buf(l))
	if(a.gt.amax) amax=a
	enddo
      return
      end
	subroutine i2dgtodec(idg,idec)
	logical*1 iadgl(2),iadecl(2)
	integer*2 iadec,iadg,idg,idec
	equivalence (iadgl(1),iadg),(iadecl(1),iadec)
		iadg=idg
		iadecl(1)=iadgl(2)
		iadecl(2)=iadgl(1)
		idec=iadec
	return
	end
	subroutine i4dgtodec(idg,idec)
	logical*1 iadgl(4),iadecl(4)
	integer*4 iadec,iadg,idg,idec
	equivalence (iadgl(1),iadg),(iadecl(1),iadec)
		iadg=idg
		iadecl(1)=iadgl(4)
		iadecl(2)=iadgl(3)
		iadecl(3)=iadgl(2)
		iadecl(4)=iadgl(1)
		idec=iadec
	return
	end
	function	rdgtodec (rdgdata)

	!*********************************************
	!** 	Convert DG-formatted Real Data	    **
	!**		into DEC-formatted Real Data**
	!**			88/04/13	    **
	!**			ver.0.0		    **
	!**			by T.Matsuzawa	    **
	!*********************************************

	byte		bt(4)
	real*4		dgdt
	integer*2	ibt(4)
	equivalence	(dgdt,bt)

	dgdt	= rdgdata

	do i = 1,4
		if (bt(i).ge.0)	then
			ibt(i)	= bt(i)
		else
			ibt(i)	= bt(i) + 256
		endif
	enddo

	fr	= 0.
	do i = 4,2,-1
		fr	= (fr + ibt(i))/256.
	enddo

	iex	= ibt(1)
	if (iex.ge.128)	then
		fr	= - fr
		iex	= iex - 128
	endif
	IF(IEX.GE.95) IEX=95
	ex		= 16.**(iex-64)
	rdgtodec	= ex*fr

	return
	end
        function atan3(s,c)
        sg=abs(s)+abs(c)
        if(sg.lt.0.000001) c=-0.000001
        atan3=atan2(s,c)
        return
        end
	function sqrt2(s)
	sq=abs(s)
	sqrt2=sqrt(sq)
        return
        end

      subroutine     wbzone(gorphi,gorram,ordep,dst,ind)                
c                                                                       
      dimension      ram(6,11),phinet(6),ramnet(6),xphi(3),xram(3),     
     &               xdep(3),x(3),y(3),z(3)
!
!      data ram/54*0.0,148.40,142.90,140.20,138.00,135.85,133.35,
!     *         150.20,146.25,143.00,139.00,136.45,133.70/
!
      data ram/141.95,140.25,139.25,138.00,136.50,135.00,               
     *         142.15,139.95,138.65,137.20,135.45,133.80,               
     *         142.70,139.55,138.00,136.45,134.65,132.75,               
     *         143.50,139.90,137.95,136.35,134.50,132.50,               
     *         143.95,140.50,138.55,136.70,134.70,132.55,               
     *         144.15,140.80,139.00,137.00,134.95,132.65,               
     *         144.30,141.00,139.25,137.25,135.15,132.85,               
     *         144.70,141.10,139.45,137.50,135.35,133.00,               
     *         146.75,141.35,139.60,137.65,135.55,133.20,
     *         148.40,142.90,140.20,138.00,135.85,133.35,
     *         150.20,146.25,143.00,139.00,136.45,133.70/
      data phinet/37.0, 38.0, 39.0, 40.0, 41.0, 42.0/                   
      data ramnet/139.8,140.1,140.8,140.9,140.95,141.2/                 
c                                                                       
      phi0    = 34.0                                                    
      if(ordep.lt.50.0)                                go to 100        
      jlat    = gorphi-phi0+1.0                                         
      if(jlat.le.0.or.jlat.ge.11)                      go to 100        
      hasu    = gorphi-int(gorphi)                                     
      ramtr   = ram(1,jlat)+(ram(1,jlat+1)-ram(1,jlat))*hasu            
      if(gorram.ge.ramtr.and.ordep.le.150.0)           go to 100        
      ramwt   = ram(4,jlat)+(ram(4,jlat+1)-ram(4,jlat))*hasu            
      if(gorram.le.ramwt.and.ordep.le.150.0)           go to 100        
      if(gorphi.ge.42.0.and.ordep.lt.100.0)            go to 100        
      if(gorphi.le.36.5.and.ordep.lt.100.0)            go to 100        
      ramct    = ram(2,jlat)+(ram(2,jlat+1)-ram(2,jlat))*hasu           
      ramct    = (9*ramct+ramtr)/10.0                                   
      if(gorram.gt.ramct.and.ordep.lt. 50.0)           go to 100        
      go to 102                                                         
  100 dst      = 0.0                                                    
      ind      = 0                                                      
      return                                                            
  102 continue                                                          
      do 103 i=1,6                                                      
      ramcr    = ram(i,jlat)+(ram(i,jlat+1)-ram(i,jlat))*hasu           
      if(gorram.ge.ramcr)                              go to 104        
  103 continue                                                          
  104 idep     = i-1                                                    
      if(idep.le.0)            idep = 1                                 
      ramc     = (ram(idep,jlat)+ram(idep+1,jlat+1))/2.0                
      xphi(1)  = gorphi                                           
      xphi(3)  = xphi(1)+1.0                                            
      xram(1)  = ram(idep,jlat)                                         
      xram(3)  = ram(idep+1,jlat+1)                                     
      xdep(1)  = (idep-1)*100                                           
      xdep(3)  = xdep(1)+100.0                                          
      if(gorram.ge.ramc)                                  go to 105     
      xphi(2)  = xphi(1)                                                
      xram(2)  = ram(idep+1,jlat)                                       
      xdep(2)  = xdep(1)+100.0                                          
      go to 106                                                         
  105 xphi(2)  = xphi(1)+1.0                                            
      xram(2)  = ram(idep,jlat+1)                                       
      xdep(2)  = xdep(1)                                                
  106 continue                                                          
      pideg    = 3.141593/180.0                                         
      r        = 6369.5                                                 
      do 107 i=1,3                                                      
      s        = (90.0-xphi(i))*pideg                                   
      t        = xram(i)*pideg                                          
      u        = r-xdep(i)                                              
      x(i)     = u*sin(s)*cos(t)+3000.0                                 
      y(i)     = u*sin(s)*sin(t)                                        
      z(i)     = u*cos(s)
107     continue      
      ap       = y(1)*(z(2)-z(3))+y(2)*(z(3)-z(1))+y(3)*(z(1)-z(2))     
      bp       = x(1)*(z(3)-z(2))+x(2)*(z(1)-z(3))+x(3)*(z(2)-z(1))     
      cp       = x(1)*(y(2)-y(3))+x(2)*(y(3)-y(1))+x(3)*(y(1)-y(2))     
      dp       = x(1)*(z(2)*y(3)-z(3)*y(2))+x(2)*(y(1)*z(3)-y(3)*z(1))  
     &          +x(3)*(y(2)*z(1)-y(1)*z(2))                             
      ep       = sqrt(ap*ap+bp*bp+cp*cp)                                
      if(dp.gt.0.0)            ep = -ep                                 
      plx      = ap/ep                                                  
      pmx      = bp/ep                                                  
      pnx      = cp/ep                                                  
      ppx      =-dp/ep                                                  
      call     dist_plane(gorphi,gorram,ordep,plx,pmx,pnx,ppx,dst)    
      dst      = dst+15.0                                               
      pai2     = 1.570797                                               
      p11      = pai2-gorphi*pideg                                      
      r11      = gorram*pideg                                           
      delmn    = 600.0                                                  
      do 108 i=1,6                                                      
      p22      = pai2-phinet(i)*pideg                                   
      r22      = ramnet(i)*pideg                                        
      call                     conaz(r11,p11,r22,p22,del,az)            
      if(del.lt.delmn)         delmn = del                              
  108 continue                                                          
      haba     = 50.0+(delmn-100.0)/5                                   
      if(delmn.le.100.0)       haba = 50.0                              
      habad    = 35.0+(delmn-100.0)/5                                   
      if(delmn.le.100.0)       habad=35.0                               
      ind      = 5                                                      
      if(ordep.lt. 50.0)       ind = 1                                  
      if(abs(dst).le.haba)     ind = 1                                  
      if(abs(dst).le.habad)    ind = 0                                  
      return                                                            
      end                                                               


	subroutine dist_plane(gorphi,gorram,ordep,plx,pmx,pnx,ppx,pdist)
!
!	Distance between Hypocenter and Projection Plane
!
	deg	= 180.0/3.14159265
	r	= 6369.5
	s	= (90.0-gorphi)/deg
	t	= gorram/deg
	u	= r-ordep
	x	= u*sin(s)*cos(t) + 3000.0
	y	= u*sin(s)*sin(t)
	z	= u*cos(s)
	pdist	= x*plx+y*pmx+z*pnx-ppx
	return
	end
	subroutine wbzone_bosai(phi,ram,dep,dis,ind)		!dis>0 upper
	dimension ramg(4,7), phig(7),ca(11),ansa(11),depg(4)
      dimension cmat(10,11)
	data phig /34.0,34.5,35.0,35.5,36.0,36.5,37.0/
	data depg /50.0,100.0,150.0,200.0/
      data ramg/141.00,139.95,139.65,138.80,           !34.0             
     *         140.80,140.12,139.42,138.68,            !34.5
     *         140.70,139.95,139.17,138.40,            !35.0
     *         141.00,139.88,138.85,137.80,            !35.5
     *         141.00,139.85,138.90,138.00,            !36.0
     *         141.25,140.00,139.00,138.15,            !36.5
     *         141.50,140.10,139.05,138.20 /           !37.0
	data ibgn /99/
	if(ibgn.ne.99) go to 10
	ibgn=0
	do k=1,10
	do m=1,11
	cmat(k,m)=0.0
	enddo
	enddo
	do k=1,7
	do m=1,4
	call fun_ca(ca,phig(k),ramg(m,k),depg(m))
	do kk=1,10
	do mm=1,11
	cmat(kk,mm)=cmat(kk,mm)+ca(kk)*ca(mm)
	enddo
	enddo
	enddo
	enddo

	call liner3(10,11,cmat)
	do k=1,10
	ansa(k)=cmat(k,11)
	enddo
   10   continue
	call fun_ca(ca,phi,ram,dep)
	dis=-dep
	do k=1,10
	dis=dis+ca(k)*ansa(k)
	enddo
	return
	end

      subroutine liner3(m,nn,cmat)
      dimension cmat(m,nn)
      do 10 k=1,m
      w=cmat(k,k)
      if(w.eq.0.0) w=0.00001
      kk=k+1
      do 20 j=kk,nn
      cmat(k,j)=cmat(k,j)/w
 20   continue
      do 30 i=1,m
      if(i.eq.k) go to 30
      w=cmat(i,k)
      do 40 j=kk,nn
      cmat(i,j)=cmat(i,j)-w*cmat(k,j)
   40 continue
   30 continue
   10 continue
      return
      end
	subroutine fun_ca(ca,phi,ram,dep)
	dimension ca(11)
	pq=phi-35.0
	rq=ram-140.0
	ca(1)=1.0
	ca(2)=pq
	ca(3)=rq
	ca(4)=pq*pq
	ca(5)=pq*rq
	ca(6)=rq*rq
	ca(7)=pq*pq*pq
	ca(8)=pq*pq*rq
	ca(9)=pq*rq*rq
	ca(10)=rq*rq*rq
	ca(11)=dep
	return
	end





	subroutine get_minimum_order(dela,nnn,mmm,num_ansa)
	dimension dela(nnn),num_ansa(mmm)
	dimension qqqa(20000)
	do k=1,nnn
	qqqa(k)=dela(k)
	enddo
	do m=1,mmm
	delmin=1.0*100000.0*100000.0*100000.0
	do k=1,nnn
	if(delmin.ge.qqqa(k)) then
	delmin=qqqa(k)
	num_ansa(m)=k
	endif
	enddo
	mm=num_ansa(m)
	qqqa(mm)=0.9*100000.0*100000.0*100000.0
	enddo
	return
	end
	subroutine get_maximum_order(dela,nnn,mmm,num_ansa)
	dimension dela(nnn),num_ansa(mmm)
	dimension qqqa(20000)
	do k=1,nnn
	qqqa(k)=dela(k)
	enddo
	do m=1,mmm
	delmax=-1.0*100000.0*100000.0*100000.0
	do k=1,nnn
	if(delmax.le.qqqa(k)) then
	delmax=qqqa(k)
	num_ansa(m)=k
	endif
	enddo
	mm=num_ansa(m)
	qqqa(mm)=-0.9*100000.0*100000.0*100000.0
	enddo
	return
	end

	subroutine seismo_cor(f0,h,w,amp)
	complex  amp,c
        c=(0.0,1.0)
ccc	input f0,h,w; frequency of seismometer, h constant, angular frequency.
ccc	output amp; a complex value of amplitude for output of seismometer. 
	cf=2.0*3.14159265*f0
	cf2=cf*cf
ccc	d2a/dt2 + 2hn da/dt + n*n*a = - d2x/dt2
	amp=w*w/(w*w-cf2-2.0*h*cf*w*c)
	return
	end
	subroutine amp_p_inc(thsw,vp,vs,upx,upz)
ccc	calculation of amplitude at surface by P-wave incidence.
	pideg=3.14159265/180.0
	ths=thsw
	if(ths.gt.90.0) ths=180.0-ths
	th=ths*pideg
	sp=sin(th)
	cp=cos(th)
	p=sp/vp
	q=cos(th)/vp
	ss=p*vs
	cs=sqrt(1.0-ss*ss)
	r=cs/vs
	qm=1.0-2.0*vs*vs*p*p
	vs2=vs*vs
	vs4=vs2*vs2
	p2=p*p
	qmat=4.0*vs4*p2*q*r+qm*qm
	B=(4.0*vs4*p2*q*r-qm*qm)/qmat
	C=-4.0*vs2*p*q*qm/qmat
	upx=(p*(1.0+B)-C*r)*vp		!upx_in=sin(th)
	upz=-((B-1.0)*q+C*p)*vp		!upz_in=cos(th) upward +
ccc	write(*,*)  'B,C=', B,C
ccc	write(*,*) 'p,cs,vs,vp=',p,cs,vs,vp
ccc	write(*,'(2x,f6.1,4f7.2)') ths,B,C,upx,upz
	return
	end
	subroutine amp_s_inc(thsw,vp,vs,upx,upxi,upz,upzi)
ccc	calculation of amplitude at surface by S-wave incidence.
	complex*8 q,qmat,A,B,C,a1,a2,a3,a4,b1,b2,upxim,upzim
	pideg=3.14159265/180.0
	ths=thsw
	if(ths.gt.90.0) ths=180.0-ths
	th=ths*pideg
	sp=sin(th)
	cp=cos(th)
	p=sp/vs
	r=cos(th)/vs
	ss=p*vp
	sa=1.0-ss*ss
	if(sa.gt.0.0) then
	q=sqrt(sa)/vp
	else 
	q1=sqrt(-sa)
	q=cmplx(0.0,q1)/vp
	endif
cc	write(*,*)  'sa=',sa
	vs2=vs*vs
	vs4=vs2*vs2
	g=1.0-2.0*vs2*p*p
	a1=(2.0*p)*q
	a2=p*p-r*r
	a3=1.0-2.0*vs2*p*p
	a4=2.0*vs2*p*r
	b1=-a2
	b2=2.0*vs2*p*r
	qmat=a1*a4-a2*a3
	B=(b1*a4-a2*b2)/qmat	!complex
	C=(a1*b2-a3*b1)/qmat
	A=1.0
	upxim=(p*B+(A-C)*r)*vs		!upx_in=sin(th)
	upzim=-(B*q+(A+C)*p)*vs		!upz_in=cos(th) upward +
	upx=real(upxim)
	upxi=imag(upxim)
	upz=real(upzim)
	upzi=imag(upzim)
cc	write(*,*)  'txy,txy1,txy2,C=',txy,txy1,txy2,C
ccc	write(*,100) ths,txz,tzz,upx,upxi,upz,upzi
ccc  100   format(2x,f7.3,4f7.4,5x,2(2f6.2,2x))
	return
	end
	subroutine day_year(idd,iy,im,id)
	if(im.le.2) then
	idd=(im-1)*31+id
	go to 1
	endif
	iuruu=iy-(iy/4)*4
	id3=31+28
	if(iuruu.eq.0) id3=id3+1
	idd=id3+(im-3)*31+id		!4,6,9,11
	if(im.ge.5) idd=idd-1
	if(im.ge.7) idd=idd-1
	if(im.ge.10) idd=idd-1
	if(im.ge.12) idd=idd-1
    1   continue
	idd=idd+366*iy
	return
	end

             
   	subroutine iyimid_idd(idd,iy,im,id)
ccc	input iy,im,id,   output idd
	if(im.le.2) then
	idd=(im-1)*31+id
	go to 1
	endif
	iuruu=iy-(iy/4)*4
	id3=31+28
	if(iuruu.eq.0) id3=id3+1
	idd=id3+(im-3)*31+id		!4,6,9,11
	if(im.ge.5) idd=idd-1
	if(im.ge.7) idd=idd-1
	if(im.ge.10) idd=idd-1
	if(im.ge.12) idd=idd-1
    1   continue
        iyq=iy
cc        if(iyq.gt.1900) iyq=iyq-2000
        n_uru=(iyq-1)/4
       
	idd=idd+365*iy+n_uru

	return
	end

	subroutine idd_iyimid(idd,iy,im,id)
ccc	 input; idd   output iy,im,id 
	iy=idd/365
        iyq=iy
cc        if(iyq.gt.1900) iyq=iyq-2000
        n_uru=(iyq-1)/4
	id=idd-n_uru
	iuruu=iy-(iy/4)*4
	id=id-iy*365
	im=1
cc        write(*,*) idd,iy,im,id,iuruu,n_uru

	if(id.le.0) then
	   id=32+id-1
	   im=12
	   iy=iy-1
	   return
	   endif
	if(id.le.31) return
        im=2
	id=id-31
	iq=28
	if(iuruu.eq.0) iq=29
	if(id.le.iq) return
	id=id-iq
	im=3
	if(id.le.31) return
	id=id-31	!April
	im=4
	if(id.le.30) return
	id=id-30
	im=5
	if(id.le.31) return
	id=id-31
	im=6
	if(id.le.30) return
	id=id-30
	im=7
	if(id.le.31) return
	id=id-31
	im=8
	if(id.le.31) return
	id=id-31
	im=9
	if(id.le.30) return
	id=id-30
	im=10
	if(id.le.31) return
	id=id-31
	im=11
	if(id.le.30) return
	id=id-30
	im=12
	return
	end
	subroutine del_order(nnn,dela,ich_ord)
ccc	order according to small distances
	dimension dela(nnn),ich_ord(nnn)
	dimension delq(20000)
	do k=1,nnn
	delq(k)=dela(k)
	enddo
	del_minw=10.0**30.0
	do m=1,nnn
	del_min=del_minw
	do k=1,nnn
	if(delq(k).lt.del_min) then
	del_min=delq(k)
	nk=k
	endif
	enddo
	delq(nk)=del_minw*0.999
	ich_ord(m)=nk
	enddo
	return
	end
	subroutine plot_order(nnn,icha,dela,ich_ord)
ccc	input icha,dela; ch number and its distance
ccc	output ich_ord; ch number in the order of small distance.
	dimension dela(nnn),ich_ord(nnn),icha(nnn)
	dimension delq(20000)
	delmax=10.0**20.0
	do k=1,nnn
	delq(k)=dela(k)
	if(delq(k).gt.delmax) delq(k)=delmax
	enddo
	del_minw=10.0**20.0
	do m=1,nnn
	del_min=del_minw
	nk=m
	do k=1,nnn
	if(delq(k).lt.del_min) then
	del_min=delq(k)
	nk=k
	endif
	enddo
	delq(nk)=del_minw*0.999
	ich_ord(m)=icha(nk)
	enddo
	return
	end
  	subroutine small_order(nnn,nord,dela,ich_ord)
ccc	input dela;  distance etc
ccc	output ich_ord; ch number in the order of small distance.
	dimension dela(nnn),ich_ord(nnn),icha(nnn)
	dimension delq(100000)
	delmax=10.0**20.0
	do k=1,nnn
	delq(k)=dela(k)
	if(delq(k).gt.delmax) delq(k)=delmax
	enddo
	del_minw=10.0**20.0
	do m=1,nord
	del_min=del_minw
	nk=m
	do k=1,nnn
	if(delq(k).lt.del_min) then
	del_min=delq(k)
	nk=k
	endif
	enddo
	delq(nk)=del_minw*0.999
	ich_ord(m)=nk
	enddo
	return
	end
      subroutine rev_matrix(a,n,y)
      dimension a(n,n),y(n,n),indx(n)
      np=n
      do i=1,n
         do j=1,n
            y(i,j)=0.0
         enddo
         y(i,i)=1.0
      enddo
c decompose A = LU
      call ludcmp(a,n,np,indx,d)
c inverse matrix A
      do j=1,n
         call lubksb(a,n,np,indx,y(1,j))
c        Note that FORTRAN stores two-dimensional matricses
c        by column, so y(1,j) is the address of the jth 
c        column of y.
      enddo
      return
      end

      SUBROUTINE ludcmp(a,n,np,indx,d)
      INTEGER n,np,indx(n),NMAX
      REAL d,a(np,np),TINY
      PARAMETER (NMAX=500,TINY=1.0e-20)
c Given a matrix a(1:n,1:n), with physical dimension np by np,
c this routine replaces it by the LU decomposition of a rowwise
c permutation of itself. a and n are input.
c a is output, arranged as in equation (2.3.14) above; index(1:n)
c is an output vector that records the row permutation effected by
c partial povoting; dis output as +/- 1 depending on whether
c the number of row interchanges was even or odd, respectively.
c This routine is used in combination with lubksb to solve
c linear equations or ivert a matrix.
 
      INTEGER i,imax,j,k
      REAL aamax,dum,sum,vv(NMAX)
 
      d=1
      do i=1,n
         aamax=0
         do j=1,n
            if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
         enddo
      if (aamax.eq.0) stop
      vv(i)=1./aamax
      enddo
       
      do j=1,n
         do i=1,j-1
           sum=a(i,j)
           do k=1,i-1
              sum=sum-a(i,k)*a(k,j)
           enddo
           a(i,j)=sum
         enddo
         aamax=0
         do i=j,n
           sum=a(i,j)
           do k=1,j-1
              sum=sum-a(i,k)*a(k,j)
           enddo
           a(i,j)=sum
           dum=vv(i)*abs(sum)
           if (dum.ge.aamax) then
               imax=i
               aamax=dum
           endif
         enddo
         if (j.ne.imax) then
             do k=1,n
               dum=a(imax,k)
               a(imax,k)=a(j,k)
               a(j,k)=dum
             enddo
             d=-d
             vv(imax)=vv(j)
         endif
         indx(j)=imax
         if(a(i,j).eq.0.) a(i,j)=TINY
         if(j.ne.n) then
            dum=1./a(j,j)
            do i=j+1,n
               a(i,j)=a(i,j)*dum
            enddo
         endif
      enddo
      return
      END
  
c#############################################################
 
      SUBROUTINE lubksb(a,n,np,indx,b)
      INTEGER n,np,indx(n)
      REAL a(np,np),b(n)
c solve the set of n linear equations Ax=b. Here is a input,
c not at the matrix A but rather as its LU decompsition,
c determined by the routine ludcmp. indx is input as the
c permutation vector returned by ludcmp. b(1:n) is input
c as the right-hand side vectror b, and returns with the
c solution vector x.
c a, n, np and indx are not modified by this routine
c and be left in place for successive calls with different
c right hand sides b. This routine takes into account
c the possibility that b will begin with many zero elements,
c so it is efficient for use in matrix inversion.
 
      INTEGER i,ii,j,ll
      REAL sum
 
      ii=0
      do i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if (ii.ne.0) then
            do j=ii,i-1
               sum=sum-a(i,j)*b(j)
            enddo
        else if (sum.ne.0) then
            ii=i
        endif
        b(i)=sum
      enddo
      do i=n,1,-1
        sum=b(i)
        do j=i+1,N
        sum=sum-a(i,j)*b(j)
        enddo
        b(i)=sum/a(i,i)
      enddo
      return
      END
	subroutine idd_iyimid2(idd,iy,im,id)
ccc	 input; idd,iy  output im,id
	iyear=iy+2000
	id=idd
	iuruu=iyear-(iyear/4)*4
	im=1
	if(id.le.31) return
        im=2
	id=id-31
	iq=28
	if(iuruu.eq.0) iq=29
	if(id.le.iq) return
	id=id-iq
	im=3
	if(id.le.31) return
	id=id-31	!April
	im=4
	if(id.le.30) return
	id=id-30
	im=5
	if(id.le.31) return
	id=id-31
	im=6
	if(id.le.30) return
	id=id-30
	im=7
	if(id.le.31) return
	id=id-31
	im=8
	if(id.le.31) return
	id=id-31
	im=9
	if(id.le.30) return
	id=id-30
	im=10
	if(id.le.31) return
	id=id-31
	im=11
	if(id.le.30) return
	id=id-30
	im=12
	return
	end
	subroutine iyimid_idd2(idd,iy,im,id)
ccc	input iy,im,id,   output idd
	if(im.le.2) then
	idd=(im-1)*31+id
	go to 1
	endif
	iuruu=iy-(iy/4)*4
	id3=31+28
	if(iuruu.eq.0) id3=id3+1
	idd=id3+(im-3)*31+id		!4,6,9,11
	if(im.ge.5) idd=idd-1
	if(im.ge.7) idd=idd-1
	if(im.ge.10) idd=idd-1
	if(im.ge.12) idd=idd-1
    1   continue
       
	return
	end
  	subroutine make_dkwrt_dir(wave_dir2,iy,im,id,dir,ndir)
	character*80 file_wave,text2,file,dir,text,wave_dir2,wave_dir
        character*2 iyc,imc,idc
	character*13 wavename
	if(iy.le.0) return
        write(iyc,'(i2.2)') iy
        write(imc,'(i2.2)') im
        write(idc,'(i2.2)') id
        wave_dir=wave_dir2
        do k=1,80
           if(wave_dir(k:k).eq.' ') then
           nk=k-1
           go to 30
           endif
           enddo
 30        continue
           if(wave_dir(nk:nk).ne.'/') then
              nk=nk+1
              wave_dir(nk:nk)='/'
              endif
           file=wave_dir(1:nk)//iyc//'/test.dat'
           open(7,file=file,status='unknown',err=1)
           close(7)
           go to 2
 1         continue
           file='mkdir '//wave_dir(1:nk)//iyc
           call system(file)
 2         continue
           if(im.le.0) then
	      dir=wave_dir(1:nk)//iyc//'/'
	      ndir=nk+3
	      return
	      endif
          file=wave_dir(1:nk)//iyc//'/'//imc//'/test.dat'
           open(7,file=file,status='unknown',err=3)
           close(7)
           go to 4
 3         continue
           file='mkdir '//wave_dir(1:nk)//iyc//'/'//imc
           call system(file)
 4         continue
	   if(id.le.0) then
ccc     year and month directory only  ***************
	     dir=wave_dir(1:nk)//iyc//'/'//imc//'/'
	     ndir=nk+6
	     return
	     endif

           dir=wave_dir(1:nk)//iyc//'/'//imc//'/'//idc
           ndir=nk+8
         file=dir(1:ndir)//'/test.dat'
cc         WRITE(*,*) FILE
           open(7,file=file,status='unknown',err=5)
           close(7)
           go to 6
 5         continue
           file='mkdir '//dir
c          write(*,*) file
           call system(file)
 6         continue
           ndir=ndir+1
           dir(ndir:ndir)='/'
           return
           end
       subroutine count_text(file1,n)
       character*80 file1
cc	write(*,*) 'count text file=',file1
	n=0
      do j=2,1200
	 n=j-1
         if(file1(j:j).eq.' ') return
         enddo
            return
            end
  
	subroutine idate_iyimid(idate,iy,im,id,ih,imin,sec)
	dimension idate(6)
	iy=idate(1)
	im=idate(2)
	id=idate(3)
	ih=idate(4)
	imin=idate(5)
	sec=float(idate(6))
	return
	end
    	subroutine date_correction(iy,im,id,ih,imin,sec)
cc	write(*,*) 'date correction=',iy,im,id,ih,imin,sec
	 do kk=1,5
	if(sec.lt.-0.0001) then
	 sec=sec+60.0
	 imin=imin-1
	 endif
	 if(sec.ge.60.0) then
	    sec=sec-60.0
	    imin=imin+1
	    endif
	    enddo
cc            write(*,*) 'imin=',imin

CC	    if(imin.ge.0.and.imin.le.59) return
c            write(*,*) imin,sec

	    if(imin.lt.0) then
	       imin=imin+60
	       ih=ih-1
	       endif
	 	    if(imin.ge.60) then
	       imin=imin-60
	       ih=ih+1
	       endif  

c               write(*,*) 'ih=',ih
	       if(ih.lt.0) then
		  ih=ih+24
		  id=id-1
		  endif		  
                  if(ih.ge.24) then
		  ih=ih-24
		  id=id+1
		  endif
		  call iyimid_idd(idd,iy,im,id)
  	          call idd_iyimid(idd,iy,im,id)
		  return
		  end
	subroutine jpoint_iuda(ih,imin,sec,nmem,nsamp,jp)
	if(nsamp.le.0) return
	nn=(3600*ih+imin*60)*nsamp+ifix((sec+1.0)*float(nsamp)+0.1)
	n=nn/nmem
	jp=nn-n*nmem-nsamp


c	isec_share=nmemt/nsa
c	iq=iminsec/isec_share
c	isec_jp=iminsec-iq*isec_share
c	j1=nsa*(isec_jp-1)


        return
	end
	subroutine get_file_in_dir_lib(direct,nnn,file_wave)
	character*80 file_wave(200000),text,text2,file	
	character*80 direct,directa(10000),direct_1(1000)
ccc     subroutine to get file names of win raw, pick,or trg files 	
cc      input directory name eg. ./trg 
cc      out nnn: number of win files
cc      file_wave; fine names together with directory name 
	nc_dir=1
	nc=0
	directa(1)=direct
	do 10 kkk=1,10
	   newd=0
	   do 1 kk=1,nc_dir
	      direct=directa(kk)
	      if(direct(1:1).eq.' ') go to 1
	      do 11 j=1,80
		 if(direct(j:j).eq.' ') then
		    jd=j-1
 		       if(direct(jd:jd).ne.'/') then
		       jd=jd+1
		       direct(jd:jd)='/'
		       endif
	   text='ls '//direct(1:j)//'>print.dat'
cc	   write(*,*) text
	   call system(text)
	   go to 2
	   endif
 11	   continue
 2	   continue
	   open(7,file='print.dat')
	   do 4 k=1,10000
	   jsra=0
	   jper=0
	   je=0
	      read(7,'(a80)',end=44) file
	          do 12 j=1,80
		 if(file(j:j).eq.'.'.and.j.ge.2) jper=j
		 if(file(j:j).eq.'/') jsra=j
		 if(file(j:j).eq.' ') then
		    je=j
		    go to 3
		 endif
 12		 continue
 3		 continue
		 if(jper.eq.0) then
ccc               directory *************
		  nc_dir=nc_dir+1
		  directa(nc_dir)=direct(1:jd)//file(1:je)
		  newd=newd+1
cc	          write(*,*) kkk,kk,newd,file
	          endif
		  if(je.lt.10) go to 4
		  nc=nc+1
		  file_wave(nc)=direct(1:jd)//file(1:je)
cc		  write(*,*) nc,file_wave(nc)
 4	      continue
 44	      continue
	      close(7)
               directa(kk)='  '
 1	      continue
	      if(newd.le.0) go to 15
c	      write(*,*) 'kkk,news=',kkk,newd,nc
c	      read(*,*) iok 
 10	      continue
 15	      continue
	      nnn=nc
cc	      write(*,*) 'nnn=',nnn
	      return
	      end
	      
	subroutine get_jma_win_file(iy,im,id,ih,imin,isec,isec_dur,
     1 dir,trg_dir,file2)
	byte ibytea(300000000),ibyteb(300000000),ibuf_w(2048),ic(4)
	equivalence (ic,ic4)
	character*18 datec
	character*13 text13
  	character*80 dir,file,trg_dir,file2
cc	write(*,*) iy,im,id,ih,imin,isec,isec_dur
	do j=2,80
	   jetrg=j-1
        if(trg_dir(j:j).eq.' ') then
	   if(trg_dir(j-1:j-1).ne.'/') then
	     jetrg=j
	     trg_dir(j:j)='/'
	     endif
	     go to 7
	     endif
	     enddo
 7	     continue
	iyw=iy
	imw=im
	idw=id
	ihw=ih
	iminw=imin
	isecw=isec
	iminsec1=isec+imin*60
	iminsec2=iminsec1+isec_dur

	jt=0
	do 10  mmm=1,4        !read one min file
	call get_win_file_name(iyw,imw,idw,ihw,iminw,dir,file)
cc	write(*,*) iyw,imw,idw,ihw,iminw
cc	write(*,*) 'file=',file
        open(7,file=file,status='unknown',access='direct',
     1  recl=2048,iostat=iostat)
	if(iostat.ne.0) go to 11
ccc   store one minute win format wave data to ibytea(jj)
	jj=0
	do 20 irec=1,10000
	read(7,rec=irec,iostat=iostat) ibuf_w   !read  2048 data
	if(iostat.ne.0) go to 21
	do j=1,2048
	   jj=jj+1
	ibytea(jj)=ibuf_w(j)
	enddo
 20	continue			!irec 
 21	continue
	close(7)
cc	write(*,*) 'one min wave data size=',jj
	nnn=jj
ccc     store trigger wave file to ibyteb(1)-
	j1=1
	do 30 kk=1,1000000
	j2=j1+9
	ic(4)=ibytea(j1)
	ic(3)=ibytea(j1+1)
	ic(2)=ibytea(j1+2)
	ic(1)=ibytea(j1+3)
	if(ic4.le.0.or.ic4.gt.10000000) go to 31
	write(datec,'(6z3)') (ibytea(j),j=j1+4,j2)
cc	write(*,*) 'date=',datec
	read(datec,'(6i3)',err=31) iy,im,id,ih,imin,isec
cc	write(*,'(20z3)')  (ibytea(j),j=1,20)
cc	write(*,'(6i3,2i8)') iy,im,id,ih,imin,isec,ic4
	ims=imin*60+isec
	isa=ims-iminsec1
	if(isa.lt.-1800) isa=isa+3600
	if(isa.gt.1800) isa=isa-3600
	if(isa.ge.0.and.isa.lt.isec_dur) then
	   if(jt.le.0) then  !ftrom trigger time
ccc      define trigger file name 
	   write(text13,110) iy,im,id,ih,imin,isec
 110	   format(5i2.2,'.',i2.2)
	file2=trg_dir(1:jetrg)//text13
cc	write(*,*) 'trg file=',file2
c	write(*,*) iyw,imw,idw,ihw,iminw,isecw, iy,im,id,ih,imin,isec
c	write(*,*) ims,iminsec1,'isa=',isa
	endif
ccc     store triggered wavefor to memory ibyteb********
	   j3=j1+ic4-1
	   do j=j1,j3
	      jt=jt+1   !Pointer of trigger wave
	      ibyteb(jt)=ibytea(j)
	      enddo
cc	write(*,'(6i3,3i8)') iy,im,id,ih,imin,isec,ic4,isa,jt
	      endif
	      j1=j1+ic4
	if(j1.gt.nnn) go to 31
 30	continue
 31	continue
ccc     read next one min file *****************
	ims=iminw*60+60
	isa=iminsec2-ims
	if(isa.lt.-1800) isa=isa+3600
	if(isa.gt.1800) isa=isa-3600
cc	write(*,*) iminsec1,iminsec2,ims,'isa=',isa
	if(isa.lt.0) go to 11
	iminw=iminw+1
	sec=float(isecw)
	call date_correction(iyw,imw,idw,ihw,iminw,sec)
c	WRITE(*,*) 'Date cor',iyw,imw,idw,ihw,iminw,sec
 10	continue
 11	continue
       open(7,file=file2,status='unknown',access='direct',
     1  recl=1)
c	write(*,*) 'dkwrt trig file',file2
c	write(*,*) 'size=',jt
	do irec=1,jt
	   write(7,rec=irec) ibyteb(irec)
	   enddo
	   close(7)
	   go to 33
ccc      check 
	   jj=0
        open(7,file=file2,status='unknown',access='direct',
     1  recl=2048)
	do  irec=1,10000
	read(7,rec=irec,iostat=iostat) ibuf_w   !read  2048 data
	if(iostat.ne.0) go to 51
	do j=1,2048
	   jj=jj+1
	ibytea(jj)=ibuf_w(j)
	enddo
        enddo   !irec 
 51	continue
	close(7)
	write(*,*) 'one min wave data size=',jj
	nnn=jj
ccc     store trigger wave file to ibyteb(1)-
	j1=1
	do  kk=1,isec_dur
	j2=j1+9
	ic(4)=ibytea(j1)
	ic(3)=ibytea(j1+1)
	ic(2)=ibytea(j1+2)
	ic(1)=ibytea(j1+3)
	write(datec,'(6z3)') (ibytea(j),j=j1+4,j2)
cc	write(*,*) 'date=',datec
	read(datec,'(6i3)') iy,im,id,ih,imin,isec
cc	write(*,'(20z3)')  (ibytea(j),j=1,20)
	write(*,'(6i3,2i8)') iy,im,id,ih,imin,isec,ic4
	j1=j1+ic4
	enddo
 33	continue

	return
	end
	subroutine get_win_file_name(iy,im,id,ih,min,dir,file)
  	character*80 dir,file,trg_dir
	character*11 text11
	do j=2,80
	   je=j-1
	   if(dir(j:j).eq.' ') then
	      if(dir(je:je).ne.'/') then
		 je=je+1
		 dir(je:je)='/'
		 endif
		 go to 1
		 endif
		 enddo
 1		 continue
		 write(text11,100) iy,im,id,ih,min
 100		 format(4i2.2,'.',i2.2)
		 file=dir(1:je)//text11
		 return
		 end
        subroutine cal_av_rm_err(nn,av,resa,sd2)
        dimension resa(nn)
        av=0.0
        if(nn.le.0) return
           nc=0 
            z=0.0
            zz=0.0
            do j=1,nn
        z=z+resa(j)
        zz=zz+resa(j)*resa(j)
        enddo
        av=z/float(nn)
        avk=av
cc        write(*,*) 'av=',av
        zz=zz/float(nn)
        sd=sqrt(zz-av*av+0.00000001)
cc        write(*,*) 'sd=',sd
        nc=0
        av=0.0
	sd2=0.0
        do j=1,nn
           z=abs(resa(j)-avk)
           if(z.lt.sd) then
              av=av+resa(j)
              nc=nc+1
	      sd2=sd2+z*z
              endif
           enddo
cc           write(*,*) 'av,nc=',av,nc
           if(nc.ge.1) then
              av=av/float(nc)
	      sd2=sqrt(sd2/float(nc))
              else
                 av=avk
                 endif
                 return
                 end
