      program extpianiaWCmm
C .........................................................June 2021
C Add weights ro=log((np+1)/(nn+1)) to produce dominant DU, DL index instead of WU, WL
C One output file providing DU or DL 
C No output of DE = DU-DE indices< because DE=DU or DE=DL
C..........................................................May 2021
C Output if date + 4 lines for each day 
C..........................................................Apr.2021
C 1h UT indices WU, WL, WE quantifying ion activity indices 
C JPL IONEX GIM-TEC extracted from JPLI for 00, 15, 30, 45 min.
C.........................................................Dec. 2020
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),iwnp(0:96)
C-	+,iwpos(0:96)
C-	+,iwran(0:96)
	CHARACTER*128 INFILE,OUTFILE
C	+,indate !global weighted WN,WL,WU = range index
	CHARACTER*2 AYR,AUT,AMN,ADY,AMM
C+	+,AMN_cur,DD2,AMN_pre
C+	+,ADY_pre
	CHARACTER*4 YEAR,AUTOUT(0:95),TTMM
C	+,PYR,PMN,PDY,AYRI,AMNI,YEAR_pre
	INTEGER ivsr(0:72) !  set of Wf 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+	indate='d:/web/graf/date'    !!! Tamara !!!!!!!	 
C+	open(10,file=indate)							 
C+	read(10,799) YEAR,AMN_cur,DD2,YEAR_pre,AMN_pre,ADY_pre,lda,lda_pre  !FORMAT NUMBER = 799
C+	WRITE(*,798) YEAR,AMN_cur,DD2,YEAR_pre,AMN_pre,ADY_pre,lda,lda_pre !FORMAT NUMBER = 798
C+	close(unit=10)
C+  798	format(' YEAR, AMN, DD2 = ',A4,2(1X,A2),1X,A4,2(1X,A2),2(1X,I3))
C+  799	format(A4,2(1X,A2),1X,A4,2(1X,A2),2(1X,I3))	
C
	nn=0
	do iut=0,23
	call blet2(iut,aut)
	TTMM(1:2)=AUT
	do mm=0,45,15
	call blet2(mm,amm)
	TTMM(3:4)=AMM
	autout(nn)=TTMM
	nn=nn+1
	enddo
	enddo
C      infile='d:\web\dlg\YY\MM\wcMMDDUTMM.jYY'   ! current mn	  Observed 0,1`5,30,45min
C      infile='d:\web\gtc\YY\MM\wsMMDDUTMM.sYY'   ! current mn Shub. W(GTEC) index
      infile='d:\web\graf\dwu\YY\MM\ut\wcMMDDUTMM.uYY'   ! uqrg_15min source
C
C      outfilen='d:\web\graf\datw\wlnYYmm.txt'   !  Cnt(W-) 	 WL	Observed 0,1`5,30,45min
C      outfilep='d:\web\graf\datw\wupYYmm.txt'   !  Cnt(W+) 	 WU	Observed 0,1`5,30,45min
C      outfiler='d:\web\graf\datw\werYYmm.txt'   !  Cnt(Wrange) WE	Observed 0,1`5,30,45min
c      outfilen='d:\web\graf\dats\wlnYYms4.txt'   !  Cnt(W-) 	 WL	Observed 0,1`5,30,45min
c      outfilep='d:\web\graf\dats\wupYYms4.txt'   !  Cnt(W+) 	 WU	Observed 0,1`5,30,45min
c      outfiler='d:\web\graf\dats\werYYms4.txt'   !  Cnt(Wrange) WE	Observed 0,1`5,30,45min
      outfile='d:\web\graf\datu\dulYYmm.txt'   !  Cnt(W-) 	DU&DL	Observed 0,1`5,30,45min
C
   17 format(A2)
   31 FORMAT('YYMMDD',96(1X,A4),3X,'DAY')
   32 FORMAT(73(1X,I4))
C
C Yearly cycle
C
C      DO 500 iyear=1996,2021  ! 
CREM      DO 500 iyear=1997,2021  ! 
C      iyear=2021  ! Temp
      iyear=2014
	call blet4(iyear,YEAR)
C
      AYR=YEAR(3:4)  ! for current day  	
C+	IYR=int(ryr)   ! current yr
      iyr=iyear-iyear/100*100
	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(38:39)=AYR
	outfile(21:22)=AYR
C
	OPEN(11,file=outfile,access='APPEND')
C
C month-to-month cycle
C
       imn1=1
c-
c-       imn1=5	   ! Temp 2021
c-	imn2=12
	imn2=12
	if (iyear.eq.1996) imn1=12
	if (iyear.eq.2021) imn2=8
C
	DO 51 imn=imn1,imn2
	call blet2(imn,AMN)
	mndy=im(imn)
	infile(20:21)=amn
	infile(28:29)=amn
C
C  Day-to-day cycle
C
      DO 52 idy=1,mndy
       call blet2(idy,ADY)
	infile(30:31)=ady
C
C?  330	CONTINUE
C
	IF ((AMN.eq.'01').and.(ADY.eq.'01')) THEN
	write(*,31) (autout(k),k=0,95)
