program Apostorm15s
C.......................................................Apr 2023
C Apo proxy (7h smoothed median centered on 4th day) Apo peak > 100 nT
C
C.......................................................July 2022
C Storms from 1h Hpo/Apo data at Ap > 15 nT, Apmax >=48 nT (Kp >= 5.0)
C DETECT STORM ONSET as soon as Apo>15 nT; Storm-end when Dstaver<15 nT
C This version to detect storm periods using ONLY input Hpo/Apo-indices
C T.L. Gulyaeva                                                
c DD1,IIT1,DDM,IIUTMAX,AAMAX,DD2
      DIMENSION IM(12),DA(366,0:24)
	DIMENSION XDEV(0:23)
	CHARACTER*30 INFILE
	CHARACTER*4 AYEAR,BYEAR
		CHARACTER*2 BYR,BMN,BDY
	CHARACTER*25 OUTFILE
	CHARACTER*10 YRINP,RE(366),D1,D2,DM,DD1,DD2,DDM
	+,DD1_pred,DDM_pred,DD2_pred
	REAL YEAR
	COMMON DA
      DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
C-	INFILE='D:\APKP\Hpo\Hpo19952023.txt'
	INFILE='D:\APKP\Hpo\ApoYEARs.txt'
Cre	OUTFILE='Apostorm15.txt'
	OUTFILE='Apostorm15s.txt'   ! Smoothed Apo max >100 nT
C
	icnt=0
	icnt100=0
 100	 DO 101 iyear=1995,2023
C 100	 DO 101 iyear=2023,2023
C      WRITE(*,'(A\)') ' ENTER YEAR    :'
C      READ (*,*) SYEAR
C	READ (SYEAR,*) YEAR
       call blet4(iyear,AYEAR) 
	infile(16:19)=AYEAR
	BYEAR=AYEAR
	OPEN(11,FILE=INFILE,action='READ')
C
        YEAR=float(iyear)
	YEAR1=YEAR
C	iyear=int(YEAR)
C
      IF ((YEAR/4.).EQ.IFIX(YEAR/4.)) THEN 
	IM(2)=29 
	IDNR=366
	ELSE 
	IM(2)=28
	IDNR=365
	ENDIF
C Put initials: D1,IT1,DM,IDAY,AVER,AMAX,IUTMAX,D2,IT2,IUTT,ADST
	DD1=''
	DD1_pred=DD1
	DD2=''
	DD2_pred=DD2
	DDM=''
	DDM_pred=DDM
      AAMAX_pred=0.
	IIT1=24
	IIT2=24
	IIDAY=0
	AAVER=0.
	AAMAX=0.
	IIUTMAX=24
	IIUTT=0
	AADST=0.
	IFL=0
C End of initials
   31  FORMAT (I4)
C-    1	read(11,31,end=2) kyear
C-      if (kyear.ne.iyear) then
C-	goto 1
C-	else
C-	backspace(11)
C-	endif
C
C Start of data input for current year
C
	DO 15 J=1,IDNR
C Calculation of daily mean:
	DAV=0.
C-	DO 14 K=0,23
C-	READ (11,32,END=2,ERR=99) YRINP,XDEV(k)
	YRINP='YEAR MN DD'
		READ (11,32,END=2,ERR=99) BYR,BMN,BDY,(XDEV(k),k=0,23)
	BYEAR(3:4)=BYR
	YRINP(1:4)=BYEAR
	YRINP(6:7)=BMN
	YRINP(9:10)=BDY
C  
C   32  FORMAT (I4,3(1X,I2),36X,F4.0)
C- 32  FORMAT (A10,43X,F4.0)
  32  FORMAT (3A2,24(1X,F4.0))

	RE(J)=YRINP
	do 14 k=0,23
	DA(J,K)=xdev(K)
	DAV=DAV+XDEV(K)
   14 CONTINUE
	DA(J,24)=DAV/24.
   15 CONTINUE

   2   CONTINUE
       CLOSE (UNIT=11)
C
		AAMAX_print=0.
C
		OPEN(12,FILE=OUTFILE,ACCESS='APPEND')	
		      DO 810 J=1,IDNR
C ++++++ Record results for preceding solution:
C	D1,IT1,DM,IDAY,AVER,AMAX,IUTMAX,D2,IT2,IUTT,ADST
	IIFL=IFL
	       IF ((J.GT.1).AND.(IFL.EQ.1)) THEN
	DD1=D1
	DD2=D2
	DDM=DM
	IIT1=IT1
	IIT2=IT2
	IIDAY=IDAY
	AAVER=AVER
	AAMAX=AMAX
	IIUTMAX=IUTMAX
	IIUTT=IUTT
	AADST=ADST
	       ENDIF

      AMAX=0.
	CNT=0
	
      DO K=0,23
      IF (DA(J,K).GT.AMAX) THEN 
	AMAX=DA(J,K) 
	IUTMAX=K
	ENDIF
	ENDDO

