      program tcvs24m2
C.......................................................July 2015
C Vs effectiveness at EQ-quiet and EQ-storm 
C Moving window: = 24 grids around EQ epicenter
C Select area +-10 deg: 4 lats>glat,+4lats<glat; 2longs>glong+2 longs <glong
C
C Source Vs-map 'd:\web\drs\YR\MN\tiMNDDUT.jYR'  
	dimension ify(71,0:72),IM(12)
	character*64 eqlist,infile,outfile
	character*4 AYEAR
	character*2	AYR,AMN,ADY,AUT,BYR,BMN,BDY
	character*1 sq
	real ut1
	common/const/	umr
	DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/ 
C

   33	format(73(1X,I4))
   35	format(1X,A4,2A2,1X,F5.2,1X,3(1X,F6.2),2(1X,F5.2),1X,F5.1)
   36	format(A1)
   37	format(F2.0)
C
	 PI=ATAN(1.0)*4.
        UMR=PI/180.
   98 write(*,*) ' Enter s or q for EQ-storm or EQ-quiet: '
	read(*,36) sq
	outfile='EQvspnx2.txt'	  ! cnt Vs-pos, Vs-neg in EQ-24 grids +1 h
	outfile(7:7)=sq
	OPEN(15,file=outfile,access='APPEND')
C=   	write(*,*) ' Enter delta_t <dhr>: '
C=	read(*,37) dhr
	dhr=0.0				   ! const
C   	sq='s' ! storm EQ	 !TEMP
c	if (sq.eq.'s') then
c	write(*,*) ' EQ storm'
c	write(15,*) ' EQ storm'
c	else
c	write(*,*) ' EQ quiet'
c	write(15,*) ' EQ quiet'
c	endif
	eqm=6.0 ! starting M of EQ
C
C+  100	eqlist='exteq6s.txt'
  100	eqlist='eqvspnq.txt'
	eqlist(7:7)=sq
	mcntp=0
	mcntn=0
	mcntm=0
	mmm=0
	mrest=0
c=	if (eqm.lt.8.0) then
c=	eqmup=eqm+0.5
c=	          else
	eqmup=10.0
c=	             endif
	open(13,file=eqlist,action='READ')
  101	read(13,35,end=27,err=26) AYEAR,BMN,BDY,UT1
 	+,glat,glon,depth,eqmag,hlt,xhi
	mmm=mmm+1
	if (AYEAR.eq.'2016') goto 27   ! end of EQs
Crem	if (xhi.gt.90.) goto 101  ! exclude night ; results for day
Crem	if (xhi.lt.90.) goto 101  ! exclude day ; results for night
C
c	mcntm=mcntm+1
	BYR=AYEAR(3:4)
	read(BYR,*) ryr
	iyr=int(ryr)
		if(int(ryr/4.)*4.eq.iyr) THEN 
		 im(2)=29
	ndnyr=366
	                          ELSE
	IM(2)=28
	ndnyr=365
	ENDIF
C
	read(BMN,*) rmn
	imn=int(rmn)
	read(BDY,*) rdy
	idy=int(rdy)
	IF (ut1.lt.0.) THEN
	kut1=0
	                ELSE
C	ut2=ut1+0.005   ! nearest h UT
c	ut2=ut1+1.005   ! nearest+1 h UT
	ut2=ut1+0.005+dhr   ! nearest+dhr, h UT
	kut1=int(ut2)+1

	rlt=ut2+glon/15.
	if (rlt.ge.24.) rlt=rlt-24.
	if (rlt.lt.0.) rlt=rlt+24.

C Add calculation of XHI to exclude night or day:
	lda=ndoy(iyr,imn,idy)
C=	call soco(lda,rlt,glat,glon,dec,xhi,sr,ss)
Crem		if (xhi.gt.90.) goto 101  ! exclude night ; results for day
C=	 	if (xhi.lt.90.) goto 101  ! exclude day ; results for night
	if ((eqmag.ge.eqm-0.001).and.(eqmag.lt.eqmup)) then
	goto 177
	            else
	mrest=mrest+1
c+	goto 101
	endif
  177	mcntm=mcntm+1
	if (kut1.ge.24) then
	kut1=kut1-24
	idy=idy+1
	   if (idy.gt.im(imn)) then
	   idy=1
	   imn=imn+1
      	if (imn.gt.12) then
      	imn=1
      	iyr=iyr+1
      	endif
	   endif
      endif
	ENDIF
