	subroutine subwp3nc(AYR,AMN,IDY1,IDY2)
C.....................................................Dec  2023
C
C.....................................................Nov. 2012
C JPL results
C.....................................................Oct.2012
C Update monthly Wp, Np, Pp files
C
C.....................................................Apr. 2012
C
C IONEX maps
C for glat=-87.5:2.5:87.5, glong=-180:5:175  (5112)
C Correction of (W+max-W-max)*(1+cntstrm/5112)
C.....................................................Feb. 2012
C
C.....................................................Aug.2011
C
C........................................................Apr.2008
C...... Count RANGE only when W+=>3, and/or W-<=-3 ..............
C
c........................................................Oct.2007
C cntstrm=total number of W+=3,W+=4, W-=-3,W-=-4 for the IONEX-map
C
c .......................................................Sep. 2007
C Range of TEC indices from 35 latitude results TEC.fYR.........
C
c from hourly IONEX maps
c
	DIMENSION RESP(0:71),RESN(0:71),DHR(0:71),IM(12)
	+,cnts(71,3,0:23),xind(71,3,0:23)
     +,RES(0:71),RUP(0:71),RLO(0:71)
	+,clat(71,3,0:23),relat(71,3,0:23)
	CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
	CHARACTER*64,infile,outfile,outfilep,outfilen,outinfo,title
	CHARACTER*135 wpast(0:365),wwpos(0:365),wwneg(0:365),txt
	CHARACTER*2 AYR,AMN,AUT,ADY
	+,aa(24),bb
	CHARACTER*6 YRMNDY
	CHARACTER*4 AYEAR,AMEAN
      DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
C Controlling Selected Latitudes 1=Yes, 0 = No:
	DATA aa/'00','01','02','03','04','05','06','07','08','09','10'
	+,'11','12','13','14','15','16','17','18','19','20','21','22','23'/
	bb='UT'
	amean='mean'
	title='YYMMDD UT CNTNEG CNTPOS CNTTOT   TOT%    Wp'
	CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT

C start
	AYEAR(1:2)='20'
	AYEAR(3:4)=ayr
	read(ayr,*) year
	iyr=int(year)

	read(amn,*) rmn
	mn=int(rmn)
	if(int(year/4.)*4.eq.iyr) THEN 
	ndnr=366
	im(2)=29
	                          ELSE
	ndnr=365
	im(2)=28
	ENDIF

 13   FORMAT(I2)

CC CALCULATION OF DAYS OF YEAR	
      infile='c:\web\graf\dwc\YR\MN\ut\wcMNDYUT.jYR' !W(glat & glong map)!Tamara
  	infile(17:18)=ayr									   !Tamara
	infile(20:21)=amn									   !Tamara
  	infile(36:37)=ayr									   !Tamara
	infile(28:29)=amn									   !Tamara
C      outfile='c:\web\graf\datw\wptecYR.txt'	!W-index from W-maps geogr.coord.!Tamara
	outfile='c:\web\graf\datw\wptecYR.txt'	!              !Tamara
	outfile(23:24)=AYR									   !Tamara
 	outfilep=outfile
	outfilen=outfile
	outfilep(18:18)='p'									   !Tamara
	outfilen(18:18)='n'									   !Tamara
	outinfo='wjinfoYR.txt'
	outinfo(7:8)=ayr									   !Tamara
	OPEN(10,FILE=OUTINFO,ACCESS='APPEND')
	write(10,14) title
   14	format(1X,A43)
C
  	OPEN(12,FILE=OUTFILE)
	OPEN(13,FILE=OUTFILEP)
	OPEN(14,FILE=OUTFILEN)
 198	format(4X,A2,24(3X,A2),3X,A5)
C
	lda1=ndoy(iyr,mn,idy1)  ! day-of-year
	mndy=im(mn)
  130	format(A135)

	if((mn.eq.1).and.(idy1.eq.1)) then 

	write(12,198) bb,aa,amean
	write(13,198) bb,aa,amean
	 write(14,198) bb,aa,amean
      write(*,198) bb,aa,amean
	goto 197 ! avoid read-re-write
	endif
