      program dstphaseav
C.........................................................Oct. 2021
C
C Specify storm pgase: Fi=0 (t0=SSC), Fi=1 (storm peak), Fi=0(t0+47h)
C
      CHARACTER*64 infile,outfile
	CHARACTER*2 AMN,ADY,AUT,AMM,fi
	CHARACTER*4	AYEAR
	DIMENSION idat(9,0:47),res(9,0:47) ! 1h dst for 9 storms SC23
     +,thr1(11),thr2(11),avdst(20),dst_sum(20),dst_cnt(20)
C
	infile='dstphase.txt'
      outfile='dstphasav.txt'
	AMM='00'
C
        thr1(1)=0.0
	thr2(1)=-1.0
        do k=1,10
	thr1(k+1)=thr1(k)+0.1
      thr2(k+1)=thr2(k)+0.1
	  enddo
	thr2(11)=-0.00001
C
C
      OPEN(10,file=infile,action='READ')
      OPEN(11,file=outfile)  ! Phase FI of storm & Dst from to(SSC) to to+47
   33	FORMAT(6X,I4,2I2,1X,I2)
   34 FORMAT(1X,A2,I1)
   35	FORMAT(72(1X,I4))
   36 FORMAT(1X,A4,2A2,1X,2A2,1X,F6.3,1X,I4)
   37 FORMAT(1X,F5.2,1X,F6.1)
C
C
      DO nt=1,9
C Select peak DST
	read(10,34) fi,iev
	do ii=0,47
	read(10,36,end=40) AYEAR,AMN,ADY,AUT,AMM,res(nt,ii),idat(nt,ii)
	enddo
      enddo
   40	close(unit=10)
C
       DO 123 nn=1,9      ! Training set
C
        do k=1,20
	avdst(k)=0.0
	dst_sum(k)=0.
	dst_cnt(k)=0.
	enddo
C Phase fi=[0,1]:
C
      DO nt=0,47 ! aver (fi)					   !>>>>>>>>>>>>>>>>>>>>>
C
        IF (res(nn,nt).ge.0.) THEN
	do k=1,10
	if ((res(nn,nt).ge.thr1(k)).and.(res(nn,nt).le.thr1(k+1))) then
	dst_sum(k)=dst_sum(k)+idat(nn,nt)
	dst_cnt(k)=dst_cnt(k)+1.0
	endif
	enddo
	dst_sum(11)=res(nn,nt)  ! fi=-1.0
	dst_cnt(11)=1.0
	                     ELSE
	do k=1,10
	if ((res(nn,nt).gt.thr2(k)).and.(res(nn,nt).le.thr2(k+1))) then
	dst_sum(k+10)=dst_sum(k+10)+idat(nn,nt)
	dst_cnt(k+10)=dst_cnt(k+10)+1.0
	endif
	enddo
	                   ENDIF
	ENDDO                    !>>>>>>>>>>>>>>>>>
C
C Output results
C
	write(*,34) fi,nn
	write(11,34) fi,nn
        DO 75 ifi=1,20	  !Phase
	if (ifi.lt.11) then
	thr=(thr1(ifi)+thr1(ifi+1))/2.0
	else
	thr=(thr2(ifi-10)+thr2(ifi-9))/2.0
	endif
      if (dst_cnt(ifi).gt.0.) then
	avdst(ifi)=dst_sum(ifi)/dst_cnt(ifi)
	endif
c
	write(*,37) thr,avdst(ifi)
	write(11,37) thr,avdst(ifi)
   75	CONTINUE		 !1h data
  123 CONTINUE
	close(unit=11)
	pause ' '
	STOP
	END

      	
	 