      program addeqjap
C.............................................................................Oct 2022
C  set of EQs over Japan
C....................................................
C..Add EQ daily Mmax, Dkm, Mav, glat,glon set, for the day of Kpmax...........May 2022
C
      CHARACTER*64 outfile,infileq
	+,infssn
		CHARACTER*4 AYEAR
	DIMENSION eqs(1500),glats(1500),glons(1500),kmds(1500),kms(1500)
	+,kuts(1500),kmms(1500)
C==     +,irnd(39)
C==	 DATA irnd/ 3, 5, 33, 35, 38, 45, 48, 52, 55, 68
C==     +, 72, 78, 94, 107, 121, 130, 154, 157, 172, 183
C==     +, 197, 200, 209, 226, 231, 233, 238, 241, 248, 256
C==     +, 262, 271, 285, 292, 303, 306, 337, 343, 345/

C
      infssn='SSNd.txt' ! check /exclude SSN=0
      outfile='cateqjap.txt'
	infileq='d:\Apkp\Equak\Obninsk\eqYEAR.txt'
C
       OPEN(10,file=infssn,action='READ')
       OPEN(11,file=outfile,access='APPEND')
   30 format(I4,2(1X,I2),1X,F8.3,1X,I4)  !SSN2
   33 format(1X,F8.3,1X,I3,2(1X,I2),1X,F5.3,1X,I2,1X,F5.3,1X,I4)
   31 format(I4,A64)
C
	ndy=0
C
	DO 41 iyear=1995,2020
C-      DO 41 iyear=1997,2020

C
	iyr=iyear-iyear/100*100
		ryr=float(iyr)
		if(int(ryr/4.)*4.eq.iyr) THEN 
	ndr=366
	                          ELSE
	ndr=365
	ENDIF

C++
CC Cycle on lda:
C
          DO 42 lda=1,ndr
C==	lda=irnd(iev)
C==	call submmdd(lda,iyr,imn,idy)
   21	read(10,30,end=50) jyear,jmn,jdy,ryear,issn
      jyr=jyear-jyear/100*100
      ldaj=ndoy(jyr,jmn,jdy)
      if ((jyear.eq.iyear).and.(ldaj.eq.lda)) then
	goto 22
	                  else
	goto 21
	               endif
C
   22	 CONTINUE
C==      if (issn.eq.0) GOTO 42  ! exclude SSN =0
C
	do k=1,90
	eqs(k)=0.0
	glats(k)=0.0
	glons(k)=0.0
	kmds(k)=0
	kms(k)=0
	enddo
	eqav=0.0
	aeqmax=0.0
	knum=0
	kk=0
C-	if (iyear.eq.2021) goto 2 ! no EQ file for 2021
C
C 
       call blet4(iyear,AYEAR)
	infileq(25:28)=AYEAR
 	OPEN(12,file=infileq,action='READ')
	kk=0
  1	read(12,*,end=55) kyear,kmn,kdy,dut,glat,glon,kmd,eqmag,knum
      kyr=kyear-kyear/100*100
	ldak=ndoy(kyr,kmn,kdy)
C==      if (kmn.gt.imn) goto 55
       if (ldak.gt.lda) then 
	backspace(12)
	 goto 55
	               endif
C==	      if ((kmn.eq.imn).and.(kdy.gt.idy)) goto 55
      if ((kmn.eq.jmn).and.(kdy.eq.jdy)) then
	if ((glat.lt.30.).or.(glat.gt.47.)) goto 1
      if ((glon.lt.134.).or.(glon.gt.151.)) goto 1
	kutmmss=nint(dut)
	kuth=int(kutmmss/10000.)
	rmmss=kutmmss-kuth*10000.
	kmm=int(rmmss/100.)
	rss=rmmss-kmm*100.
      rmm=kmm+rss/60.
	kmm=nint(rmm)
	kk=kk+1
	glats(kk)=glat
	glons(kk)=glon
	kmds(kk)=kmd
	eqs(kk)=eqmag
	kms(kk)=knum
	kuts(kk)=kuth
	kmms(kk)=kmm
	eqav=eqav+eqmag
	endif
C
	 goto 1
C
   55	close(12)
      if (kk.gt.0) then 
	eqav=eqav/kk
	               else
	goto 42
	             endif
C
C select daily peak
	aeqmax=0.
        do jj=1,kk
	if (aeqmax.lt.eqs(jj)) then 
	aeqmax=eqs(jj)
	knum=kms(jj)
	jm=jj
	endif
	 enddo

C Output
   2  ndy=ndy+1
      write(*,32) ryear,lda,jmn,jdy,kuts(jm),kmms(jm),glats(jm)
	+,glons(jm),kmds(jm),aeqmax,eqav,kk,knum,issn,ndy
	write(11,32) ryear,lda,jmn,jdy,kuts(jm),kmms(jm),glats(jm)
	+,glons(jm),kmds(jm),aeqmax,eqav,kk,knum,issn,ndy
   42 continue
C
   41 continue
   32 format(1X,F8.3,1X,I3,2(1X,I2),2(1X,I3),2(1X,F6.2),1X,I4,2(1X,F5.1)
     +,4(1X,I5))
   50	close(unit=11)
	pause ' '
	STOP
	END