program whstrms5i
C......................................................................Dec. 2021
C Med(-15d) mag coord
C......................................................................Aug. 2018
C From WptemYR.txt - GIM-TEC in magn. coords
C
C Wpt-STORM SELECT STORM pwriods= Wp-max>3.0 WITH 12-HRS Wp-indices
C ......Interpolate from 12 2-hrsUT to 24 1=hrsUT Wp values...................
C Wp>=3.0 = storm onset, storm end. Power p=Wpmax*Tdur>=200.
C T.L. Gulyaeva                                                        May, 2008
      DIMENSION IM(12),DA(366,0:24)
	DIMENSION XDEV(0:24)
	CHARACTER*35 INFILE
	CHARACTER*24 OUTFILE
	CHARACTER*128 txt
	CHARACTER*6 BL,YRINP,RE(366),D1,D2,SYEAR,AME
	+,D1J(30),D2J(30),DMJ(30)
	CHARACTER*2 AYR
	REAL YEAR
C D1J,IT1J,DMJ,IUTMAXJ,AMAXJ,D2J,IT2J,IUTTJ,ADSTJ,POWERJ
	DIMENSION IT1J(30),IT2J(30),IUTMAXJ(30),AMAXJ(30)
	+,IUTTJ(30),ADSTJ(30),POWERJ(30)
	COMMON DA
	integer irun
c	COMMON irun
	common /bl1/D1J,IT1J,DMJ,IUTMAXJ,AMAXJ,D2J,IT2J,IUTTJ,ADSTJ,POWERJ
      DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
		BL=" "
	irun=0
 100	WRITE(*,'(A\)') ' ENTER YEAR    :'
      READ (*,*) SYEAR
	READ (SYEAR,*) YEAR
	YEAR1=YEAR

cc	INFILE='wptec99.tau'   ! tau
cc    INFILE='wptec99.txt'   ! txt
cc	INFILE='wptj99.tau'	   ! mag-lat tau
CmagILE='wptj99.txt'	   ! mag-lat txt

ct      INFILE='wptk99.txt'	   ! maglat & MLT txt
crem      INFILE='wtuh99.txt'	   ! maglat & MLT txt
c///      INFILE='wutec99.txt'	   ! glat -85:5:85, glon -180:15:165
C//	INFILE='c:\msdev\projects\Wtec\wjtec99.txt'	   ! glat -85:5:85, glon -180:15:165
	INFILE='d:\web\grif\datw\wptem14.txt'
C1	OUTFILE='wptcstau.out'	! tau
cc	OUTFILE='wptcsta6.out'  ! mag-lat  tau
CC	OUTFILE='wptcstxt.out'	! txt
cc	OUTFILE='wptgstx6.out'	! txt Wpmax>6.0 duration >=4h GLAT
cc	OUTFILE='wptjstau.out'  ! mag-lat  tau
vv		OUTFILE='wptjsta6.out'  ! mag-lat  tau
C-	OUTFILE='wptjstxt.out'  ! mag-lat  txt
Cmag	OUTFILE='wptjstx6.out'  ! mag-lat  txt Wp>6.0 duration >=4h
crem	OUTFILE='wptkstx6.out'  ! mag-lat  txt
	OUTFILE='whstrms5i.out'  ! mag-lat  txt

	AYR=SYEAR(3:4)
	AME(1:)=AYR
c///      INFILE(6:7)=AYR		   ! mag-lat  tau,txt
c//      INFILE(29:30)=AYR		   ! mag-lat  tau,txt
      INFILE(23:24)=AYR		   ! JPL IONEX-lat  tau,txt
c/		INFILE(5:6)=AYR		! tau, txt GLAT
      IF ((YEAR/4.).EQ.IFIX(YEAR/4.)) THEN 
	IM(2)=29 
	ELSE 
	IM(2)=28
	ENDIF

      IDNR=365
      IF (IM(2).EQ.29) IDNR=366
C	PAUSE ' '
	do k=0,24
	xdev(k)=0.
	enddo
	OPEN(111,FILE=INFILE)
	READ (111,17,END=18,ERR=99) TXT ! extra lone
   17 format(A128)
