program calpredict45ac
C.....................................................Jan  2024
C Produce ssn1(tau27) and F10.7(tau27) for current 
C
C,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,Mar. 2022
C Include outfile predict45.txt
C Include update of SSN2 and F10.7 for preceding day
C
C MOD45d model F107(fi) and SSN2
C using -45d(SC25) ~ -45d(SC24) with polifit coefficients
C add weight w
C weight w: yproxy=[yi(Eq.4a)-ymin(-45dSC25)]*rcoef
C rcoef=(ymax-ymin)/(ymax_pre-ymin_pre)
C
      dimension x(45),y(45)
	+,tyearp(0:90),jmnp(0:90),jdyp(0:90),txip(0:90),txindp(0:90) ! training set -45d:45d
	+,xyearc(0:90),imnc(0:90),idyc(0:90),xfic(0:90),xindc(0:90)  !modelling set -45d:45d
	CHARACTER*2 ANN,ANN_2 ! CS number
	+,AYR,AMN,ADY,AMN_pre,ADY_pre
	 	 CHARACTER*4 AYEAR,AYEAR_pre
	 CHARACTER*64 infile1,infile2,outfile,indate
C	 ,outpoli
C..	+ ,infiledf
	CHARACTER*1 SF		  ! SF='s' ~ SSN2; 'f' ~ F107.0
      CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
	 DIMENSION IM(12)
	DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
C
      CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
  799	format(A4,2(1X,A2),1X,A4,2(1X,A2),2(1X,I3))				 !FEB2015 NEW COMMAND NUMBER = 799
  798	format(' YEAR, AMN, DD2 = ',A4,2(1X,A2),1X,A4,2(1X,A2),2(1X,I3))
	indate='c:/web/graf/date'    !!! Tamara !!!!!!!	 !FEB2015
C#	indate='/var/www/izmiran/ionosphere/weather/graf/date' !!! LIUBA !!FEB2015
	open(10,file=indate)							   !FEB2015			
	read(10,799) AYEAR,AMN,DD2,AYEAR_pre,AMN_pre,ADY_pre
	+,lda_cur,lda_pre  ! FORMAT NUMBER = 799
	WRITE(*,798) AYEAR,AMN,DD2,AYEAR_pre,AMN_pre,ADY_pre
	+,lda_cur,lda_pre ! FORMAT NUMBER = 798
	close(unit=10)
C
	AYR=AYEAR(3:4)
       ADY=DD2
	read(AYR,*) ryr
	kyr=int(ryr)
	kyear=2000+kyr
	read(ADY,*) rdy
	kdy=int(rdy)
C
	z1=kyr/4.0
      jz=int(z1)*4
      IF(jz.EQ.kyr) THEN
	IM(2)=29
		  ELSE
                IM(2)=28
		       ENDIF
C
	read(AMN,*) rmn
	kmn=int(rmn)
	read(ADY,*) rdy_cur
 	kdy_cur=int(rdy_cur)
C+
	read(AMN_pre,*) rmn_pre
	kmn_pre=int(rmn_pre)
	read(ADY_pre,*) rdy_pre
 	kdy_pre=int(rdy_pre)
C++ Vsw, Eby
      call subbyvsw(kyear,kmn_pre,kmn,kdy_pre,kdy)
      call subbyvsw(kyear,kmn,kmn,kdy,kdy)
C	kdy=kdy_cur-2
	kdy=kdy_cur-1
C
	IF (kdy.le.0) THEN
      kmn=kmn-1
	if (kmn.eq.0) then
	kmn=12
	 kyr=kyr-1
           endif
	kdy=im(kmn)+kdy
	       ENDIF
C
C
C+      write(*,*) ' Enter NN of SC:'	! solar cycle for prediction
C+	read(*,*) nn
C
C     write(*,*) ' Update SSN2-phase COV-phase for prec. day:'
C     pause ' '
C        call subcursac
        call subcursbc
C
       nn=25         ! SC25
C+	nn_2=nn-2      ! Norm SC-2
 	nn_2=nn-1	 !Norm SC-1
	call blet2(nn,ANN)
	call blet2(nn_2,ANN_2)
C
C+      write(*,*) ' Enter <s> for SSN2 or <f> for F10.7:'
C+	read(*,*) SF
C          
C Start  using ' s':
	sf='s'
C++      write(*,*) ' Enter kyr:'
C++	read(*,*) kyr
	call blet2(kyr,AYR)
