      program extpianiaQ3NS
C.........................................................FEB. 2021
C 3h UT indices WU, WL, WE quantifying ion activity indices for North/South hemispheres
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,OUTFILEN,OUTFILES !global weighted QU,QL,QE = range index
	CHARACTER*2 AYR,AUT,AMN,ADY,AUTOUT(0:7),AMN_cur,DD2,AMN_pre
	+,ADY_pre
	CHARACTER*4 AYEAR,YEAR_pre
C	+,PYR,PMN,PDY,AYRI,AMNI
	INTEGER ivsr(0:72) !  set of 3h-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
	do iut=0,7
	iut3=iut*3
	call blet2(iut3,aut)
	autout(iut)=aut
	enddo
C
      infile='d:\web\graf\dwc\YY\MM\ut\wcMMDDUT.jYY'   ! current mn
C
      outfilen='d:\web\graf\datw\qtcYYN.txt'   !   	 WL North
	outfiles='d:\web\graf\datw\qtcYYS.txt'   !   	 WL	South
   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(73(1X,I4))
C
      write(*,*) 'ENTER Year0-start,YEAR1-end:'
	 read(*,*) jyear0,jyear1
C
      DO 330 iyear=jyear0,jyear1
	call blet4(iyear,AYEAR)
	AYR=AYEAR(3:4)
	read(AYR,*) ryr
	IYR=int(ryr)   ! 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
	infile(17:18)=ayr
	infile(36:37)=ayr
	outfilen(21:22)=ayr
	outfilea(21:22)=ayr
C
	OPEN(11,file=outfiles,access='APPEND')
	OPEN(12,file=outfilen,access='APPEND')
C
C month-to-month cycle
C
      imn2=12
      DO 51 imn=1,imn2           ! mn-to-mn
	call blet2(imn,amn)
	infile(20:21)=amn
	infile(28:29)=amn
	mndy=im(imn)
C
C  Day-to-day cycle
C
      DO 52 idy=1,mndy
	call blet2(idy,ADY)
	infile(30:31)=ady
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
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		  South hemisphere
	sumn1=0.		 ! count of Wf <-1
	sump1=0.		 ! count of Wf >+1
	icntp1=0
	swegp1=0.
      icntn1=0
	swegn1=0.
C
C		  North hemisphere
	sumn2=0.		 ! count of Wf <-1
	sump2=0.		 ! count of Wf >+1
	icntp2=0
	swegp2=0.
      icntn2=0
	swegn2=0.
C
C 3h UT period:
C
	DO 153 iut3=iut1,iut1+2
	call blet2(iut3,AUT)
      infile(32:33)=aut
C
   	OPEN(10,file=infile,action='READ')
C
C/ cycle for lats:
       glat=-90.0	                  ! AD
Crem	  do 54 ilat=1,71
C
	  do 54 ilat=1,36       ! South hemisphere
	glat=glat+2.5				  ! AD
	elat=1.0+0.0001235*glat**2	  ! AD
	  read(10,32) (ivsr(k),k=0,72) 
	if (ilat.eq.36) exit  !===============> exclude glat=0
C
	DO 333 ilon=0,72
          if (ivsr(ilon).gt.1) then
        iad=ivsr(ilon)
c	sump=sump+weg(iad)*ivsr(ilon)
	sump1=sump1+elat*weg(iad)*ivsr(ilon)  ! AD
	icntp1=icntp1+1
	swegp1=swegp1+weg(iad)
	goto 333
	endif
          if (ivsr(ilon).lt.-1) then
        iad=-ivsr(ilon)
c	sumn=sumn+weg(iad)*ivsr(ilon)
	sumn1=sumn1+elat*weg(iad)*ivsr(ilon)   ! AD
	icntn1=icntn1+1
	swegn1=swegn1+weg(iad)
	goto 333
	endif
  333	CONTINUE
   54 CONTINUE
C
	  do 154 ilat=37,71       ! North hemisphere
	glat=glat+2.5				  ! AD
	elat=1.0+0.0001235*glat**2	  ! AD
	  read(10,32) (ivsr(k),k=0,72) 
C
	DO 1333 ilon=0,72
          if (ivsr(ilon).gt.1) then
        iad=ivsr(ilon)
c	sump=sump+weg(iad)*ivsr(ilon)
	sump2=sump2+elat*weg(iad)*ivsr(ilon)  ! AD
	icntp2=icntp2+1
	swegp2=swegp2+weg(iad)
	goto 1333
	endif
          if (ivsr(ilon).lt.-1) then
        iad=-ivsr(ilon)
c	sumn=sumn+weg(iad)*ivsr(ilon)
	sumn2=sumn2+elat*weg(iad)*ivsr(ilon)   ! AD
	icntn2=icntn2+1
	swegn2=swegn2+weg(iad)
	goto 1333
	endif
 1333	CONTINUE
  154 CONTINUE
	close(unit=10)
  153 CONTINUE
      xpos1=0.
      if (icntp1.gt.0.) then
C		xpos=sump/swegp		 ! /sum(weights)
		xpos1=sump1/icntp1		 ! /cnt(Wf+)
	endif
      xneg1=0.
      if (icntn1.gt.0.) then
C			xneg=sumn/swegn	  ! /sum(weights)
			xneg1=sumn1/icntn1	  ! /cnt(Wf-)
	endif
	xrange1=xpos1-xneg1
	wdyp1=wdyp1+xpos1
	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
      close(unit=11)
  330	CONTINUE
      pause ' '
	STOP
	END
