C............................................................Jan. 2021
C Add polyfit between day_pre & day_cur to fill missed foF2 data
C
C
C       "MEDIAN" 
      program cntfstatp
C
      DIMENSION im(12),CC(15,0:23),CINP(0:15,0:23),icnut(0:23)
     +,xmed(0:23),iw(0:23),x24(24),y24(24),xmon(31,0:23),imedmn(0:23)
	+,ires(0:23),icnt(0:23),XMEDMN(0:23)
      character*2 ayr,amn,ady,st,AYR_pre,c2(0:15,0:23),aa,CMON(31,0:23)
	character*45 infile1,infile2,outfile,outfile2 ! foF2+missed _ proxy foF2
	character*128 txt
		character*1 a1,a2
	character*4 AYEAR
		character*5 blank,blanks(0:23)
		COMMON /BL1/CC,XMED
		COMMON /BL2/XMON,XMEDMN
	DATA  IM/31,28,31,30,31,30,31,31,30,31,30,31/
c
C
   97 format(1X,A2)
      blank='     '
	do i=0,23
	blanks(i)=blank
	enddo
C
       write(*,*) 'ENTER st:'
	 read(*,*) st
       write(*,*) 'ENTER Year0-start,YEAR1-end:'
	 read(*,*) jyear0,jyear1
       write(*,*) 'ENTER imn0:'
	 read(*,*) imn0
Crem	st='mo'			  !Moscow
	infile2='d:\data\IWG\st\stYR'
	infile2(13:14)=st
	infile2(16:17)=st
	infile1=infile2
	outfile2=infile2
	outfile2(20:)='f.txt'
C-----------------------------------
	write(*,97) st
	
C
Crem	DO 400 jyr=1946,1998
C+	DO 400 jyr=jyear,2015
	DO 400 jyyyy=jyear0,jyear1
	call blet4(jyear,AYEAR)
	iyr=jyyyy-(jyyyy/100)*100 
	if (jyyyy.eq.jyear0) then
	imn1=imn0
	else
	imn1=1
	endif

C
CCC      jyr_pre=iyr-1  !
      jyr_pre=iyr           ! Starting from median for 15 days of starting month
	call blet2(jyr_pre,AYR_pre)
	infile1(18:19)=AYR_pre
	OPEN(11,file=infile1,action='READ')
C
  111	FORMAT(A120)
  107	FORMAT(24(F3.1,A2))
  117	FORMAT(24(I3,2X))
  127	FORMAT(24(I3,A2))
  137	FORMAT(24(A5))
C Read Last 14 days of prec.year
C-      DO 106 ii=1,12
	read(11,111,end=109) txt ! title line
C-	idmn1=im(ii)
C
	    do kk=1,15
	read(11,107,end=109) (CINP(15,k),C2(15,k),k=0,23)
C move data one line up:
	do n=1,15
	do m=0,23
	CINP(n-1,m)=CINP(n,m)
	C2(n-1,m)=C2(n,m)
      CC(n-1,m)=CINP(n-1,m)
	enddo
	enddo
      	enddo  
  106 CONTINUE
  109 close(unit=11)
       call rmedian15(15)
c
c Current year	
      call blet4(jyyyy,AYEAR)
	AYR=AYEAR(3:4)
	yr=float(jyyyy)
	if(int(yr/4.)*4.eq.jyyyy) THEN 
	im(2)=29
	ndnr=366
	                          ELSE
	IM(2)=28
	ndnr=365
	ENDIF
	infile2(18:19)=AYR
	OPEN(12,file=infile2,action='READ')
C
	outfile2(18:19)=AYR
	OPEN(13,file=outfile2)
C
	outfile=infile2
	outfile(20:)='wf.txt'
	OPEN(10,file=outfile)
C
C monthly cycle
C
       DO 410 imn=imn1,12
C
	read(12,111,end=108) txt ! title line
	write(13,111) txt ! title line
	call blet2(imn,AMN)
C
         ndmn=im(imn)
	    do 444 kk=1,33		   !+++++++++++++++++++++++++++++++++++++++++
	if (kk.le.ndmn) then		 !!!!!!!!!!!!!!!!!!
	read(12,107,end=108) (CINP(15,k),C2(15,k),k=0,23)
	idy=kk
	call blet2(idy,ADY)
C 
      do k=0,23
	icnut(k)=0
	enddo
	do m=0,23
      CC(15,m)=CINP(15,m)   ! Current day
      xmon(idy,m)=CINP(15,m)   ! Current day
	CMON(idy,m)=C2(15,k)     ! Current day
	if (CC(15,m).eq.0.) icnut(m)=icnut(m)+1 
	enddo
	mflag=0
	      do k=0,23
 	if (icnut(k).eq.15) mflag=1
	enddo
C
	if (mflag.eq.0) call rmedian15(15)
C Wf index
       iflag=0
      do 442 ll=0,23					!============
	if ((CC(15,ll).eq.0.).and.(iflag.eq.0)) then         !>>>>>>>>>>>>
	iflag=1
C
C Fill missed foF2 by regression
       DO mm=0,23					 !.............
