      program extwcpiaeqr
C......................................................................Oct 2022
C Random EQs set at SSN > 0
C
C..Add peak Wc(pia)..for EQs at SSN0 (Kp0=1183, Kp3=609)...............June 2022
C JPL Wc indices 
C |Wc| peak from -1d:tEQ:+1d (SSN = 0) =from 4 W indices: at dlati=2.5o, dlong=5o around EQ=Mmax
C Wc = combined  with weights from 4 W-indices surrounding EQ daily peak
C 
C 
C
      CHARACTER*64 outfile,infileq
	CHARACTER*2 AYR,AMN,ADY,AUT,AMM,SMN,SDY,SUT
	DIMENSION wp(0:71)
C
C
      infileq='cateqrnd.txt' ! Random set for SSN>0  1995-2020
	outfile='addwcpiaeqr.txt'
C
       OPEN(11,file=infileq,action='READ')
       OPEN(12,file=outfile,access='APPEND')
   30 format(1X,F8.3,1X,I3,4(1X,I2),2(1X,F6.2),1X,I4,2(1X,F4.1)
     +,4(1X,I4))
   33 format(3I2,24(1X,F4.1))
C
	DO k=0,71
	wp(k)=0.
	ENDDO
	ndy=0
C
Cnorm	DO 41 iev=1,786	   ! random EQ(SSN>0) norm
	DO 41 iev=1,787	     ! TEMP
    	read(11,30,end=50) ryear,lda,imn,idy,iuteq,imm,glats,glons,kmds
  	+,aeqmax,eqav,kk,knum,issn,ndyeq
      iyear=int(ryear)
	if (iev.lt.787) goto 41  !TEMP +Tohoku EQ
C
      DO j=0,71
	wp(j)=0.0
	enddo
C
C
C 	 produce lag for day 
C++ READ 1h Wp
		z1=iyear/4.0
      jz=int(z1)*4
      IF(jz.EQ.iyear) THEN
	ndnr=366
		  ELSE
                ndnr=365
		       ENDIF

	iyr=iyear-iyear/100*100
	call blet2(iyr,AYR)
	call blet2(imn,amn)
	call blet2(idy,ady)
	call blet2(iuteq,aut)
	call blet2(imm,amm)
C
	lda0=lda-1
	if (lda0.lt.1) GOTO 43
	call submmdd(lda0,iyr,imn0,idy0)
	call blet2(imn0,smn)
	call blet2(idy0,sdy)
C
  43    DO iut=0,23
	call blet2(iut,sut)
	call subinterpw(glats,glons,ayr,smn,sdy,sut,res)
	wp(iut)=res
	ENDDO
C
      DO iut=0,23
      call blet2(iut,sut)
      call subinterpw(glats,glons,ayr,amn,ady,sut,res)
	wp(iut+24)=res
	ENDDO
C
      lda1=lda+1
	if (lda1.gt.ndnr) GOTO 44
	call submmdd(lda1,iyr,imn1,idy1)
	call blet2(imn1,smn)
	call blet2(idy1,sdy)
      DO iut=0,23
	call blet2(iut,sut)
	call subinterpw(glats,glons,ayr,smn,sdy,sut,res)
	wp(iut+48)=res
	ENDDO
  44    niteq=iuteq+24
       wpeq=wp(niteq)
C Wc-peak_pre
      wc_pre=0.
	awc_pre=0.
	it_pre=0
       kt1=niteq-24
	kt2=niteq+23
C
	do k=niteq-1,kt1,-1
      abswc=abs(wp(k))
	if ((abswc.gt.awc_pre)) then
	awc_pre=abswc
	wc_pre=wp(k)
	it_pre=k
            endif
	enddo
	jtw_pre=it_pre-niteq
c-	if ((wc_pre.eq.0.0).and.(it_pre.eq.0)) jtw_pre=0
C
C Wc-preak_nex
      wc_nex=0.0
	awc_nex=0.0
	it_nex=0
      	do k=niteq+1,kt2
      abswc=abs(wp(k))
	if ((abswc.gt.awc_nex)) then
	awc_nex=abswc
	wc_nex=wp(k)
	it_nex=k
            endif
	enddo
	jtw_nex=it_nex-niteq
C Output
      write(*,32) ryear,lda,amn,ady,aut,amm,glats,glons,kmds
 	+,aeqmax,eqav,kk,knum,wpeq,jtw_pre,wc_pre,jtw_nex,wc_nex,issn,ndyeq

	write(12,32) ryear,lda,amn,ady,aut,amm,glats,glons,kmds
 	+,aeqmax,eqav,kk,knum,wpeq,jtw_pre,wc_pre,jtw_nex,wc_nex,issn,ndyeq

C
   41 continue
   32 format(1X,F8.3,1X,I3,4(1X,A2),2(1X,F6.2),1X,I4,2(1X,F4.1)
     +,2(1X,I4),1X,F5.2,2(1X,I3,1X,F5.2),2(1X,I4))
   50	close(unit=11)
	close(unit=12)
	pause ' '
	STOP
	END