      program dstphase22
C.........................................................Jul. 2021
C
C Specify storm pgase: Fi=0 (t0=SSC), Fi=1 (storm peak), Fi=0(t0+47h)
C
      CHARACTER*64 infile,infssc,outfile
	CHARACTER*2 AMN,ADY,AUT,AMM
	CHARACTER*4	AYEAR
	DIMENSION idat(9,0:47),res(9,0:47),idst(0:71) ! 1h dst for 9 storms SC23
     +,iYRs(9),iMNS(9),iDYS(9),iUTs(9),IM(12)
       DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/

C
      infile='dst22d.txt'
	outfile='dstphase22.txt'
	infssc='storms22c.txt'
	AMM='00'
	OPEN(12,file=infssc,action='READ')
	DO ii=1,8
	read(12,33) iyear,imn,idy,iut
	iyrs(ii)=iyear
	imns(ii)=imn
	idys(ii)=idy
	iuts(ii)=iut
	ENDDO
	close(unit=12)
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)
   35	FORMAT(72(1X,I4))
C
      DO iev=1,8
	read(10,35) (idst(k),k=0,71)
	do k=0,47
	idat(iev,k)=idst(k+24)
	enddo
	ENDDO
	close(unit=10)
C
      DO 50 iev=1,8
C Select peak DST
      mindst=100.
Crem	do k=0,47
	do k=0,23
	if (idat(iev,k).lt.mindst) then
	mindst=idat(iev,k)
	k1=k
	endif
	enddo
C Phase fi=[0,1]:
	dmin=float(k1)
       DO k=0,k1+1
	FI=float(k)/dmin  ! Main Phase
	res(iev,k)=FI
	 enddo
	DO k=k1+1,47
	FI=(float(k)-48.)/abs(dmin-48.)	! Recovery phase
	res(iev,k)=FI
       ENDDO
   50	CONTINUE
C
C Output results
C
        DO iev=1,8
	write(*,34) iev
	write(11,34) iev
	iyr=iyrs(iev)
      imn=imns(iev)
	idy=idys(iev)
	iut=iuts(iev)
	call blet4(iyr,AYEAR)
	ryear=float(iyear)
c
 	z1=iyr/4.0
      jz=int(z1)*4
      IF(jz.EQ.iyr) THEN
               IM(2)=29
	ndnyr=366
        ELSE
                IM(2)=28
	  ndnyr=365
	       ENDIF

	  DO 75 ii=0,47
      call blet2(imn,AMN)
	call blet2(idy,ADY)
	call blet2(iut,AUT)
	write(*,36) AYEAR,AMN,ADY,AUT,AMM,res(iev,ii),idat(iev,ii)
	write(11,36) AYEAR,AMN,ADY,AUT,AMM,res(iev,ii),idat(iev,ii)
C+	if (res(iev,ii).eq.1) then 
C+	fi=-1.
C+	write(*,36) AYEAR,AMN,ADY,AUT,AMM,fi,idat(iev,ii)
C+	write(11,36) AYEAR,AMN,ADY,AUT,AMM,fi,idat(iev,ii)
C+      endif
	iut=iut+1
	IF (iut.gt.23) THEN
	iut=0
	idy=idy+1
	if (idy.gt.im(imn)) then
	idy=1
	imn=imn+1
	endif
	ENDIF
   75	CONTINUE
	ENDDO
	close(unit=11)
   34 FORMAT(1X,'fi',I1)
   36 FORMAT(1X,A4,2A2,1X,2A2,1X,F6.3,1X,I4)
	pause ' '
	STOP
	END

      	
	 