C
C read-re-write outfile for days from 1 to lda1-1:
	ncnt=0
	do 149 ipast=0,lda1-1
	read(12,130,END=151,ERR=2) txt
	wpast(ipast)=txt
	YRMNDY=txt(1:6)
  151	read(13,130,END=152,ERR=2) txt
 	wwpos(ipast)=txt
  152	read(14,130,END=153,ERR=2) txt
 	wwneg(ipast)=txt
	if (yrmndy.eq.'      ') goto 149
	ncnt=ncnt+1
  149	continue
  153 close(unit=12)
      close(unit=13)
	close(unit=14)
C
  	OPEN(12,FILE=OUTFILE)
	OPEN(13,FILE=OUTFILEP)
	OPEN(14,FILE=OUTFILEN)
C re-write past data:
	do ipast=0,ncnt-1
	write(12,130) wpast(ipast)
cc	write(*,130) wpast(ipast)  ! test
	write(13,130) wwpos(ipast)
	write(14,130) wwneg(ipast)
	enddo
      close(unit=12)
      close(unit=13)
	close(unit=14)
CC  
  197 continue
  	OPEN(12,FILE=OUTFILE,ACCESS='APPEND')
	OPEN(13,FILE=OUTFILEP,ACCESS='APPEND')
	OPEN(14,FILE=OUTFILEN,ACCESS='APPEND')

C
C  day-to-day cycle
C ========================================

       DO 500 idy=idy1,idy2
	lda=ndoy(iyr,mn,idy)  ! day-of-year
	call blet2(idy,ady)
	infile(30:31)=ady						  !Tamara
C
C	cnts(71,3,0:23),xind(71,3,0:23)

	do i=0,23		 ! ut=0,1,...,23
	do k=1,3
	do n=1,71
	cnts(n,k,i)=0.
	xind(n,k,i)=0.
	clat(n,k,i)=0.        !clat(71,3,0:23)
 	relat(n,k,i)=0.
 	enddo
	enddo
	enddo
C UT 
      do k=0,23
	res(k)=0.
	rup(k)=0.
	rlo(k)=0.
	resp(k)=0.
	resn(k)=0.
	enddo
C
 	sum=0.
	sumn=0.
	sump=0.

C cycle on UT:
C -----------------------------------------------------
      DO 801 iut=0,23		! cycle on UT
	call blet2(iut,aut)
	infile(32:33)=aut							 !Tamara
	OPEN(11,FILE=INFILE)
c
	cnttot=0.       !!!new for W-map
	cntneg=0.       !!!new for W-map
	cntpos=0.       !!!new for W-map
C ===== Input datafile 
 1311	DO 1799 jj=1,71	  ! glats= -87.5:2.5:87.5
c	pause ' '
    1	READ (11,12,END=3,ERR=2) (DHR(k),k=0,71) ! W(longs=-180:5:175)
  12	FORMAT(72(1X,F4.0))
C
	DO k=0,71				! Cycle on Longi	 <<<<<<<<
c==	if (dhr(k).gt.90.) dhr(k)=0.
C++ Count stormy indices W-:
	if (dhr(k).lt.-2) then 
	dcnt=1.					  ! weight=1 for ind=-3, -4
	cnts(jj,2,iut)=cnts(jj,2,iut)+dcnt	 !cnt(Wneg<-2) at lat jj
	cnts(jj,1,iut)=cnts(jj,1,iut)+dcnt	 !cnt(Wstr) at lat jj
C Counts at each Lati (divided by 72 glongs):
	clat(jj,2,iut)=clat(jj,2,iut)+dcnt/72. !glat,Ind,UT_iut
c!!!	cntneg=clat(jj,2,iut)
      cntneg=cntneg+1.  !!!new for W map
	clat(jj,1,iut)=clat(jj,1,iut)+dcnt/72.    ! Range
c!!!	cnttot=clat(jj,1,iut)
	cnttot=cnttot+1.  !!!new for W map
	 endif
C++	Count stormy indices W+:
		if (dhr(k).gt.2) then 
	dcnt=1.					  ! weight=1 for ind=+3, +4
 	cnts(jj,3,iut)=cnts(jj,3,iut)+dcnt	!cnt(Wpos>2) at lat jj
	cnts(jj,1,iut)=cnts(jj,1,iut)+dcnt	!cnt(Wstr) at lat jj
