      program extpianiaWf
C.........................................................Nov. 2020
C Weighted PIA -pos ion activity, NIA - neg ion activity
C extracted from global Wf index derived from GIM-fcF2 (IRI-Plas-TEC conversion)
C
	DIMENSION weg(4),IM(12),iwneg(0:24),iwpos(0:24),iwran(0:24)
	CHARACTER*48 INFILE,OUTFILEN,OUTFILEP,OUTFILER !global weighted Wfp,Wfn,Wfr = range index
	CHARACTER*2 AYR,AUT,AMN,ADY,AUTOUT(0:23)
C	+,PYR,PMN,PDY,AYRI,AMNI
	INTEGER ivsr(0:72) !  set of Wf indices for long=-10: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
      infile='d:\web\dwf\YY\MM\wfMMDDUT.jYY'   ! current mn
	outfilen='wfnYR.txt'   ! Cnt(Wf-) => /Sum_pre(weights)
	outfilep='wfpYR.txt'   ! Cnt(Wf+) => /Sum_pre(weights)
	outfiler='wfrYR.txt'
	do iut=0,23
	call blet2(iut,aut)
	autout(iut)=aut
	enddo
   17 format(A2)
   31 FORMAT('YYMMDD',1X,24(2X,A2,1X),2X,'DAY')
   32 FORMAT(73(1X,I4))
C
  330	WRITE(*,*)' ENTER YR OF INPUT FILE=	'
	read(*,17) AYR
	read(AYR,*) ryr
	IYR=int(ryr)   ! current yr
	if (iyr.lt.90) then
	iyyyy=2000+iyr
	              else
	iyyyy=1900+iyr
	endif
	z1=iyyyy/4.0
      jz=int(z1)*4

      IF(jz.EQ.iyyyy) THEN
               IM(2)=29
	dnr=366.
        ELSE
                IM(2)=28
	  	dnr=365.
	       ENDIF
	infile(12:13)=ayr
	infile(28:29)=ayr
	outfilep(4:5)=ayr
	outfilen(4:5)=ayr
	outfiler(4:5)=ayr
C
	OPEN(11,file=outfilep,access='APPEND')
	OPEN(12,file=outfilen,access='APPEND')
	OPEN(13,file=outfiler,access='APPEND')

	write(*,31) (autout(k),k=0,23)
	write(11,31) (autout(k),k=0,23)
	write(12,31) (autout(k),k=0,23)
	write(13,31) (autout(k),k=0,23)
	 
C
C month-to-month cycle
C
      imn2=12
	if (iyr.eq.20) imn=10  ! TEMP
      DO 51 imn=1,imn2
	call blet2(imn,amn)
	infile(15:16)=amn
	infile(20:21)=amn
	mndy=im(imn)
C
C  Day-to-day cycle
C
      DO 52 idy=1,mndy
	call blet2(idy,ADY)
	infile(22:23)=ady
C
C  UT cycle
C
      DO k=0,24
	iwneg(k)=0
	iwpos(k)=0
	iwran(k)=0
	ENDDO
	 wdyp=0.
	wdyn=0.
	wdyr=0.
	DO 53 iut=0,23
	call blet2(iut,AUT)
      infile(24:25)=aut
C
   	OPEN(10,file=infile,action='READ')
	sumn=0.		 ! count of Wf >=1
	sump=0.		 ! count of Wf <=-1
	icntp=0
	swegp=0.
      icntn=0
	swegn=0.
C
C/ cycle for lats:
	  do 54 ilat=1,71
	  read(10,32) (ivsr(k),k=0,72) 
C
	DO 333 ilon=0,72
          if (ivsr(ilon).gt.1) then
        iad=ivsr(ilon)
	sump=sump+weg(iad)*ivsr(ilon)
	icntp=icntp+1
	swegp=swegp+weg(iad)
	goto 333
	endif
          if (ivsr(ilon).lt.-1) then
        iad=-ivsr(ilon)
	sumn=sumn+weg(iad)*ivsr(ilon)
	icntn=icntn+1
	swegn=swegn+weg(iad)
	goto 333
	endif
  333	CONTINUE
   54 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
      iwneg(iut)=nint(xneg*100.)
	iwpos(iut)=nint(xpos*100.)
	iwran(iut)=nint(xrange*100.)
	close(unit=10)
   53	CONTINUE           ! end UT cycle
       iwpos(24)=nint(wdyp/24.*100.)
	iwneg(24)=nint(wdyn/24.*100.)
	iwran(24)=nint(wdyr/24.*100.)
	write(*,33) AYR,AMN,ADY,(iwpos(k),k=0,24)
	write(11,33) AYR,AMN,ADY,(iwpos(k),k=0,24)
	write(*,33) AYR,AMN,ADY,(iwneg(k),k=0,24)
	write(12,33) AYR,AMN,ADY,(iwneg(k),k=0,24)
	write(*,33) AYR,AMN,ADY,(iwran(k),k=0,24)
	write(13,33) AYR,AMN,ADY,(iwran(k),k=0,24)
   33 format(3A2,24(1X,I4),2X,I4)
   52 CONTINUE  ! daily cycle
   51 CONTINUE  ! Monthly cycle
      close(unit=11)
	close(unit=12)
	close(unit=13)
      if (iyr.ne.20) GOTO 330
      pause ' '
	STOP
	END
