      program addeqrnd
C.............................................................................Oct 2022
C Random set of EQs (except SSN = 0)
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(90),glats(90),glons(90),kmds(90),kms(90)
	+,kuts(90),kmms(90),irnd(39)
	 DATA irnd/ 3, 5, 33, 35, 38, 45, 48, 52, 55, 68
     +, 72, 78, 94, 107, 121, 130, 154, 157, 172, 183
     +, 197, 200, 209, 226, 231, 233, 238, 241, 248, 256
     +, 262, 271, 285, 292, 303, 306, 337, 343, 345/

C
      infssn='SSNd.txt' ! check /exclude SSN=0
      outfile='cateqrnd.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
C	DO 41 iyear=1995,2020
       iyear=2021
C-      DO 41 iyear=1997,2020

C
	iyr=iyear-iyear/100*100
C++
CC Cycle on random events:
C
          DO 42 iev=1,39
	lda=irnd(iev)
	call submmdd(lda,iyr,imn,idy)
   21	read(10,30,end=50) jyear,jmn,jdy,ryear,issn
      if ((jyear.eq.iyear).and.(jmn.eq.imn).and.(jdy.eq.idy)) then
	goto 22
	                  else
	goto 21
	               endif
C
   22	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
      if (kmn.gt.imn) goto 55
	      if ((kmn.eq.imn).and.(kdy.gt.idy)) goto 55
      if ((kmn.eq.imn).and.(kdy.eq.idy)) then
	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,imn,idy,kuts(jm),kmms(jm),glats(jm)
	+,glons(jm),kmds(jm),aeqmax,eqav,kk,knum,issn,ndy
	write(11,32) ryear,lda,imn,idy,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,4(1X,I2),2(1X,F6.2),1X,I4,2(1X,F4.1)
     +,4(1X,I4))
   50	close(unit=11)
	pause ' '
	STOP
	END