C........................................................Sep 2023
C........................................................July 2023
C -VBz and Vsw, for Catalogue of AE>1000 + Apo . 
C 
C 
C
C  Search new t0 = t(AE>1000) peak 
C retrive new t0-24h ... t0(AE>1000) peak ... t0+47 h ===total 72 hrs
C
C Extract Apo index and GEC  for -12h<t0< +36h for 223 Apo > 100 nT events 1995-2021
C or test storms 224-230 events 2022-2023
C
C........................................................Mar. 2023
C 
C t0 at Apo smoothed>=100nT storm peak
C
C Extract Apo smoothed-index and GEC  for -24h<t0< +48h for Apo events 
C 
C Select period of 72hrs starting from t0-24h, ending t0+48h; t0 @ AEmax_hour
C 
      PROGRAM devbzmn12
C 
C----------------------------------------------------------------
C
	DIMENSION vsw(0:72) ! 1h indices for 72h file 
	+,vbz(0:72),ivsw(0:71),ivbz(0:71)
	+,jDYR(0:72),jDMN(0:72),jDDY(0:72),jDUT(0:72)
C
	CHARACTER*4 AYEAR,CYEAR
      CHARACTER*128 INFILE,INFHSS,outvsw,outvbz,txt 
c--	+,outfile  !<= add GECave to infhss
	CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
      CHARACTER*2 AYR,AMN,ADY,AUT
c--	+,CMN,CDY,CUT
C
      ICALLS=0
	CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT
C
      infile='d:\Apkp\OMNI\swYEAR.txt'  ! OMNI data
C  	infhss='catapo6mn.txt'  ! Apo(tau6) monthly peak
  	infhss='cataeapo3d.txt'  ! Apo(AE1000) peak
C
c--	outfile='catgecavmn90.txt'  !Catalogue for 1995-2023 + GECav at pre-event 24hrs 
	outvsw='vswapoae3d.txt' ! Vsw profiles: t0-24...t0+47 (72 hourly )
	outvbz='vbzapoae3d.txt' ! VBz profiles: t0-24...t0+47 (72 hourly)
C
C Start 
C
    1	CONTINUE
C
C    7	OPEN(11,FILE=OUTFILE,access='APPEND')
C      OPEN(12,FILE=OUTGEC,access='APPEND')
C	OPEN(13,FILE=OUTDEV,access='APPEND')
C--    7	OPEN(11,FILE=OUTFILE)
    7	CONTINUE
      OPEN(12,FILE=OUTVSW)
	OPEN(13,FILE=OUTVBZ)
C
C Ae storms from 1995 to 2017
C
   28 format(1X,F8.3,1X,I3,3(1X,I2),1X,F5.0)
   29 format(1X,F8.3,1X,I3,3(1X,A2),1X,F5.0,1X,F5.2,1X,I4)
  197 format(3X,2A2,1X,A2,10X,24(1X,I4))
C
	open(10,file=infhss,action='READ')
C
C-	open(15,file=infile,action='READ')	 !+++++++++++++++++++++++++!!!!!!!!!!
C
C Cycle on events ====================================================
C
Cnorm	nev=945 1995- 2023.02
C       	nev=342 ! total Apo monthly peak	1995-2023.06
              	nev=276 ! total Apo(AE1000) peak	1995-2017
C
	DO 177 iev=1,nev
C
      read(10,28,end=27) xyear,lda1,imn1,idy1,iut1,aep
C	if (iev.lt.290) GOTO 177  ! test outputTEMPPPPPPPPP
C
        call blet2(imn1,AMN)
      call blet2(idy1,ADY)
	call blet2(iut1,AUT)

	iutcnt=0
	iyear=int(xyear)
	call blet4(iyear,AYEAR)
	infile(16:19)=AYEAR
C
	OPEN(15,file=infile,action='READ')	 !+++++++++++++++++++++++++!!!!!!!!!!
	do k=1,9
	read(15,*) txt  ! extra lines
	enddo
C
	iyr=iyear-iyear/100*100
	iut0=iut1-12
	lda0=lda1-1 ! pre-event day
C
C
       nhr=0
	khr=0
C
C--------------------------------------------------
  222	CONTINUE
C
		if(int(xyear/4.)*4.eq.iyear) THEN 
	ndnyr=366
	                          ELSE
	ndnyr=365
	ENDIF
C
C------------------------------------------------------
C
   36 format(I4,1X,I3,1X,I2,1X,F4.2)	 ! hourly gec
   37	format(A4,2(1X,A2),7X,F5.1,16X,F5.0)
C
C
C read all years:
C
C
1135       call blet4(iyear,AYEAR)
        AYR=AYEAR(3:4)
C
C============================================
       iflag=0
	napo=72
C
        bz_pre=0.
	vsw_pre=0.
C
C-  333 read(15,36,end=102) jyear,jDMN(72),jDDY(72),jDUT(72),gec(72)
CC  333 read(15,36,end=102) jyear,lda,jDUT(72),bz,vsw(72)
  333 read(15,*,end=102) jyear,lda,jDUT(72),bh,bz,tp,enp,vsw(72)
C       
	call blet4(jyear,CYEAR)
	jyr=jyear-jyear/100*100
	call submmdd(lda,jyr,jmn,jdy)
	jDYR(72)=jyr
	jDMN(72)=jmn
	jDDY(72)=jdy
 
C
C Check missing data:
C
	if (bz.gt.990.) then
	bz=bz_pre
	endif
	bz_pre=bz
C
	if (vsw(72).gt.9900.) then
	vsw(72)=vsw_pre
	endif
	vsw_pre=vsw(72)
C
      vbz(72)=-1.*bz*vsw(72)
C
C Move data up 1 line
        do kk=0,71
	vsw(kk)=vsw(kk+1)
	vbz(kk)=vbz(kk+1)
	jDYR(kk)=jDYR(kk+1)
	jDMN(kk)=jDMN(kk+1)
	jDDY(kk)=jDDY(kk+1)
	jDUT(kk)=jDUT(kk+1)
	  enddo
C
C
C continue read event profile:
C
  334	IF ((jDYR(24).eq.iYR).and.(jDMN(24).eq.iMN1).and.
     +(jDDY(24).eq.iDY1).and.(jDUT(24).eq.iut1)) THEN		 !==
         GOTO 335
	                else				 !==
              	GOTO 333
	               endif				 !==
C
  335	CONTINUE                     !
      do kk=0,71
	ivbz(kk)=nint(vbz(kk)/10.)
	if (ivbz(kk).gt.9999) then
	ivbz(kk)=ivbz(kk-1)+100
	endif
	if (ivbz(kk).lt.-999) ivbz(kk)=-999
	ivsw(kk)=nint(vsw(kk))
	if (ivsw(kk).eq.9999) then
	write(*,*) iyear,lda1,iev
c	pause ' '
	ivsw(kk)=ivsw(kk-1)
	endif
C	if (ivsw(kk).gt.9999) ivsw(kk)=9999
	enddo
C
   33 FORMAT(1X,3A2,1X,A2,1X,F7.2,1X,I3)
   38 format(72(1X,I4))
C Output 
 	write(*,38) (ivsw(k),k=0,71)
	write(12,38) (ivsw(k),k=0,71)
C
 	write(*,38) (ivbz(k),k=0,71)
	write(13,38) (ivbz(k),k=0,71)
  102 close(15)
  177	CONTINUE
   27 close(unit=10)
      close(unit=11)
	close(unit=12)
	close(unit=13)

C
C
   35      write(*,*) infile
      pause ' ' 
      STOP
       END

C***********************************************************