C iut0 = starting UT for 24 hrs interpolation
c	READ (111,11,END=2,ERR=99) YRINP,iut0,txt ! Line of UT hrs
c   11  FORMAT (A6,1X,I4,A62)
C day-to-day cycle
   18	iend=0
	iflag=0
	DO 15 J=1,IDNR
		wp0=xdev(11) ! remember for 24 hrs UT interpolation
	READ (111,12,END=13,ERR=99) YRINP,(XDEV(k),k=0,24)
   12  FORMAT (A6,24(1X,F4.1),3X,F4.1)
C	IF (YRINP(3:3).EQ.BL) YRINP(3:3)="0"
C	IF (YRINP(5:5).EQ.BL) YRINP(5:5)="0"
	goto 11
   13	iflag=1
   11   RE(J)=YRINP
cc	da(j,24)=xdev(12)
C Interpolation from 12 to 24 hourly values
cc	if (j.eq.1) wp0=xdev(0)
C
	DO 14 K=0,24
C iut=0:
	DA(j,k)=xdev(k)
   14 CONTINUE
C
	iend=iend+1
	if (iflag.eq.1) exit
   15 CONTINUE
C
   2   CONTINUE
       CLOSE (UNIT=111)
C
C	OPEN(112,FILE=OUTFILE,ACCESS='APPEND',STATUS='OLD')	
		OPEN(112,FILE=OUTFILE,ACCESS='APPEND')	
C
	jj=0
      DO 810 J=1,IEND		! ================================
	AMAX=0.
	CNT=0
	IFL=0
      DO K=0,23
	IF (DA(J,K).GT.AMAX) THEN
	AMAX=DA(J,K)
	IUTMAX=K
	ENDIF
	ENDDO
CC      IF (AMAX.LT.5.) goto 810 
	 if(j.eq.iend) then
	call subout(1,irun)
	if (jj.eq.2) then
	call subout(2,irun)
	endif
	if (jj.eq.3) then
	 call subout(3,irun)
	endif
	endif

Temp	      IF (AMAX.LT.6.) goto 810  ! Collect all cases of Wp>4
c+	      IF (AMAX.LE.4.) goto 810  ! Collect all cases of Wp>4
	      IF (AMAX.LE.3.0) goto 810  ! Collect all cases of Wp>3.0
C GO TO SUBROUTINE DETECTING <STORM-START;STORM-END> USING Wp>=4.
  710 CONTINUE
      D1='' 
      D2=''
	IT1=0
	IT2=0
	ADST=0.
	IDAY=J

	CALL DETECT(IDAY,RE,D1,IT1,D2,IT2,IUTT,ADST,IUTMAX,POWER)
C	IF (ADST.EQ.0.) GOTO 810
	AVER=DA(J,24)
Temp	if (power.lt.100.) goto 810
Temp	if (iutt.lt.3) goto 810 ! collect storms and sub-storms with DUR>3
C	if (AMAX.lt.5.0) goto 810 ! collect storms and sub-storms
	jj=jj+1
C D1J,IT1J,DMJ,IUTMAXJ,AMAXJ,D2J,IT2J,IUTTJ,ADSTJ,POWERJ
  807	D1J(jj)=D1
	IT1J(jj)=IT1
	DMJ(jj)=RE(j)
	IUTMAXJ(jj)=IUTMAX
	AMAXJ(jj)=AMAX
	D2J(jj)=D2
	IT2J(jj)=IT2
	IUTTJ(jj)=IUTT
	ADSTJ(jj)=ADST
	POWERJ(jj)=POWER
	do l=jj+1,30
	powerj(l)=0.
	enddo
	if (jj.eq.1) goto 810
