	program extrpch4
C....................................................Oct. 2019
C Pch Catalogue100
C
C Wp ~pmtocYR.txt, Wn ~ nmtocYR.txt = magn coord. map
C....................................................Jan 2017
C outfile='extrpch4.txt' for Pch d0+5 days
C....................................................Dec 2016
C
C Extract geom. indices for 2011-2019 at Pch-storms1 (Dst_min <= -50 nT, Dst <=-30 nT) 
C                    Include AE,ap(tau),aa,ap,Wp (potocYR.txt) indices extraction
C select EXTREME values for time to=nint(to_ICME):to+47h
C outfile='extrgicme1.txt'
c
	integer kyr1,kmn1,kdy1,kut1,IM(12)
	dimension iaa(0:366,0:7),iap(0:366,0:7)
c	+,apt(366,0:7)
     +,iae(0:366,0:23),iap2(0:24)
c     +,apt2(0:24)
     +,iaa2(0:24)
     +,iae2(0:24),idst(366,0:23),idst2(0:24),iwp(366,0:23)
C     +,ivn(366,0:23)
	+,iwp2(0:24)
C     +,ivn2(0:24)
C//    +,idst2(0:47)
	CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
	CHARACTER*65 infaa,infpch,infae,inpdst,infwp
C	+,infvn,infapt
	+,outfile,txt
	CHARACTER*2 AMN,AYR,ADY,AUT
	CHARACTER*4 AYEAR
C//	LOGICAL aavar,apvar,aptvar,aevar,dstvar !'TRUE'('FALSE') for distirbed (quiet)
	DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/ 
C
	CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT
C

	infaa='d:/apkp/AA/aaYEAR' !aa, ap
	inpdst='d:/apkp/dst/dstYEAR'	 ! Dst
	infapt='d:/apkp/APT/aptYEAR' ! ap(tau)
	infae='d:/apkp/AE/aeYEAR'	 ! AE
C	infvp='c:/msdev/projects/vsigma/potocYR.txt' ! Vp
C	infvn='c:/msdev/projects/vsigma/notocYR.txt' ! Vn
	infvp='c:/msdev/projects/wtec/pmtocYR.txt' ! Vp
	infvn='c:/msdev/projects/wtec/nmtocYR.txt' ! Vn
	write(*,*) ' Enter YEAR:'
	read(*,*) IYEAR
CREM	outfile='extrgicme1.txt'  ! 	total file for 1999-2016 =376 ICME List t0+47 hrs
	outfile='extrgicme4.txt' !total file for 1999-2016 =376 ICME List t0+47hrs
C Enter ICME list for 1999-2016

   36 format(1X,I4,A1)
   37 format(20X,24(I4))	 ! AE
  137 format(13X,2(1X,I3)) ! ap,aa
  237 format(6X,8(1X,F5.1))  ! ap(tau)
   34 format(6X,24(1X,I4))  ! Vp, Vn
      icntot=0
C Cycle on years 1999-2016
	iyear=1998   !TEMP
C	iyear=2015	 !Temp
	call blet4(iyear,AYEAR)
c//	outfile(6:9)=AYEAR
		OPEN(16,FILE=OUTFILE,access='APPEND')
	nicmes=0	   ! total ICME events 199-2016/04
C
C start cycle on year: =====================================
C
C+	do 400 KK=1,52	! 1964-2015
	do 400 KK=1,18	! 1999-2016
C-	kk=1 !Temp
C//	if (kk.gt.1) ndyr0=ndnyr
	iyear=iyear+1
Crem	iyear=2015 !TEMP
	call blet4(iyear,AYEAR)
	ryear=float(iyear)
	AYR=AYEAR(3:4)
	read(ayr,*) ryr
	iyr=int(ryr)
		if(int(ryr/4.)*4.eq.iyr) THEN 
		 im(2)=29
	ndnyr=366
	                          ELSE
	IM(2)=28
	ndnyr=365
	ENDIF
C       !==============================================
	do nn=1,ndnyr
	do k=0,24
	iap(nn,k)=0
	iaa(nn,k)=0
	apt(nn,k)=0.
	iae(nn,k)=0
	enddo
	enddo
C Reading aa & ap file
	infaa(14:17)=AYEAR
	OPEN(10,FILE=infaa)
	do idd=1,ndnyr
	do ii=0,7
	read(10,137,end=139) iap(idd,ii),iaa(idd,ii)
	enddo
	enddo
  139	close(unit=10)
C
C Reading apt file
	infapt(16:19)=AYEAR
	OPEN(10,FILE=infapt)
	do idd=1,ndnyr
	read(10,237,end=238) (apt(idd,k),k=0,7)
	enddo
  238	close(unit=10)
C
C Reading AE file
C
	infae(14:17)=AYEAR
	OPEN(10,FILE=infae)
	do idd=1,ndnyr
	read(10,37,end=339) (iae(idd,ii),ii=0,23)
	enddo
  339	close(unit=10)