C Counts at each Lati:
		clat(jj,1,iut)=clat(jj,1,iut)+dcnt/72. ! Range
c!!!	cnttot=clat(jj,1,iut)
      cnttot=cnttot+1.  !!!new for W map
 	clat(jj,3,iut)=clat(jj,3,iut)+dcnt/72.	  ! W+
c!!!	cntpos=clat(jj,1,iut)
	cntpos=cntpos+1.  !!!new for W map
	 endif

C MAX index at 71 lati:
	if (rup(iut).lt.dhr(k)) rup(iut)=dhr(k)	! at current lat
	if (xind(jj,3,iut).lt.dhr(k)) xind(jj,3,iut)=dhr(k) !at lat=jj
	
C MIN index at 71 lati:
	if (rlo(iut).gt.dhr(k)) rlo(iut)=dhr(k) ! at current lat
	if (xind(jj,2,iut).gt.dhr(k)) xind(jj,2,iut)=dhr(k)	!at lat=jj
	enddo		    ! Cycle on Longi	 >>>>>>>>>>>>>
C   Range for given Lati, UT:
	amin=0.	 ! W<0.
	bmax=0.		!W>0.
	do k=0,71				   ! longs
	if ((dhr(k).lt.0.).and.(amin.gt.dhr(k))) amin=dhr(k) ! at lat=jj
	if ((dhr(k).gt.0.).and.(bmax.lt.dhr(k))) bmax=dhr(k) ! at lat=jj
	enddo

cc	wpos=bmax*(1.+clat(jj,3,iut)/35.) ! 	W+	/35?????????????
	wpos=bmax                            ! at lat=jj
	relat(jj,3,iut)=wpos	!W+			 ! at lat=jj
cc	wneg=amin*(1.+clat(jj,2,iut)/35.) !     W-	/35?????????????
      wneg=amin							 ! at lat=jj
	relat(jj,2,iut)=wneg    !W-			 ! at lat=jj
	if (clat(jj,1,iut).gt.0.) then
ctemp	wrange=(bmax-amin)*(1.+clat(jj,1,iut)) !  Range
cc	wrange=(bmax-amin)*(1.+clat(jj,1,iut)/35.) !  Range	/35?????????????
c!!!	wrange=(bmax-amin)*(1.+clat(jj,1,iut)/10.) !  Range	
      wrange=(bmax-amin)             		   ! at lat=jj
	                          else
	wrange=amax1(bmax,abs(amin))		   ! at lat=jj
	                          endif
	relat(jj,1,iut)=wrange  ! W Range	   ! at lat=jj

 1799	CONTINUE				   ! Cycle jj on glats (total W map)
  3   CLOSE(unit=11)

C----------------------------------------------------------
C average for 71 Lati => Basic Wp for each ut
C
		cneg=0.
 	cpos=0.
	ctot=0.
	 dneg=0.
	dpos=0.
	dtot=0.
	lat1=1    !cnew
	lat2=71	  !cnew
	alats=71.  !cnew				   ! 71 glati -87.5:2.5:87.5
c<<	grids=alats*72.          ! alats*24 hrs UT
      	grids=alats*72.          ! alats*glongs (=0:71 glongs)?????????????/
	do kk=lat1,lat2				   ! 
	dneg=dneg+relat(kk,2,iut)
	dpos=dpos+relat(kk,3,iut)
      dtot=dtot+relat(kk,1,iut)
	ctot=ctot+cnts(kk,1,iut)
	cneg=cneg+cnts(kk,2,iut)
	cpos=cpos+cnts(kk,3,iut)
C ++++++++++++
	enddo					   !  glati
	res(iut)=(1.+ctot/grids)*dtot/alats
	resn(iut)=(1.+cneg/grids)*dneg/alats
	resp(iut)=(1.+cpos/grids)*dpos/alats
	if (res(iut).ge.5.0) then
	totper=ctot/grids*100.
	write(*,73) AYR,AMN,ADY,iut,cneg,cpos,ctot,totper,res(iut)
	write(10,73) AYR,AMN,ADY,iut,cneg,cpos,ctot,totper,res(iut)
	endif
c   
  801 continue            ! end cycle on UT 
  73	format(1X,3A2,1X,I2,3(2X,F5.0),2(2X,F5.1))