C++      write(*,*) ' Enter kmn:'
C++	read(*,*) kmn
	call blet2(kmn,AMN)
C++      write(*,*) ' Enter kdy:'
C++	read(*,*) kdy
		call blet2(kdy,ADY)
C..	infiledf='YRMMDD45DF.dat'
C	outpoli='calpolisp.txt'
  123	continue
C      OPEN(13,file=outpoli,access='APPEND')
C..       infiledf='YYMMDD45DF.dat'
C
C Continue first 's' then 'f': ----------------------------------
C
  	if (SF.eq.'f') then
       infile1='solphasecovNN.txt'  ! F10.7
	infile1(12:13)=ANN_2      ! training set
	infile2=infile1
	infile2(12:13)=ANN       ! modelling set
      outfile='covp45YRMMDDw.txt'    
		                else
	 infile1='solphasessnNN.txt'  ! SSN2
	infile1(12:13)=ANN_2      ! training set
	infile2=infile1
	infile2(12:13)=ANN       ! modelling set
	      outfile='ssnp45YRMMDDw.txt'    
      endif
C
C
	OPEN(11,file=infile1,action='READ')
	OPEN(12,file=infile2,action='READ') ! modelling set
C
C
	outfile(7:8)=AYR
	outfile(9:10)=AMN
	outfile(11:12)=ADY
	OPEN(10,file=outfile)
C
	lda_pre=ndoy(kyr,kmn,kdy)
C?	lda1=lda_pre+1			  !??????????
C?	lda2=lda_pre+45			  !??????????
C++
	lda_45=lda_pre-44  ! start of -45 prec days for modeling set
	kyr_pre=kyr
	if (lda_45.le.0) then
	kyr_pre=kyr_pre-1
C
	z1=kyr_pre/4.0
      jz=int(z1)*4
      IF(jz.EQ.kyr_pre) THEN
      ndnr_pre=366
		  ELSE
                IM(2)=28
	  		ndnr_pre=365
		       ENDIF
	lda_45=ndnr_pre+lda_45
	                 endif 
C++
	z1=kyr/4.0
      jz=int(z1)*4

      IF(jz.EQ.kyr) THEN
               IM(2)=29
	dnr=366.
      ndnr=366
		  ELSE
                IM(2)=28
	  	dnr=365.
	ndnr=365
		       ENDIF
C
      kyr_nex=kyr
	lda1=lda_pre+1
	kyr1=kyr
	if (lda1.gt.ndnr) then
	lda1=1
	      kyr_nex=kyr+1
	kyr1=kyr+1
	endif
C	lda2=lda_pre+45
	lda2=lda1+44

C
C	lda2_nex=lda2
      if (lda2.gt.ndnr) then
	lda2_nex=lda2-ndnr
	lda2=ndnr
	kyr_nex=kyr+1
	endif
C
	nd=0
      txi_pre=0.
      txind_pre=0.
	IF (kyr_pre.eq.kyr) GOTO 111
C
C ++++++++++++++++++++++++++++++Prec year
C
        iflag_pre=0
C
	   DO 49 lda=lda_45,ndnr_pre
	call submmdd(lda,kyr_pre,imn_pre,idy_pre)
C find day of modeling set: 
   47	read(12,*) xyearc(45),imnc(45),idyc(45),xfic(45),xindc(45) ! modelling set
      iyear=int(xyearc(45))
	imn=imnc(45)
	idy=idyc(45)
	iyr=iyear-iyear/100*100
	xfi=xfic(45)
	ldax=ndoy(iyr,imn,idy)
	if ((iyr.eq.kyr_pre).and.(ldax.eq.lda)) then
	goto 45
	                              else
C move data up
       do k=1,45
	xyearc(k-1)=xyearc(k)
	imnc(k-1)=imnc(k)
	idyc(k-1)=idyc(k)
	xfic(k-1)=xfic(k)
	xindc(k-1)=xindc(k)
	 enddo
	goto 47
	endif
   45	if (iflag_pre.eq.1) GOTO 49
C find days of training set
   46      read(11,*,end=51) tyearp(45),jmnp(45),jdyp(45),txip(45)	! training set
     +,txindp(45)
	 txi=txip(45)
C move data up
       do k=1,45
	tyearp(k-1)=tyearp(k)
	jmnp(k-1)=jmnp(k)
	jdyp(k-1)=jdyp(k)
	txip(k-1)=txip(k)
	txindp(k-1)=txindp(k)
	 enddo
