      program extpianiaD3st
C.........................................................Jan. 2021
C 3h UT indices DU, DL, DE quantifying ion activity indices for station
C
C Include weight e(glat)
C
C.........................................................Nov. 2020
C Weighted PIA -pos ion activity, NIA - neg ion activity
C extracted from global Wc(TEC) index derived from JPL GIM-TEC (IONEX)
C for ADY_pre and forecast for DD2
C
	DIMENSION weg(4),IM(12),iwneg(0:8),iwpos(0:8),iwran(0:8)  ! 3h UT 
	CHARACTER*128 INFILE,OUTFILE !local weighted DU,DL,DE = range index
	CHARACTER*2 AYR,AUT,AMN,ADY,AUTOUT(0:7),st,CYR,CMN,CDY
	CHARACTER*4 AYEAR
C	+,PYR,PMN,PDY,AYRI,AMNI
	INTEGER iw(0:23) !  set of 1h-W indices for 0:23 UT
		DATA WEG/0.,0.6931,1.0986,1.3863/  ! =ln1, ln2, ln3, ln4
            DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
C
       write(*,*) 'ENTER st:'
	 read(*,*) st
       write(*,*) 'ENTER glat:'
	 read(*,*) glat
	elat=1.0+0.0001235*glat**2	  ! AD
C
      write(*,*) 'ENTER Year0-start,YEAR1-end:'
	 read(*,*) jyear0,jyear1
C
	infile='d:\data\IWG\st\stYRwf.txt'
	infile(13:14)=st
	infile(16:17)=st
C
	do iut=0,7
	iut3=iut*3
	call blet2(iut3,aut)
	autout(iut)=aut
	enddo
C
       DO 200 jyear=jyear0,jyear1
	call blet4(jyear,AYEAR)
      AYR=AYEAR(3:4)  ! 
		infile(18:19)=AYR
	OPEN(10,file=infile,action='READ')
	outfile=infile
	outfile(21:21)='3'
	OPEN(11,file=outfile)
C
	jyr=jyear-(jyear/100)*100 
C
	yrp=float(jyr_pre)
	if(int(yrp/4.)*4.eq.jyr) THEN 
	im(2)=29
	ndnr1=366
	                          ELSE
	im(2)=28
	ndnr1=365
	ENDIF
C
C
   17 format(A2)
C   31 FORMAT('YYMMDD',1X,24(2X,A2,1X),2X,'DAY')
   31 FORMAT('YYMMDD',1X,8(2X,A2,1X),8(3X,A2,1X),8(2X,A2,1X))
   32 FORMAT(3A2,1X,24(I4,1X))
C
  330	CONTINUE
C
C month-to-month cycle
C
      DO 51 imn=1,12           ! mn-to-mn
Crem      imn=imn_pre
	call blet2(imn,amn)
	mndy=im(imn)
C
C  Day-to-day cycle
C
      DO 52 idy=1,mndy
Crem       idy=idy_pre
	call blet2(idy,ADY)
Crem      ADY=ADY_pre
C
	IF ((AMN.eq.'01').and.(ADY.eq.'01')) THEN
	write(*,31) (autout(k),k=0,7),(autout(k),k=0,7),(autout(k),k=0,7)
	write(11,31) (autout(k),k=0,7),(autout(k),k=0,7),(autout(k),k=0,7)
	ENDIF
		read(10,32,end=108) CYR,CMN,CDY,(iw(k),k=0,23)		! daily/UT W-index

C  UT cycle
C
      DO k=0,7
	iwneg(k)=0
	iwpos(k)=0
	iwran(k)=0
	ENDDO
	 wdyp=0.
	wdyn=0.
	wdyr=0.
	DO 53 iut=0,7
	iut1=iut*3
C
	sumn=0.		 ! count of Wf <-1
	sump=0.		 ! count of Wf >+1
	icntp=0
	swegp=0.
      icntn=0
	swegn=0.
C
C 3h UT period:
C
	DO 153 iut3=iut1,iut1+2
	call blet2(iut3,AUT)
C
C
          if (iw(iut3).gt.1) then
        iad=iw(iut3)
	sump=sump+elat*weg(iad)*iw(iut3)  ! AD
	icntp=icntp+1
	swegp=swegp+weg(iad)
	goto 153
	endif
          if (iw(iut3).lt.-1) then
        iad=-iw(iut3)
	sumn=sumn+elat*weg(iad)*iw(iut3)   ! AD
	icntn=icntn+1
	swegn=swegn+weg(iad)
	goto 153
	endif
  153 CONTINUE
      xpos=0.
      if (icntp.gt.0.) then
		xpos=sump/icntp		 ! /cnt(Wf+)
	endif
      xneg=0.
      if (icntn.gt.0.) then
			xneg=sumn/icntn	  ! /cnt(Wf-)
	endif
	xrange=xpos-xneg
	wdyp=wdyp+xpos
	wdyn=wdyn+xneg
      wdyr=wdyr+xrange
      iwneg(iut)=nint(xneg*100.)
	iwpos(iut)=nint(xpos*100.)
	iwran(iut)=nint(xrange*100.)
   53	CONTINUE           ! end UT cycle
       iwpos(8)=nint(wdyp/8.*100.)  ! Daily mean not to output
	iwneg(8)=nint(wdyn/8.*100.)	   ! Daily mean not to output
	iwran(8)=nint(wdyr/8.*100.)	   ! Daily mean not to output
	write(*,33) AYR,AMN,ADY,(iwpos(k),k=0,7),(iwneg(k),k=0,7)
	+,(iwran(k),k=0,7)
	write(11,33) AYR,AMN,ADY,(iwpos(k),k=0,7),(iwneg(k),k=0,7)
	+,(iwran(k),k=0,7)
   33 format(3A2,8(1X,I4),8(1X,I5),8(1X,I4))
   52 CONTINUE  ! daily cycle
   51 CONTINUE  ! Monthly cycle
 200  CONTINUE  ! Annual cycle
  108	close(unit=10)
      close(unit=11)
      pause ' '
	STOP
	END