C      IF (AMAX.LE.-50.) THEN 
C-      IF (AMAX.GT.48.) THEN 
      IF (AMAX.GT.100.) THEN 
	  IFL=1
	  GOTO 710
	ELSE
	  IFL=0
        IF (IIFL.EQ.1) THEN
	  GOTO 790
	            ELSE
	  GOTO 805
	           ENDIF
	ENDIF
	
C GO TO SUBROUTINE DETECTING <STORM-START;STORM-END> USING KPYR.D27
  710 CONTINUE
      D1='' 
      D2=''
	IT1=0
	IT2=0
	ADST=0.
	IDAY=J

	CALL DETECT(IDAY,RE,D1,IT1,D2,IT2,IUTT,ADST,IUTMAX)
	AVER=DA(J,24)
	DM=RE(J)
C

	IF ((D1.EQ.DD1).AND.(D2.EQ.DD2).AND.(IIUTT.EQ.IUTT)) THEN
	IF (J.EQ.IDNR) GOTO 800
	   IF (AMAX.GT.AAMAX) THEN

	       GOTO 805
	                      ELSE
	D1=DD1
	D2=DD2
	DM=DDM
	IT1=IIT1
	IT2=IIT2
	IDAY=IIDAY
	AVER=AAVER
	AMAX=AAMAX
	IUTMAX=IIUTMAX
	IUTT=IIUTT
	ADST=AADST
	IFL=IIFL
	      
	   ENDIF
	                                                ELSE
	if (iifl.eq.1) then
	goto 790
	 else
	      IF (J.EQ.IDNR) THEN 
		  GOTO 800
	                     ELSE
	      GOTO 805
	                     ENDIF
	endif
	                                                ENDIF
	GOTO 805
  790	CONTINUE
      if (iiutt.lt.3) goto 178
crem      WRITE(*,170) DD1,IIT1,DDM,IIUTMAX,IIDAY,AAVER,AAMAX,DD2
crem 	&,IIT2,IIUTT,AADST
crem 	WRITE(12,170) DD1,IIT1,DDM,IIUTMAX,IIDAY,AAVER,AAMAX,DD2
crem	&,IIT2,IIUTT,AADST
	       if (DD1.eq.DD1_pred) then              !<<<<<<<<<<<<<<<<<<<<<<
      	if ((AAMAX_pred.gt.AAMAX).and.(AAMAX_pred.gt.0.)) then !++++++++
	goto 178
	      else									  
	DD1_pred=DD1											!++++++++
	DD2_pred=DD2
	DDM_pred=DDM
	iit1_pred=iit1
	iit2_pred=iit2
	IIUTMAX_pred=IIUTMAX
	IIUTT_pred=IIUTT
	AAMAX_pred=AAMAX
	goto 178
	    endif                                                 !++++++++
                                 else					!<<<<<<<<<<<<<<<<<<<<<<
	    if (AAMAX_print.ne.AAMAX_pred) then				   !-----------
C	if (iiutt_pred.lt.3) goto 181
		icnt=icnt+1
C?	if (AAMAX.le.-100.) icnt100=icnt100+1
C-	if (AAMAX.ge.48.) icnt100=icnt100+1
	if (AAMAX.ge.100.) icnt100=icnt100+1
	     AAMAX_print=AAMAX_pred
      WRITE(*,170) ICNT,DD1_pred,IIT1_pred,DDM_pred,IIUTMAX_pred
     +,AAMAX_pred,DD2_pred,IIT2_pred,IIUTT_pred
 	WRITE(12,170) ICNT,DD1_pred,IIT1_pred,DDM_pred,IIUTMAX_pred
	+,AAMAX_pred,DD2_pred,IIT2_pred,IIUTT_pred
  181 continue   
	     else												!-------------
	DD1_pred=DD1
	DD2_pred=DD2
	DDM_pred=DDM
	iit1_pred=iit1
	iit2_pred=iit2
	IIUTMAX_pred=IIUTMAX
	IIUTT_pred=IIUTT
	AAMAX_pred=AAMAX
	goto 177
	       endif											!-----------
	                                    
                                endif                       !<<<<<<<<<<<<<<<<<<<<<
  177	continue
      if (iiutt_pred.lt.3) goto 178
      icnt=icnt+1