C
	IF (xfi.ge.0) THEN   
C  Fi >0
	   if(txi.ge.xfi) then
	iflag_pre=1 
	goto 49
	               else
	goto 46
	endif
	               ELSE
C  Fi < 0
	   if(xfi.le.txi) then
	iflag_pre=1 
	goto 49
	               else
	goto 46
	endif
                    ENDIF
   49	continue	   ! End Prec year
C
C   +++++++++++++++++++++++++++++ Current year
  111 iflag=0     
C-      DO 50 lda=lda1,lda2
C-	call submmdd(lda,kyr,imn,idy) 
C?      lda=lda1
C find day of modeling set: 
    1	read(12,*,end=2) xyearc(45),imnc(45),idyc(45),xfic(45),xindc(45)	   ! modelling set
      iyear=int(xyearc(45))
	imn=imnc(45)
	idy=idyc(45)
	iyr=iyear-iyear/100*100
	xfi=xfic(45)
	ldax=ndoy(iyr,imn,idy)
	if ((iyr.eq.kyr1).and.(ldax.eq.lda1)) then
C	if ((iyr.eq.kyr_nex).and.(ldax.eq.lda)) then
	goto 2
	                              else
C move data up
       do k=1,45
	xyearc(k-1)=xyearc(k)
	imnc(k-1)=imnc(k)
	idyc(k-1)=idyc(k)
	xfic(k-1)=xfic(k)
	xindc(k-1)=xindc(k)
	 enddo
	goto 1
	endif
   55	if (iflag.eq.1) GOTO 50
C find days of training set
   2      read(11,*,end=51) tyearp(45),jmnp(45),jdyp(45),txip(45)		 ! training set
     +,txindp(45)
	 txi=txip(45)
C
	IF (xfi*txi.lt.0) goto 2  ! opposite signs of XHI mod/tfain
C
	IF (xfi.ge.0) THEN
C  Fi > 0
C	   if(xfi.ge.txi) then
	   if(txi.ge.xfi) then
	iflag_pre=1 
	goto 50
	               else
C move data up
       do k=1,45
	tyearp(k-1)=tyearp(k)
	jmnp(k-1)=jmnp(k)
	jdyp(k-1)=jdyp(k)
	txip(k-1)=txip(k)
	txindp(k-1)=txindp(k)
	 enddo
	goto 2
	endif
	               ELSE
C  Fi < 0
	   if(xfi.le.txi) then
	iflag_pre=1 
	goto 50
	               else
C move data up
       do k=1,45
	tyearp(k-1)=tyearp(k)
	jmnp(k-1)=jmnp(k)
	jdyp(k-1)=jdyp(k)
	txip(k-1)=txip(k)
	txindp(k-1)=txindp(k)
	 enddo
	goto 2
	endif
                    ENDIF
C
   50	CONTINUE
C
C========================================================
C read data of modeling set:
      	jj=46
       DO lda=lda1,lda2
	read(12,*) xyearc(jj),imnc(jj),idyc(jj),xfic(jj),xindc(jj)
	jj=jj+1
	 ENDDO
C Check next year output:
       IF (lda2_nex.eq.lda2) GOTO 51
C
C read param for the next year: +++++++++++++++++++
        DO lda=1,lda2_nex
C find day of modeling set: 
  	read(12,*,end=52) xyearc(jj),imnc(jj),idyc(jj),xfic(jj),xindc(jj)
      jj=jj+1
	  ENDDO
   51	close(unit=12)
C read next 45d data of training set:
         DO kk=46,90
      read(11,*,end=51) tyearp(kk),jmnp(kk),jdyp(kk),txip(kk)
     +,txindp(kk)
	 	 	   ENDDO
   52	close(unit=11)
C
C Find coefficient of polyfit-45
C
        do k=1,45
	x(k)=txindp(k)		  !-45d
	y(k)=xindc(k)		  ! -45d
C..	y(k)=DF(k)
	  enddo
C First polyfit before weight w
       call subpoly45(x,y,p1fo,p2fo,p3fo,rmseo,xaveo,yaveo)
	 xw=txindp(46)
	fnmo=xw*(p1fo*xw+p2fo)+p3fo  ! Polyval
	weig1=xindc(45)/fnmo    ! Weight for polyfit
	xw2=txindp(90)
		weig2=xw2/yaveo
C+
CNEW      weig3=yaveo/xaveo
       weig3=1.0
	dweig=(weig1-weig3)/45.