C
	call blet2(iyr,AYR)
	call blet2(imn,AMN)
	call blet2(idy,ADY)
	call blet2(kut1,AUT)
C
C
C
	infile='d:\web\drs\YR\MN\tiMNDDUT.jYR' ! EQ-day
C
	infile(12:13)=AYR
	infile(15:16)=AMN
	infile(20:21)=AMN
	infile(28:29)=AYR
C
	infile(24:25)=AUT
	infile(22:23)=ADY

	OPEN(12,FILE=infile,action='READ') ! EQ day 
	qlat=-90.0
	    do ilat=1,71  !read  for EQ-day
	qlat=qlat+2.5
	read(12,33,end=233,err=26) (ify(ilat,k),k=0,72) ! Vs index
	if ((glat.ge.qlat).and.(glat.lt.qlat+2.5)) iqlat=ilat
	enddo
  233	close(unit=12)  !

C 
	qlon=-185.
	xlon=glon
	if (glon.gt.180.) xlon=glon-360.
	do ilon=0,71
	qlon=qlon+5.
	if (qlon.ge.180.) qlon=qlon-360.

	if ((xlon.ge.qlon).and.(xlon.lt.qlon+5.0)) iqlon=ilon
	enddo
	iqlon_1=iqlon-1
	if (iqlon_1.lt.0) iqlon_1=iqlon_1+70
	iqlon1=iqlon+1
	if (iqlon1.gt.72) iqlon1=iqlon1-72
	iqlon2=iqlon+2
	if (iqlon2.gt.72) iqlon2=iqlon2-72
C
	icntp=0
	icntn=0
	icntm=0 !total =24 grids
      	do ilat=iqlat-3,iqlat+4
	    if ((ilat.ge.iqlat-1).and.(ilat.lt.iqlat+3)) then
	icntm=icntm+4							  !ilat-1,ilat,ilat+1,ilat+2
	if (ify(ilat,iqlon_1).le.-2) icntn=icntn+1
	if (ify(ilat,iqlon).le.-2) icntn=icntn+1
	if (ify(ilat,iqlon1).le.-2) icntn=icntn+1
	if (ify(ilat,iqlon2).le.-2) icntn=icntn+1

	if (ify(ilat,iqlon_1).ge.2) icntp=icntp+1
	if (ify(ilat,iqlon).ge.2) icntp=icntp+1
	if (ify(ilat,iqlon1).ge.2) icntp=icntp+1
	if (ify(ilat,iqlon2).ge.2) icntp=icntp+1
                                else
	icntm=icntm+2							  !ilat-3,-2,+3,+4
	if (ify(ilat,iqlon).le.-2) icntn=icntn+1
	if (ify(ilat,iqlon1).le.-2) icntn=icntn+1
	if (ify(ilat,iqlon).ge.2) icntp=icntp+1
	if (ify(ilat,iqlon1).ge.2) icntp=icntp+1
                   endif
	enddo
C Output results:
Crem      write(*,34) AYEAR,BMN,BDY,UT1,glat,glon,depth,eqmag,hlt,xhi,kut1
Crem	+,icntp,icntn,icntm
Crem	write(15,34) AYEAR,BMN,BDY,UT1,glat,glon,depth,eqmag,hlt,xhi,kut1
Crem	+,icntp,icntn,icntm
	mcntp=mcntp+icntp
	mcntn=mcntn+icntn
	GOTO 101 ! Next EQ


   27	if (mcntm.gt.0.) then
      perp=float(mcntp)/(mcntm*24.)
      pern=float(mcntn)/(mcntm*24.)
	endif
C=   	write(*,31) dhr,mcntp,mcntn,mcntm,perp,pern
C= 	write(15,31) dhr,mcntp,mcntn,mcntm,perp,pern
           	write(*,31) eqm,mcntp,mcntn,mcntm,perp,pern
	write(15,31) eqm,mcntp,mcntn,mcntm,perp,pern
C=	dhr=dhr+1.0
	eqm=eqm+0.5
	close(unit=13)
C=	if (dhr.lt.13.) goto 100
C=	if (eqm.lt.8.1) goto 100
      close(unit=15)
C
   34	format(1X,A4,2A2,1X,F5.2,1X,3(1X,F6.2),2(1X,F5.2),1X,F5.1
     +,4(1X,I4))
   31	format(1X,F5.1,3(1X,I6),2(1X,F6.3))
   	goto 400
   26	write(*,*) 'Error in input file=',eqlist
  400	continue ! 
	pause ' '
	stop
	END
C