C
C Reading Dst file
C
	inpdst(16:19)=AYEAR
	OPEN(10,FILE=inpdst)
	do idd=1,ndnyr
	read(10,37,end=439) (idst(idd,ii),ii=0,23)
	enddo
  439	close(unit=10)
C
C Reading Vp file
C
	infvp(29:30)=AYR
	OPEN(10,FILE=infvp)
	read(10,*) txt ! extra line
	do idd=1,ndnyr
	read(10,34,end=539) (ivp(idd,ii),ii=0,23)
	enddo
  539	close(unit=10)
C
C Reading Vn file
C
	infvn(29:30)=AYR
	OPEN(10,FILE=infvn)
	read(10,*) txt ! extra line
		do idd=1,ndnyr
	read(10,34,end=639) (ivn(idd,ii),ii=0,23)
	enddo
  639	close(unit=10)

C 1247 DST9901*13 17  DST9901*13 23 -112.  DST9901*15 20   52 --
   35 format(I4,3(1X,I2),I2)
C
C??	DO 399 imn=1,12  !
C??	call blet2(imn,AMN)
C+
	inficme='ListICMEn1.txt'  ! 1999-2016/04
C
	OPEN(15,file=inficme,action='READ')
C
C Cycle on Dst storm events
C
	iev=0
C
	do 28 jj=1,600           ! ICME list for 1999-2016
  228 read(15,35,end=27,err=27) kyear0,kmn0,kdy0,kut0,kmin
C	+,kyr2,kmn2,kdy2,kut2,ihtot
C	 mindst=nint(dstmin)
	kyr0=kyear0-(kyear0/100)*100
	uth=float(kut0)+kmin/60.
	kut0=nint(uth)
	if (kut0.eq.24) then
	kut0=0
	kdy0=kdy0+1
	if (kdy0.gt.im(kmn0)) then
	kdy0=kdy0-im(kmn0)
	kmn0=kmn0+1
	endif
	endif
	kyr2=kyr0
	kmn2=kmn0
C
	kdy2=kdy0+2	! t0+47 hrs
CREM      kdy2=kdy0	! t0+23 hrs
	kut2=kut0-1	 	! t0+47 hrs
CREM	kut2=kut0+23	! t0+23 hrs
	ihtot=48  ! t0+47 hrs
CREM	ihtot=24 ! t0+23 hrs
	if (kut2.lt.0) then	  ! t0+47 hrs
	kut2=kut2+24		  ! t0+47 hrs
	kdy2=kdy2-1			  ! t0+47 hrs
	endif				  ! t0+47 hrs
CREM	if (kut2.ge.24) then	  ! t0+23 hrs
CREM	kut2=kut2-24		  ! t0+23 hrs
CREM	kdy2=kdy2+1			  ! t0+23 hrs
CREM	endif				  ! t0+23 hrs
	if (kdy2.gt.im(kmn2)) then
	kdy2=kdy2-im(kmn2)
	kmn2=kmn2+1
	if (kmn2.gt.12) then  ! Only for Dec.31, 2005
	kmn2=12
	kdy2=31
	kut2=23
	ihtot=kut2-kut0
	endif
	                   endif
	maxap=0
	aptmax=0.
	maxaa=0
	maxae=0
	maxap=0
	maxvp=0
	maxvn=0
	mindst=1000.
	IF (kyear0.lt.iyear) goto 28
      IF (kyear0.gt.iyear) goto 27
C	 
227	kyr1=kyr0
	kmn1=kmn0
	kdy1=kdy0
C/	kut1=kut0-6        ! -6 prec. hours
	kut1=kut0		   ! ICME start
C/	if (kut1.lt.0) then
C/	kut1=kut1+24
C/	kdy1=kdy1-1
C/	if (kdy1.lt.1) then
C/	kmn1=kmn1-1
C/	kdy1=im(kmn1)
C/	endif
C/	endif
	call blet2(kyr1,AYR)
 	call blet2(kmn1,AMN)
	call blet2(kdy1,ADY)
	call blet2(kut1,AUT)


C//	IF (kyr1.eq.iyr) goto 128
	           
C//	if (kyr1.gt.iyr) GOTO 27
C
  128	iev=iev+1
	lda1=ndoy(kyr1,kmn1,kdy1)
	lda2=ndoy(kyr2,kmn2,kdy2)
	ldatot=lda2-lda1+1       ! number of days during ICME event & t0+47 hrs
	lda=lda1
	yrlda=ryear+float(lda1)/ndnyr

	kyr=kyr1
	kmn=kmn1
	kdy=kdy1
		if (ldatot.eq.1) then ! only 1 day DSt storm
 	k1=kut1
	k2=kut2
                         else
	k1=kut1
	k2=23
	    endif 
C

 328	CONTINUE