C
	if ((POWERJ(jj).eq.POWERJ(jj-1)).and.(ADSTJ(jj).eq.ADSTJ(jj-1)))
	+ then ! ++
	 if (iday.lt.iend) then
	 	 goto 810 
	else
	     goto 809
	endif
	                           else		 ! ++
		if (jj.eq.2) then     ! 
	kkk=1
	call subout(kkk,irun)
	jj=1
	goto 807
	            else		   ! 1
	goto 809
	    endif				   ! 1
	endif								! ++
  809	continue
       jjj=jj-1
	wpmax=0.
	do k=1,jjj
	if (wpmax.lt.amaxj(k)) then
	wpmax=amaxj(k)
	kkk=k
	endif
	enddo
	call subout(kkk,irun)
	jj=1
	goto 807
C	PAUSE ' '
  810 CONTINUE		 ! ===================
C	ENDDO
C
      CLOSE (UNIT=112)
      YEAR=YEAR+1.0
      IF (YEAR.LT.2022) GOTO 100
   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,POW)
	DIMENSION DA(366,0:24)
	CHARACTER*6 D1,D2,RE(366)
	INTEGER C1,C2
	COMMON DA
	
C START-OF-STORM DETECTING================
      J1=JDAY
C
  970    C1=-1  
	K1=23
      IF (J1.EQ.JDAY) K1=IMX
      K=K1
 1000 CONTINUE
C PRINT K;DEV(J1,K)
c++      IF (DA(J1,K).GE.4.) THEN 		 ! LE.4. ???
      IF (DA(J1,K).GE.3.0) THEN 		 ! LE.3.0 ???
	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 PRINT D1,IT1 

C END-OF-STORM DETECTING =================
      J2=JDAY
 1140  C2=-1
	 K1=0
1150  IF (J2.EQ.JDAY) K1=IMX

      DO K=K1,23
C     PRINT K;DEV(J2,K)
c++       IF (DA(J2,K).GE.4.) THEN 		 ! LE.4. ???
       IF (DA(J2,K).GE.3.0) THEN 		 ! LE.3.0 ???
       C2=K 
	 GOTO 1200
	               ELSE
 1190 ISTR2=K-1  
      GOTO 1220
	             ENDIF
 1200 CONTINUE
	 ENDDO
 1210 IF (C2.EQ.23) THEN 
      J2=J2+1  
	GOTO 1140
	ENDIF
 1220 IF (ISTR2.GE.0) GOTO 1240
      J2=J2-1 
	ISTR2=23
 1240 D2=RE(J2)
      IT2=ISTR2 
C PRINT D2,IT2

C DERIVE STORM DURATION=========== 
 1290 I1=IMX-ISTR1
      IF (J1.LT.J2) I1=24-ISTR1
      I2=ISTR2+1
      IF (J1.EQ.J2) I2=ISTR2-IMX+1
C PRINT "i1=";I1,"i2=";I2
      II=I1+I2
 1350 IF ((J2-J1).GT.1) THEN
         GOTO 1360 
	   ELSE 
	   GOTO 1370
	ENDIF
 1360  II=II+(J2-J1-1)*24
 1370 IUT=II			   ! Storm duration, hrs
C
C Calculation of WP=max, WP-MEAN, POW=Wpmax*T: FOR THE STORM DURATION =========
	wpmax=0.
	do n=j1,j2
	do k=0,23
	if (wpmax.lt.da(n,k)) wpmax=da(n,k)
	enddo
	enddo
	POW=wpmax*iut    ! Power of storm
C
      SAP=0.
      L1=ISTR1 
	 L2=ISTR2 
	 JJ=J1
      IF (J2.GT.J1) L2=23
	
	      DO 116 L=L1,L2
      SAP=SAP+DA(JJ,L)
  116    CONTINUE

      IF (J2.EQ.J1) GOTO 1520