Crem	 x24(mm+1)=CC(13,mm)
Crem	y24(mm+1)=CC(14,mm)
	 x24(mm+1)=CC(14,mm)
	if(x24(mm+1).eq.0.) x24(mm+1)= xmed(mm)
	if (CC(15,mm).eq.0.) then
	y24(mm+1)=xmed(mm)
	else
	y24(mm+1)=CC(15,mm)
	endif
	 ENDDO						  !..............
C
	call subpoly3p(x24,y24,p1f,p2f,p3f,rmse)
	write(*,*) ' pi=',p1f,p2f,p3f,rmse
	endif						  !<<<<<<<<<<<<<<<<
C
	if (CC(15,ll).eq.0.) then       !+++++++++
      xx=CC(14,ll)
	yprox=xx*(p1f*xx+p2f)+p3f  ! Polyval
	if (yprox.lt.1.5) yprox=1.5
	cc(15,ll)=yprox
	aa=C2(15,ll)
	a1=aa(1:1)
	a2=aa(2:2)
	  aa(1:1)='_'
	  if (a1.ne.' ') then
      	if (a2.ne.' ') then
	    goto 123
	    else
	    aa(2:2)=a1
	    endif
	  endif
  123	C2(15,ll)=aa  
C	               
      XMON(idy,ll)=cc(15,ll)
	CMON(idy,ll)=C2(15,ll)
	CINP(15,ll)=cc(15,ll)
	endif                 !+++++++++
	ires(ll)=nint(CC(15,ll)*10.0)
	xpar=cc(15,ll)*CC(15,ll)
	xme=xmed(ll)*xmed(ll)
	delog=0.
	call windex(xpar,xme,delog,iwf)
      iw(ll)=iwf
  442 CONTINUE						   !================
	write(*,112) AYR,AMN,ADY,(iw(k),k=0,23)
      write(10,112) AYR,AMN,ADY,(iw(k),k=0,23)
C
	write(*,127) (ires(k),C2(15,k),k=0,23)
      write(13,127) (ires(k),C2(15,k),k=0,23)
	goto 441
	else
	read(12,111,end=108) txt ! title line
	if (kk.eq.32) then	  !=+=+
	call rmedmon(ndmn)
	do k1=0,23
	imedmn(k1)=nint(XMEDMN(k1)*10.)
	icnt(k1)=ndmn
	enddo
	write(*,117) (imedmn(k),k=0,23)
	write(13,117) (imedmn(k),k=0,23)
	goto 444
	endif					  !=+=+
	if (kk.eq.33) then
	write(*,117) (icnt(k),k=0,23)
	write(13,117) (icnt(k),k=0,23)
	goto 444
	endif
	write(*,137) (blanks(k),k=0,23) ! blank line
	write(13,137) (blanks(k),k=0,23) ! blank line
	goto 444
	endif	!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C move data one line up:
  441		do m=0,23
      do n=1,15
			CINP(n-1,m)=CINP(n,m)
				C2(n-1,m)=C2(n,m)
		enddo
	enddo
         do k=0,23
		do n=1,15
		CC(n,k)=CINP(n,k)
	enddo
	    enddo
  444 CONTINUE          !+++++++++++++++++++
  410 CONTINUE
  108	close(unit=12)
	close(unit=10)
	close(unit=13)
  112 FORMAT(3A2,24(1X,I4))
C
  400	continue  ! next year
C
   33	continue
      write(*,*) infile2,outfile
c-	goto 1
   34   pause ' ' 
      STOP
       END

C++
c----------------------------------------------------------------------
      subroutine subpoly3p(x,y,p1,p2,p3,rmse)
	dimension x(24),y(24) 
	nut=24  ! diurnal UT set
C  
	s0=24.
	s1=0.
	s2=0.
	s3=0.
	s4=0.
	syx2=0.
	syx=0.
	sy=0.
      DO N=1,nut
	x1=x(N)
	x2=x1*x1
	x3=x1*x2
	x4=x2*x2
	y1=y(N)
	yx=y1*x1
	yx2=y1*x2
	s1=s1+x1
	s2=s2+x2
	s3=s3+x3
	s4=s4+x4
	syx2=syx2+yx2
	syx=syx+yx
	sy=sy+y1
      ENDDO 
C
      d1=(s0*s3-s1*s2)
      c1=(s0*s4-s2*s2)/d1
      c2=(s0*syx2-sy*s2)/d1
	d2=(s0*s2-s1*s1)
      c3=(s0*s3-s1*s2)/d2
      c4=(s0*syx-sy*s1)/d2

C
C p1=a; p2=b; p3=c
C
C polyfit2 coefficients
	p1=(c2-c4)/(c1-c3)
	p2=c2-p1*c1
	p3=(sy-p1*s2-p2*s1)/s0
C RMSE
      sum=0.
	DO N=1,nut
	del=y(N)-(x(N)*(p1*x(N)+p2)+p3)
	sum=sum+del**2
	ENDDO
	rmse=sqrt(sum/nut)
C
	RETURN
	END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++