C??      if (AAMAX.le.-100.) icnt100=icnt100+1
      if (AAMAX.ge.48.) icnt100=icnt100+1
	AAMAX_print=AAMAX_pred
      WRITE(*,170) ICNT,DD1_pred,IIT1_pred,DDM_pred,IIUTMAX_pred
     +,AAMAX_pred,DD2_pred,IIT2_pred,IIUTT_pred
 	WRITE(12,170) ICNT,DD1_pred,IIT1_pred,DDM_pred,IIUTMAX_pred
	+,AAMAX_pred,DD2_pred,IIT2_pred,IIUTT_pred
crem  170	FORMAT(1X,A10,I3,2X,A10,I3,1X,I3,2F6.0,2X,A10,I3,1X,I4,F6.0)
  170	FORMAT(1X,I4,1X,A10,I3,2X,A10,I3,F6.0,2X,A10,I3,1X,I4)
  178	if (j.lt.idnr) goto 805
  800	CONTINUE
c	icnt=icnt+1
c	WRITE(*,170) ICNT,D1,IT1,DM,IUTMAX,AMAX,D2,IT2,IIUTT
c	WRITE(112,170) ICNT,D1,IT1,DM,IUTMAX,AMAX,D2,IT2,IIUTT
C	PAUSE ' '
  805	IIFL=IIFL+1
	
  810 CONTINUE
C
C
      CLOSE (UNIT=12)
C      YEAR=YEAR+1
C      IF (YEAR.LT.2020) GOTO 100
  101	 CONTINUE
	write(*,*) icnt100
	pause ' '
   99	CONTINUE
C	CLOSE (UNIT=112)
      STOP
	END
C******************************************************************************
C      SUBROUTINE TO DETECT <STORM-START;STORM-END>
	SUBROUTINE DETECT(JDAY,RE,D1,IT1,D2,IT2,IUT,AADST,IMX)
	DIMENSION DA(366,0:24)
	CHARACTER*10 D1,D2,RE(366)
	INTEGER C1,C2
	COMMON DA
	
C START-OF-STORM DETECTING================
      J1=JDAY
	J2=JDAY
	 C2=-1
       K2=IMX
C  PRINT IMX,DEX 
  970    C1=-1  
	K1=23
      IF (J1.EQ.JDAY) K1=IMX
      K=K1
 1000 CONTINUE
C      WRITE(*,*) J1,K,DA(J1,K)
      IF (DA(J1,K).GE.15.) THEN 
	C1=K  
	GOTO 1030
	                     ELSE
      ISTR1=K+1 
	GOTO 1060
	ENDIF
 1030 K=K-1
      IF (K.GE.0) GOTO 1000
      IF (C1.GE.0) THEN 
	J1=J1-1 
	GOTO 970
	ENDIF

 1060  IF (ISTR1.LT.24) GOTO 1080
      J1=J1+1 
	 ISTR1=0
 1080  D1=RE(J1)
      IT1=ISTR1 
C      WRITE(*,*) D1,IT1 
C	PAUSE ' '
C DERIVE DURATION OF THE MAIN PHASE OF STORM:=========== 
 1290 IAV=IMX-ISTR1+1+(JDAY-J1)*24
	DSTAV=0.
1100	CONTINUE
	DO JA=J1,JDAY
	   IF (JA.EQ.J1) THEN
	     I0=ISTR1
	                 ELSE
	     I0=0
	   ENDIF 

	   IF (JA.LT.JDAY) THEN
	    IK=23
C	    JB=JA+1
	                   ELSE
	    IK=IMX
C	    JB=JDAY
	   ENDIF
		dstav=dstav+DSTAVER(JA,JA,I0,IK)
	ENDDO
C DETECTING OF END-OF-STORM =================
C Check of Dstaver<-30. nT:
 1300	CONTINUE
C UT-end of storm:
 1150 IT2=K2
	JEND=J2
      K2=K2+1
      IF (K2.GT.23) THEN
	J2=J2+1
	K2=0
	ENDIF
C+ 1160	IF (DA(J2,K2).LE.-25.) THEN
 1160	IF (DA(J2,K2).GE.15.) THEN
      DSTAV=DSTAV+DA(J2,K2)
	IAV=IAV+1
	GOTO 1150
                          	ELSE
		     AADST=DSTAV/IAV
C STORM DURATION=========== 
	     IUT=IAV
	                         ENDIF
 1240 D2=RE(JEND)
    
C PRINT D2,IT2
        RETURN
	END
C
C AVERAGE FOR SELECTED DAYs AND TIMEs
	REAL FUNCTION DSTAVER(N1,N2,M1,M2)
	DIMENSION DA(366,0:24)
	COMMON DA
	DO N=N1,N2
	DO M=M1,M2
	DSTAVER=DSTAVER+DA(N,M)
	ENDDO
	ENDDO
	RETURN
	END
C*****************************************************