C  Daily mean:
C
	sum=0.
	sump=0.
	sumn=0.
		do m=0,23      ! cycle on UT		======
 175	sum=sum+res(m)
	sump=sump+resp(m)
	sumn=sumn+resn(m)
	enddo							 !=========
ct	sum=sum/12.
	sum=sum/24.
	sump=sump/24.
	sumn=sumn/24.
	yrmndy='000000'
	yrmndy(1:2)=AYR
	yrmndy(3:4)=AMN
	yrmndy(5:6)=ADY	
	write(12,72) yrmndy,(res(k),k=0,23),sum
		write(13,72) yrmndy,(resp(k),k=0,23),sump
		write(14,72) yrmndy,(resn(k),k=0,23),sumn
	write(*,72) yrmndy,(res(k),k=0,23),sum
 	write(*,72) yrmndy,(resp(k),k=0,23),sump
 	write(*,72) yrmndy,(resn(k),k=0,23),sumn
C
  500	CONTINUE ! next day
C
   72 format(A6,24(1X,F4.1),1X,F7.2)
ct   75 continue	 	  ! <<<<<<<<<<<<<<<<<<
	CLOSE(unit=10)
	CLOSE(unit=12)
	CLOSE(unit=13)
	CLOSE(unit=14)
	goto 900
    2 CONTINUE
           write(*,*) 'INPUT FILE IS NOT IN YOUR DIRECTORY '
  900 continue
c      pause ' '
	RETURN
      END

C ============================================================================
      subroutine windex(xpar,xmed,delog,iw)
C
C   Procedure to calculate the ionospheric weather W index
C
C   Input: xpar = NeF2  - instantaneous value of the F2 layer peak electron density 
C      or  xpar = TECi  - instantaneous value of GPS-derived total electron content
C    		 xmed = NmF2  - quiet reference (median) value of the F2 layer peak electron density
C      or	 xmed = TEC m - quiet reference (median) value of total electron content
C
C   Output: ionospheric weather W index characterizing ionospheric and plasmaspheric activity:
C           W = 0, +1 or -1 (quiet state)
C           W = +2 or -2 (moderate disturbance)
C           W = +3 or -3 (moderate storm)
C           W = +4 or -4 (intense storm)
C
C   Includes logarithmic scale indices -4,-3,-2,-1,0,1,2,3,4
C   according to thresholds of 
C           log10(xpar/xmed)=-0.999,-0.301,-0.155,-0.045,0,0.045,0.155,0.301,0.999  
C
C   Ref.Gulyaeva, T.L.,, I. Stanislwska, and M. Tomasik. Ionospheric weather: Cloning missed 
C       foF2 observations for derivation of variability index. Annales Geophysicae, 26, N.2, 
C       315-321, 2008.
C       Gulyaeva, T.L.,and I. Stanislwska, Derivation of a planetary ionospheric storm index. 
C       Annales Geophysicae, 26, N.9, 2645-2648, 2008 .
C
      DIMENSION TR(7)
	INTEGER IIN(7)
   		data IIN/-3,-2,-1,0,1,2,3/           ! ID-index
	data TR/-.301,-.155,-.045,0.0
	+,.045,.155,.301/
C  Temp:
C    	 if (delog.ne.0.) goto 5   ! temp if delog is entered
	      delog=alog10(XPAR)-alog10(XMED)
   5		if (delog.gt.tr(7)) then
		iw=4
	goto 2
		endif 
		if (delog.eq.0.) then
			iw=0
	goto 2
	endif
	      if(delog.le.tr(1)) then
	iw=-4                 ! W-index
	goto 2
	endif
	if (delog.ge.tr(7)) then
	iw=4				  ! W-index
 	goto 2
	endif
C
	do 1 n=1,7
	if (delog.lt.0.) then
	if ((delog.ge.tr(n)).and.(delog.lt.tr(n+1))) then
	iw=iin(n)              ! W-index
	endif
	      else
	if ((delog.gt.tr(n-1)).and.(delog.le.tr(n))) then
	iw=iin(n)              ! W-index
	endif

	endif
   1	continue
C
   2  continue
C
	RETURN
	END
C
C
C ===============================================================