C
  177	do ii=0,21,3
      i1=ii/3
	iaa2(ii)=iaa(lda,i1)
	iap2(ii)=iap(lda,i1)
	apt2(ii)=apt(lda,i1)
	enddo
	if (lda.lt.ndnyr) then
	iaa2(24)=iaa(lda+1,0)
	iap2(24)=iap(lda+1,0)
	apt2(24)=apt(lda+1,0)
	else
	iaa2(24)=iaa(lda,7)
	iap2(24)=iap(lda,7)
	apt2(24)=apt(lda,7)
	endif
C Interpolate aa, ap, aptau for 1h:
	do ii=1,22,3
	delaa=(iaa2(ii+2)-iaa2(ii-1))/3.
	iaa2(ii)=nint(iaa2(ii-1)+delaa)
      iaa2(ii+1)=nint(iaa2(ii-1)+delaa*2.)
	delap=(iap2(ii+2)-iap2(ii-1))/3.
	iap2(ii)=nint(iap2(ii-1)+delap)
      iap2(ii+1)=nint(iap2(ii-1)+delap*2.)
	delapt=(apt2(ii+2)-apt2(ii-1))/3.
	apt2(ii)=apt2(ii-1)+delapt
      apt2(ii+1)=apt2(ii-1)+delapt*2.
	enddo
C
	do ii=0,23
	iae2(ii)=iae(lda,ii)
      idst2(ii)=idst(lda,ii)
	ivp2(ii)=ivp(lda,ii)
	ivn2(ii)=ivn(lda,ii)
	enddo
	if (lda.lt.ndnyr) then
	iae2(24)=iae(lda+1,0)
	idst2(24)=idst(lda+1,0)
	ivp2(24)=ivp(lda+1,0)
	ivn2(24)=ivn(lda+1,0)
	                 else
	iae2(24)=iae(lda,23)
	idst2(24)=idst(lda,23)
	ivp2(24)=ivp(lda,23)
	ivn2(24)=ivn(lda,23)
	endif

C
C
C        
C	
C
C output for quiet conditions:
	do iut=k1,k2
C//	call blet2(iut,AUT)
	icntot=icntot+1
CREM	write(*,33) AYR,AMN,ADY,AUT,iae2(iut),idst2(iut),iaa2(iut)
CREM	+,iap2(iut),apt2(iut),icntot
CREM	write(16,33) AYR,AMN,ADY,AUT,iae2(iut),idst2(iut),iaa2(iut)
CREM	+,iap2(iut),apt2(iut),icntot
      if(maxap.lt.iap2(iut)) maxap=iap2(iut)
	if(maxaa.lt.iaa2(iut)) maxaa=iaa2(iut)
	if(maxae.lt.iae2(iut)) maxae=iae2(iut)
	if(aptmax.lt.apt2(iut)) aptmax=apt2(iut)
	if(maxvp.lt.ivp2(iut)) maxvp=ivp2(iut)
	if(maxvn.lt.ivn2(iut)) maxvn=ivn2(iut)
	if(mindst.gt.idst2(iut)) mindst=idst2(iut)
	enddo
C
  	if (lda.lt.lda2) then  !<<<<<<<<<<<<<<<<
	lda=lda+1
	kut=0
	kdy=kdy+1
C
	k1=0
	if (lda.eq.lda2) then  !------
	k2=kut2
	else				   !------
	k2=23
	endif				   !------
C
	if (kdy.gt.im(kmn)) then !===========
	kdy=1
	kmn=kmn+1
	if (kmn.gt.12) goto 27
                             endif  !=========
	 if (lda.le.ndnyr) goto 328
	endif !<<<<<<<<<<<<<
	
	write(*,33) yrlda,AYR,AMN,ADY,AUT,maxae,mindst,maxaa
	+,maxap,aptmax,maxvp,maxvn,icntot,iev
	write(16,33) yrlda,AYR,AMN,ADY,AUT,maxae,mindst,maxaa
	+,maxap,aptmax,maxvp,maxvn,icntot,iev

  28	CONTINUE
   27 close(unit=15)
c
C
  399 continue	! end of monthly cycle
   33 format(F9.3,1X,3A2,1X,A2,1X,4(1X,I5),1X,F5.1,2(1X,I5),1X,I6,1X,I4)
   31 format(1X,A4,2(1X,I5))
C
C	write(*,31) AYEAR,icnt6
C	write(16,31) AYEAR,icnt6
	nicmes=nicmes+iev

C      close(unit=16)
CC end Cycle on mn --------------------------------------
C
Ctemp	IF (iyear.gt.2015) exit
  400	CONTINUE ! next year
C	
	write(*,31) AYEAR,icntot,iev,nicmes
C-	write(16,31) AYEAR,icntot,iev
	close(unit=16)
	goto 900
    2 CONTINUE
           write(*,*) 'INPUT FILE IS NOT IN YOUR DIRECTORY '
  900 pause ' '
      STOP
      END

