      program cntpianiaD3st
C.........................................................Sep. 2022
C 3h UT indices WU, WL, WE quantifying TEC activity indices for station
C TEMP: input IONOLAB_TEC: SHUB_med; 
C TEMP: output: Windex, WU-3h, WL-3h, WE-3h
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),iwneg(0:8),iwpos(0:8),iwran(0:8)  ! 3h UT 
	+,glats(3),glons(3),teci(0:23),tmed(0:23)
	CHARACTER*128 INFILE1,INFILE2,OUTFILE,txt !local weighted WU,WL,WE = 3h range index
	CHARACTER*2 AYR,AUT,AMN,ADY,AUTOUT(0:7)
	CHARACTER*4 AYEAR,stas(3),stat
	+,PYR,PMN,PDY
        CHARACTER*8 YYYYMNDD
	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 STAS/'dyng','onsa','mobn'/
	DATA glats/38.08,57.4,55.1/
	DATA glons/23.93,11.93,36.6/

C
      infile1='multiple.txt'
C
	infile2='GTEC.out'
	OPEN(12,file=infile2,action='READ')
	read(12,*) txt
C
	outfile='D3statind.txt'
	OPEN(10,file=outfile)
   30	format(1X,A4,2(1X,F6.2))
C
C Cycle on stations
C		   +++++++++++++++++++++++++++++++++++++
         DO 123 ista=1,3
	stat=stas(ista)
	glat=glats(ista)
	glon=glons(ista)
	write(*,30) stat,glat,glon
	write(10,30) stat,glat,glon

	elat=1.0+0.0001235*glat**2	  ! AD
C
C
	do iut=0,7
	iut3=iut*3
	call blet2(iut3,aut)
	autout(iut)=aut
	enddo
C
        do k=0,23
	teci(k)=0.0
	iw(k)=0
	tmed(k)=0.0
C
C Read Shub median:
      read(12,*, end=52) kldy,klmn,klyear,rlt,kdy,kmn,kyear,rut
	+,glat,glon,ssn,tmed(k)
	  enddo
C
C
C Read ionolab data
	OPEN(11,file=infile1,action='READ')
	read(11,*) txt
C
        DO 112 kk=1,3  ! 3 days
         do 111 jut=0,23
	tcut=0.
	icnt=0
C test
  110    read(11,*,end=330) YYYYMNDD,IUTM,tc1,tc2,tc3
	if (ista.eq.1) tcj=tc1
      if (ista.eq.2) tcj=tc2
	if (ista.eq.3) tcj=tc3
	AYEAR= YYYYMNDD(1:4)
	AYR=AYEAR(3:4)
	read(AYEAR,*) ryear
	jyear=int(ryear)
	jyr=jyear-(jyear/100)*100 
	AMN= YYYYMNDD(5:6)
      ADY= YYYYMNDD(7:8)
	iutm=iutm/100 ! exclude seconds
	iut=iutm/100
	imin=iutm-iut*100
	if (iut.eq.23) then
	PYR=AYR
	PMN=AMN
	PDY=ADY
	endif
C
	IF (iut.eq.jut) THEN
	 tcut=tcut+tcj
	icnt=icnt+1
	GOTO 110
	                 ELSE
	backspace(11)
	tcut=tcut/icnt
	teci(jut)=tcut
	delog=0
	xpar=tcut
	xmed=tmed(jut)
      call windex(xpar,xmed,delog,iwind)
	iw(jut)=iwind
	    if(jut.lt.24) then
	GOTO 111
	      else
	GOTO 330
	       endif
	                 ENDIF
  111	CONTINUE
   17 format(A2)
C   31 FORMAT('YYMMDD',1X,24(2X,A2,1X),2X,'DAY')
   31 FORMAT(1X,3(A2),1X,I2,2(1X,F7.2),1X,I4)
   32 FORMAT(3A2,1X,I2,1X,24(I4,1X))
C
  330	CONTINUE
C
C output W index:
      do m=0,23
	mut=m+(ista-1)*24
      write(*,31) PYR,PMN,PDY,mut,teci(m),tmed(m),iw(m)
	write(10,31) PYR,PMN,PDY,mut,teci(m),tmed(m),iw(m)
	enddo
C/		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(10,33) AYR,AMN,ADY,(iwpos(k),k=0,7),(iwneg(k),k=0,7)
	+,(iwran(k),k=0,7)
  112	CONTINUE
   51 close(11)
 123	 CONTINUE
   33 format(3A2,8(1X,I4),8(1X,I5),8(1X,I4))
   52   close(unit=12)
  108	close(unit=10)
      pause ' '
	STOP
	END
