C*************************************************************   
CC 
      SUBROUTINE SUBTAU81(IYYYY,IMN,IDY,Ftau,Rtau)
CC-------------------------------------------------------------------
CC......................................................Apr. 2024
CC Ytau27=(1-tau)/(1-tau**n)*(F0+tau*F_1+tau**2*F_2...+tau**29*F_26) =smoothed SSN1(tau) and F10.7(tau) index
CC
C Ref. Deminov M.G. Geomagn. Aeronomy, 2022, V.62, No. 3, P.302-306 
C
CC Modified to produce sunspot number R81, and solar radio flux F81 
CC from two successive kpYEAR annual files
CC Input: 
CC       kpYEAR  (preceding year: IYYYY-1)                              
CC       kpYEAR  (given year: IYYYY)                              
c--------------------------------------------------------------------
c 
c T.L. Gulyaeva ......................................  Sep. 2001. 
c Modified for producing solar radio flux F10.7, F81, averaged for
c 81 days preceeding given day
C 
C------------------------------------------------------------------------
      DIMENSION lm(12),icurkp(0:7),iap0(8)
     &,icv(0:82),irz(0:82)
C-     ,CC(81)
      CHARACTER*4 YEAR,YEAR_pre
      CHARACTER*80 infilekp
C?      COMMON /FIKP1/iyear_cur,imn_cur,idy_cur
C-	COMMON /BL1/CC,SMED,SMIN,SMAX
      DATA LM/31,28,31,30,31,30,31,31,30,31,30,31/
C
C START: 
      ipcnt=0 ! keep cnt -28,...,-1
C	tau=0.9636
      tau=exp(-1./27.)
         do k=0,82
       irz(k)=0
       icv(k)=0
         enddo
	iyear_cur=iyyyy
	imn_cur=imn
	idy_cur=idy
C   Input file kpYEAR starts from 1948:
C      if((iyyyy.lt.1948).or.(iyyyy.gt.iyear_cur)) goto 21 !*
C      IF ((iyyyy.eq.iyear_cur).and.(imn.gt.imn_cur)) goto 21 !*
C      IF ((iyyyy.eq.iyear_cur).and.(imn.eq.imn_cur).and.    !*
C     *   (idy.gt.idy_cur)) goto 21 !*TEST
C
        iy=iyyyy-1900

 770       z1=iy/4.0
        lm(2)=28
           if(iyyyy/4*4.eq.iyyyy) lm(2)=29
      if (iy.ge.100) then
      iyr=iy-100
      else
      iyr=iy
      endif
      call blet4(iyyyy,YEAR)
      iYYYY_pre=iYYYY-1
      call blet4(iYYYY_pre,YEAR_pre)
C
      idend=idy
C
C Day-of-year:
      ldaend=ndoy(iyyyy,imn,idend)       ! current day-of-year
C
      ipre=0
       if (ldaend.lt.82) then
        ipre=1    ! Include data for preceding year
       endif
C
Cpre      infilekp='kpYEAR'  
C-       infilekp='c:\web\graf\kpYEAR' 
       infilekp='c:\MSDEV\Projects\Isomain\kpYEAR'  
C
10    FORMAT(3I2,6X,8(I2),3X,8(I3),7X,I3,I3,1X,I1)

C
C-------------------------------------------------------------------------------
C Go to read current year file kpYEAR (avoid preceding year input) if current day IDY>81 day-of-year:
C
      if (ipre.eq.0) GOTO 203 !avoid input of preceding year file kpYEAR>>>>>
C
C-------------------------------------------------------------------------------
C Include Input file 'kpYEAR' for preceding year 
C
C      infilekp(3:6)=YEAR_pre  !!! Tamara
      infilekp(29:32)=YEAR_pre  !!! Tamara

      OPEN(19,FILE=infilekp,STATUS='OLD',ERR=21)  
  201 continue

      READ(19,10,END=202) JYR,JMN,JDY,(icurkp(k),k=0,7)
     &,(iap0(i),i=1,8),ir1,if1,if2
C
C Replace missed F10.7 by regression model:
C-      if (if1.eq.0) then
C-      cov=63.8255+0.8868*ir1  ! daily F10.7(Rz) model
C-      icv(28)=nint(cov*10.)
C-      else
      icv(82)=if1*10+if2  !F107*10
C-      endif
      irz(82)=ir1
C Move data 1 line up:
      do k=0,81
      irz(k)=irz(k+1)
      icv(k)=icv(k+1)
      enddo
      GOTO 201
  202 close(unit=19)
      ipcnt=82-ldaend   ! 
C
C For current year:
C  203 infilekp(3:6)=YEAR   !
  203 infilekp(29:32)=YEAR   !
C
C      OPEN(13,FILE=infilekp,STATUS='OLD',ERR=21)
      OPEN(19,FILE=infilekp,action='READ')
C
C Read file kpyear for current year:
C
  204 continue  
C
      READ(19,10,END=30) JYR,JMN,JDY,(icurkp(k),k=0,7)
     &,(iap0(i),i=1,8),ir1,if1,if2
      irz(82)=ir1
c-      if (if1.eq.0) then
c-       cov=63.8255+0.8868*ir1  ! daily F10.7(Rz) model
c-       icv(28)=nint(cov*10.)
c-      else
       icv(82)=if1*10+if2  !F107*10
c-      endif
C Move data 1 line up:
  140 do  k=0,81
      irz(k)=irz(k+1)
      icv(k)=icv(k+1)
      enddo
C
      ipcnt=ipcnt+1

      if (ipcnt.lt.82) goto 204

      IF ((JYR.eq.IYR).and.(JMN.eq.IMN).and.(JDY.eq.idend)) then 
      goto 207
      ELSE
      ipcnt=ipcnt-1
      goto 204
      ENDIF

C data for 81 preceding days are collected                       
  207 Ftau=0.0
      Rtau=0.0
C-	fsum=icv(1)/10.
C-	rsum=irz(1)/1.
      	fsum=icv(81)/10.
		rsum=float(irz(81))
C?      fsum=0.
C?	rsum=0.
	stau=tau
	stausum=1.
C-      do k=1,26
      do k=80,0,-1
	fsum=fsum+stau*icv(k)/10.
	rsum=rsum+stau*irz(k)
	stau=stau*tau
	stausum=stausum+stau
	enddo
	coef=(1.-tau)/(1.-stau)
	Rtau=coef*rsum
	Ftau=coef*fsum
	Rtau=rsum/stausum
	Ftau=fsum/stausum
C
C
      goto 30
21      write(*,100)
100     format(1X,'Date is outside range of Kp-Ap indices file.')
      goto 30
  29  write(*,*) 'Error in input k-index file',infilekp 
C      pause ' '
      stop
  30  CLOSE(unit=19)
      RETURN
      END

C __________________________________________________________________________
C