Cpre	dweig=(weig1-weig2)/45.
C+       weig=yaveo/xaveo	   ! Weight for polyfit
	write(*,*) AYR,AMN,ADY,' p1o ',p1fo,' p2o= ',p2fo,' p3=o',p3fo
	+,' rmso=',rmseo,xaveo,yaveo,weig1,weig3,dweig
C	write(13,*) AYR,AMN,ADY,' p1o ',p1fo,' p2o= ',p2fo,' p3=o',p3fo
C	+,' rmso=',rmseo,xaveo,yaveo,weig1,weig3,dweig
C	GOTO 1001 ! Avoid sorting
C------------------------------------------
C add RMSE calculation
      sum2=0.
C-----------------------------------------------------
C produce next 45d data of modelling set:
C>> 1001     weig=weig1
 1001 CONTINUE
C?      zmin=1000.
C?	zmax=0.
        weig=weig1
 	DO jj=46,90
C
      	weig=weig-dweig
C
	xx=txindp(jj)  !SA index of training set
C	xx=xindc(jj-45)  ! SA-45d index of modeling set
C	xx=DF(jj)
	fnmo=xx*(p1fo*xx+p2fo)+p3fo  ! Polyval
C--		fnm=xx*(p1f*xx+p2f)+p3f  ! Polyval
C--	fnmres=(fnmo+2.*fnm)/3.
        fnmres=fnmo*weig
C++	weig=weig-dweig
C>	if (fnm.gt.(yave+rmse)) then
C>	fnm=yave+rmse
C>	goto 70
C>	endif
C>	if (fnm.lt.(yave-rmse)) then
C>	fnm=yave-rmse
C>	endif
C
C   70	xindc(jj)=fnm  ! => 46...90
      if ((SF.eq.'f').and.(fnmo.lt.65.)) fnmo=65.
      if ((SF.eq.'s').and.(fnmo.lt.0.)) fnmo=0.
      if ((SF.eq.'f').and.(fnmres.lt.65.)) fnmres=xindc(jj-1)
      if ((SF.eq.'s').and.(fnmres.lt.0.)) fnmres=xindc(jj-1)
	sum2=sum2+(xindc(jj)-fnmres)**2.
   70	xindc(jj)=fnmres  ! => 46...90
C//      if (zmin.gt.fnmo) zmin=fnmo
C//      if (zmax.lt.fnmo) zmax=fnmo
	ENDDO
	rmse=sqrt(sum2/45.)
C
C Output results
C       DO m=45,89
       DO m=46,90
	write(*,30) xyearc(m),imnc(m),idyc(m),xfic(m),xindc(m)
	write(10,30) xyearc(m),imnc(m),idyc(m),xfic(m),xindc(m)
	 ENDDO
C
C 81-aver for the 1st day
C
        sum=0.
	do n=6,86
	sum=sum+xindc(n)
	enddo
	ave=sum/81.
	write(*,30) xyearc(45),imnc(45),idyc(45),xfic(45),ave,weig1
	write(10,30) xyearc(45),imnc(45),idyc(45),xfic(45),ave,weig1
C
      close(unit=10)
C	      close(unit=13)
	if (SF.eq.'s') then
	SF='f'
	GOTO 123
	             else
	call subpred45d(kYR,kMN,kDY)
	endif
C
   30 format(1X,F8.3,2(1X,I2),1X,F8.4,2(1X,F6.2))
C
	pause ' '
      STOP
	END
c--------------------------------------------
	real function harmcov(x)
C............................................................Jan 2022
C Furie model N=6
C
      DIMENSION a(0:6),b(6)
		DATA a/-1.147,3.704,-2.827, 2.02,-1.222, 0.5428, -0.1046/
      DATA b/ 0.09036, 0.01687, 0.00798, -0.07502, -0.004809, 0.02875/
       w =  2.155
      sind =a(0) + a(1)*cos(x*w) + b(1)*sin(x*w) + 
     * a(2)*cos(2.*x*w) + b(2)*sin(2.*x*w) + a(3)*cos(3.*x*w) + 
     + b(3)*sin(3.*x*w) + a(4)*cos(4.*x*w) + b(4)*sin(4.*x*w) + 
     + a(5)*cos(5.*x*w) + b(5)*sin(5.*x*w) + 
     * a(6)*cos(6.*x*w) + b(6)*sin(6.*x*w)
	 harmcov=sind
	 RETURN
	END