      program extpianias3x
C.........................................................June 2021
C Proxy Wx indices
C.........................................................Apr. 2020
C Shub. GTEC indices
C.........................................................Dec. 2020
C 3h UT indices QU, QL, QE quantifying ion activity indices 
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),xwneg(0:23),xwpos(0:23),xwran(0:23)  ! 1h UT 
     +,ivsr(0:72)  !  set of 1h-Wx indices for long=-180:5:180 at given lati
	CHARACTER*128 INFILE,OUTFILE !global weighted WUx,WLx,WEx = range index
	CHARACTER*2 AYR,AUT,AMN,ADY
	CHARACTER*4 AYEAR
c	INTEGER ivsu(0:24),ivsl(0:24),ivse(0:24) !  set of 1h-W indices for long=-180:5:180 at given lati
		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
C
      infile='d:\web\graf\duc\YY\MM\ut\wxMMDDUT.sYY'   ! current mn
C
      outfile='d:\web\graf\duc\wuleYYMN.txt'   !  Cnt(W+), W-, Wrange
   17 format(A2)
c-   31 FORMAT('YYMMDD',1X,8(2X,A2,1X),8(3X,A2,1X),8(2X,A2,1X))
   32 FORMAT(73(1X,I4))
C
   1        write(*,*) ' Enter iyear <-1000> STOP: '
	read(*,*) iyear
      if (iyear.lt.0) GOTO 2
        write(*,*) ' Enter imn: '
	read(*,*) imn
        write(*,*) ' Enter idy1,idy2: '
	read(*,*) idy1,idy2

  330	CONTINUE
	call blet4(iyear,AYEAR)
	IYR=iyear-iyear/100*100   ! current yr
	z1=iyear/4.0
      jz=int(z1)*4

      IF(jz.EQ.iyear) THEN
               IM(2)=29
	dnr=366.
        ELSE
                IM(2)=28
	  	dnr=365.
	       ENDIF
      AYR=AYEAR(3:4)  ! for current day  	
	infile(17:18)=ayr
	infile(36:37)=ayr
	outfile(21:22)=ayr
C
C month
C	 
	call blet2(imn,AMN)
	infile(20:21)=AMN
	infile(28:29)=AMN
	mndy=im(imn)
	call blet2(imn,amn)
	outfile(23:24)=amn
	OPEN(11,file=outfile)
	mndy=im(imn)
C
C  Day-to-day cycle
C
      DO 52 idy=idy1,idy2		 !DAYYYYYYYYYYYYYY
	call blet2(idy,ADY)
	infile(30:31)=ady
C
C  UT cycle
C
      DO k=0,23
	xwneg(k)=0.
	xwpos(k)=0.
	xwran(k)=0.
	ENDDO
	 wdyp=0.
	wdyn=0.
	wdyr=0.
	DO 53 iut=0,23	  ! UTTTTTTTTTTTTTT
C
	sumn=0.		 ! count of Wf <-1
	sump=0.		 ! count of Wf >+1
	icntp=0
	swegp=0.
      icntn=0
	swegn=0.
C
	call blet2(iut,AUT)
      infile(32:33)=aut
C
   	OPEN(10,file=infile,action='READ')
C
C/ cycle for lats:
       glat=-90.0	                  ! AD
	  do 54 ilat=1,71
	glat=glat+2.5				  ! AD
	elat=1.0+0.0001235*glat**2	  ! AD
	  read(10,32,end=51) (ivsr(k),k=0,72) 
C
	DO 333 ilon=0,72
          if (ivsr(ilon).gt.1) then
        iad=ivsr(ilon)
c	sump=sump+weg(iad)*ivsr(ilon)
	sump=sump+elat*weg(iad)*ivsr(ilon)  ! AD
	icntp=icntp+1
	swegp=swegp+weg(iad)
	goto 333
	endif
          if (ivsr(ilon).lt.-1) then
        iad=-ivsr(ilon)
c	sumn=sumn+weg(iad)*ivsr(ilon)
	sumn=sumn+elat*weg(iad)*ivsr(ilon)   ! AD
	icntn=icntn+1
	swegn=swegn+weg(iad)
	goto 333
	endif
  333	CONTINUE
   54 CONTINUE
	close(unit=10)
  153 CONTINUE
      xpos=0.
      if (icntp.gt.0.) then
C		xpos=sump/swegp		 ! /sum(weights)
		xpos=sump/icntp		 ! /cnt(Wf+)
	endif
      xneg=0.
      if (icntn.gt.0.) then
C			xneg=sumn/swegn	  ! /sum(weights)
			xneg=sumn/icntn	  ! /cnt(Wf-)
	endif
	xrange=xpos-xneg
	wdyp=wdyp+xpos
	wdyn=wdyn+xneg
      wdyr=wdyr+xrange
      xwneg(iut)=xneg
	xwpos(iut)=xpos
	xwran(iut)=xrange
	write(*,33) AYR,AMN,ADY,AUT,xwpos(iut),xwneg(iut)
	+,xwran(iut)
	write(11,33) AYR,AMN,ADY,AUT,xwpos(iut),xwneg(iut)
	+,xwran(iut)
   51	close(unit=10) 
   53	CONTINUE           ! end UT cycle
c-       iwpos(24)=nint(wdyp/24.*100.)  ! Daily mean not to output
c-	iwneg(24)=nint(wdyn/24.*100.)	   ! Daily mean not to output
c-	iwran(24)=nint(wdyr/24.*100.)	   ! Daily mean not to output
   33 format(1X,A4,2A2,1X,A2,3(1X,F6.2))
   52 CONTINUE  ! daily cycle
      close(unit=11)
	GOTO 1
    2  pause ' '
	STOP
	END