C>>	write(11,31) (autout(k),k=0,95)
C>>	write(12,31) (autout(k),k=0,95)
C>>	write(13,31) (autout(k),k=0,95)
	ENDIF
C
C  UT cycle
C
      DO k=0,96
	iwnp(k)=0
	ENDDO
	 wdyp=0.
	wdyn=0.
	wdyr=0.
C
C
      iutmm=0
	DO 53 iut=0,23
	call blet2(iut,AUT)
      infile(32:33)=aut
C
C Minute cycle
C
      DO 200 imm=0,45,15
	call blet2(imm,AMM)
	infile(34:35)=amm
   	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:
       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=153) (ivsr(k),k=0,72) 
C
	DO 333 ilon=0,72
          if (ivsr(ilon).gt.1) then
        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)
	sumn=sumn+elat*weg(iad)*ivsr(ilon)   ! AD
	icntn=icntn+1
	swegn=swegn+weg(iad)
	goto 333
	endif
  333	CONTINUE		! end long cycle
   54 CONTINUE        ! end lat cycle
  153      xpos=0.
      ro=alog10(icntp+1.0)-alog10(icntn+1.0)
C_pre      if (icntp.gt.0.) then
             if ((ro.ge.0.).and.(icntp.gt.0.)) then
C_pre		xpos=sump/icntp		 ! /cnt(Wf+)
		xpos=sump/icntp*ro	 ! /cnt(Wf+)
		iwnp(iutmm)=nint(xpos*100.)
	endif
      xneg=0.
C_pre      if (icntn.gt.0.) then
      if ((ro.le.0.).and.(icntn.gt.0.)) then
X_pre			xneg=sumn/icntn	  ! /cnt(Wf-)
		xneg=sumn/icntn*abs(ro)	  ! /cnt(Wf-)
	      iwnp(iutmm)=nint(xneg*100.)
	endif
C-	xrange1=xpos-xneg
C-	xrange2=xneg-xpos
C-	rat=float(icntp)/float(icntn)
C-	if (rat.ge.1.0) then
C-	xrange=xrange1*(rat-1.0)
C-	else
C-	xrange=xrange2*(rat+1.0)
C-	endif
	wdyp=wdyp+xpos
	wdyn=wdyn+xneg
C-    wdyr=wdyr+xrange
C?      iwnp(iutmm)=nint(xneg*100.)
C?	iwpos(iutmm)=nint(xpos*100.)
C-	iwran(iutmm)=nint(xrange*100.)
	iutmm=iutmm+1
  	close(unit=10)
  200	continue  !end minute cycle
   53	CONTINUE           ! end UT cycle
       sum=0.
       do k=0,95
	sum=sum+iwnp(k)
	 enddo
	iwnp(96)=nint(sum/96.)
C?       iwpos(96)=nint(wdyp/96.*100.)
C?	iwneg(96)=nint(wdyn/96.*100.)
C-	iwran(96)=nint(wdyr/96.*100.)
C-	write(*,33) AYR,AMN,ADY,(iwpos(k),k=0,96)
C-	write(11,33) AYR,AMN,ADY,(iwpos(k),k=0,96)
C-	write(*,33) YEAR,AMN,ADY,(iwnp(k),k=0,96)
		write(11,33) YEAR,AMN,ADY,(iwnp(k),k=0,96)
C	write(11,34) YEAR,AMN,ADY,iwpos(96)
C	write(11,33) (iwpos(k),k=0,23)
C	write(11,33) (iwpos(k),k=24,47)
C	write(11,33) (iwpos(k),k=48,71)
C	write(11,33) (iwpos(k),k=72,95)
C-	write(*,33) YEAR,AMN,ADY,(iwneg(k),k=0,96)
C-	write(12,33) YEAR,AMN,ADY,(iwneg(k),k=0,96)
c	write(12,34) YEAR,AMN,ADY,iwneg(96)
c	write(12,33) (iwneg(k),k=0,23)
c	write(12,33) (iwneg(k),k=24,47)
c	write(12,33) (iwneg(k),k=48,71)
c	write(12,33) (iwneg(k),k=72,95)
C-	write(*,33) YEAR,AMN,ADY,(iwran(k),k=0,96)
C-	write(13,33) YEAR,AMN,ADY,(iwran(k),k=0,96)
c	write(13,34) YEAR,AMN,ADY,iwran(96)
c	write(13,33) (iwran(k),k=0,23)
c	write(13,33) (iwran(k),k=24,47)
c	write(13,33) (iwran(k),k=48,71)
c	write(13,33) (iwran(k),k=72,95)
	write(*,*) YEAR,AMN,ADY
   52 CONTINUE  ! daily cycle
   51 CONTINUE  ! Monthly cycle
      close(unit=11)
C-	close(unit=12)
C-	close(unit=13)
  500	CONTINUE  ! Yearly cycle
   34 format(A4,2(1X,A2),1X,I4)
   33 format(A4,2(A2),96(1X,I4),2X,I4)
C   33 format(24(1X,I4))
      pause ' '
	STOP
	END