1450  L1=0  
      L2=ISTR2 
	 JJ=J2
      
	      DO 117 L=L1,L2
      SAP=SAP+DA(JJ,L)
  117    CONTINUE

      IF (J2.EQ.(J1+1)) GOTO 1520
      L1=0 
	 L2=23 
	  JJ=J1
 1490  JJ=JJ+1
	      DO 118 L=L1,L2
      SAP=SAP+DA(JJ,L)
  118    CONTINUE

 1510 IF (JJ.LT.(J2-1)) GOTO 1490
 1520     AADST=SAP/II+.5  
 1560	CONTINUE
      RETURN
	END
	subroutine subout(kk,irun)
	DIMENSION IT1J(30),IT2J(30),IUTMAXJ(30),AMAXJ(30)
	+,IUTTJ(30),ADSTJ(30),POWERJ(30)
	 integer irun
	character*6 D1J(30),D2J(30),DMJ(30),YYMMDD
		character*4 AYEAR1,AYEAR2,AYEAR3
	character*2 AYR1,AMN1,ADY1,AYR2,AMN2,ADY2,AYR3,AMN3,ADY3
	COMMON /bl1/D1J,IT1J,DMJ,IUTMAXJ,AMAXJ,D2J,IT2J,IUTTJ,ADSTJ,POWERJ
C
cc	COMMON irun
C
Temp	if(iuttj(kk).lt.3) goto 171	! Wp storm duration >=3 h
	if(iuttj(kk).le.3) goto 171	! Wp storm duration >=4 h
Temp	if(POWERJ(kk).lt.100.) goto 171	! Power>=100.
c--      if(AMAXJ(kk).lt.6.) goto 171	! Wpmax>=6.
      if(AMAXJ(kk).lt.5.) goto 171	! Wpmax>=5.
C--      if(AMAXJ(kk).lt.5.3) goto 171	! Wpmax>=5.3
      irun=irun+1
	YYMMDD=D1J(kk)
	AYR1=YYMMDD(1:2)
	AMN1=YYMMDD(3:4)
	ADY1=YYMMDD(5:6)
	read(ayr1,*) xyr
	iyr=int(xyr)
	AYEAR1='0000'
	AYEAR1(3:4)=AYR1
	if (iyr.lt.30) then
	AYEAR1(1:2)='20'
	else
      AYEAR1(1:2)='19'
	endif
C
	YYMMDD=DMJ(kk)
	AYR2=YYMMDD(1:2)
	AMN2=YYMMDD(3:4)
	ADY2=YYMMDD(5:6)
	read(ayr2,*) xyr
	iyr=int(xyr)
	AYEAR2='0000'
	AYEAR2(3:4)=AYR2
	if (iyr.lt.30) then
	AYEAR2(1:2)='20'
	else
      AYEAR2(1:2)='19'
	endif
C
	YYMMDD=D2J(kk)
	AYR3=YYMMDD(1:2)
	AMN3=YYMMDD(3:4)
	ADY3=YYMMDD(5:6)
	read(ayr3,*) xyr
	iyr=int(xyr)
	AYEAR3='0000'
	AYEAR3(3:4)=AYR1
	if (iyr.lt.30) then
	AYEAR3(1:2)='20'
	else
      AYEAR3(1:2)='19'
	endif
C
C		WRITE(112,170) D1J(kk),IT1J(kk),DMJ(kk),IUTMAXJ(kk),AMAXJ(kk)
C	+,D2J(kk),IT2J(kk),IUTTJ(kk),ADSTJ(kk),POWERJ(kk),irun
		WRITE(*,170) irun,AYEAR1,AMN1,ADY1,IT1J(kk),AYEAR2,AMN2,ADY2
	+,IUTMAXJ(kk),AMAXJ(kk),AYEAR3,AMN3,ADY3,IT2J(kk),IUTTJ(kk)
     +,ADSTJ(kk),POWERJ(kk)
		WRITE(112,170) irun,AYEAR1,AMN1,ADY1,IT1J(kk),AYEAR2,AMN2,ADY2
	+,IUTMAXJ(kk),AMAXJ(kk),AYEAR3,AMN3,ADY3,IT2J(kk),IUTTJ(kk)
     +,ADSTJ(kk),POWERJ(kk)

  170	FORMAT(I4,1X,A4,2A2,I3,2X,A4,2A2,I4,F6.1,2X,A4,2A2,I3,1X,I4
     +,2(1X,F6.1))

  171	return
	end
C