C TLG....................................................July 2020
C
C Simulate missing monthly-hourly foF2, hmF2 for UT=0,2,...,22 hrs 
C
C
C 
C Output (ADD):  correlation coefficient ro, RMSE
C
      PROGRAM EXTRAHFT2pr
C---------------------------------------------------------------
C data for 24 UT=0,1,2,...23hUT daily
C
C
C*********TEST
       DIMENSION glats(20),glons(20)
	CHARACTER*128 infilef,infileh,outfile,txt
C	+,outfile_h,outfilemh ! _h (artificial data), mh (source data)
	character*2 ayr,amn,ST,Cf(0:23),Ch(0:23),CC
	+,ADY,AUT,ADY2,BYR,BMN,BDY,STS_OUT(20)
	CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
	CHARACTER*1 fht
	CHARACTER*13 STAT
	CHARACTER*4 TITf,TITh,CTIME,AYEAR
		character*5 STCOD
	+,NUM_STS(20)
	dimension im(12),RESF(0:23),RESH(0:23),x(24),y(24),ynm(24)
	+,IRESF(0:23),IRESH(0:23),xfof2(372),yfof2(372),xhmf2(372)
     +,yhmf2(372)
	DATA  IM/31,28,31,30,31,30,31,31,30,31,30,31/
      DATA STS_OUT/
	+ 'at','bc','cb','ch','db','eb','gr','jm'
     +,'jr','kk','kw','mo','nc','pq','ps','qn'
     +,'ro','ss','so','wk'/
      DATA NUM_STS/
	+ 'AT138','BC840','CB53N','RL052','DB049','EB040','GR13L','JI91J'
     +,'JR055','TO536','KJ609','MO155','NI135','PQ052','PSJ5J','THJ77'
     +,'RO041','SMJ67','SO166','WK545'/  
      DATA glats/
	+  38.0, 40.0,-35.3, 51.6, 50.1, 40.4,-33.3,-12.0
     +, 54.6, 35.7,  9.0, 55.5, 35.0, 50.0,-51.7, 77.5
	+, 41.8, 67.0, 67.4, 45.4/
      DATA glons/
     +  23.6,254.7,149.0,358.7,  4.6,  0.5, 26.5,283.2
	+, 13.4,139.5,167.2, 37.3, 33.2, 14.6,302.2,290.8
	+, 12.5,310.0, 26.6,141.7/
C
   79	FORMAT(1X,A13)
c
	      CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT
C
       CC='_'
	TITf='foF2'
	TITh='hmF2'
	CTIME='  UT'
	CTIME(1:2)=TT(1:2)
C read stations coordinates
  76	format(1X,A4,A2,2(1X,F5.1),1X,A13,1X,A2)
C
  	      WRITE(*,*)' ENTER AYEAR:'
		READ(*,15) AYEAR
 201	      WRITE(*,*)' ENTER AMN <00> STOP:'
		READ(*,16) AMN
C
	if (AMN.eq.'00') stop
   15	FORMAT(A4)
   16	FORMAT(A2)
160   FORMAT(A128)
	AYR=AYEAR(3:4)
	read(ayr,*) yr
	iyr=int(yr)
	if(int(yr/4.)*4.eq.iyr) THEN 
	im(2)=29
	ndnr=366
	                          ELSE
	IM(2)=28
	ndnr=365
	ENDIF
	read(amn,*) xmn
	mn=int(xmn)
C
	ndnr=im(mn)
	call blet2(ndnr,ADY2)
	ixf=0
	iyf=0
	ixh=0
	iyh=0
  202	continue
C
C start foF2, hmF2 or TEC
C
	infilef='d:\web\dat1\YR\styrmn.txt'
	infilef(13:14)=AYR
	infilef(18:19)=AYR
	infilef(20:21)=AMN
	infileh=infilef
	infileh(11:11)='6'
	infileh(22:)='h.txt'
C
      outfile='corr2.txt'
C
      OPEN(13,file=outfile,access='APPEND')
C
C Cycle on stations
C
       DO 51 ist=1,20
	st=STS_OUT(ist)
	stcod=NUM_STS(ist)
	glat=glats(ist)
	glon=glons(ist)
	infilef(16:17)=st
	infileh(16:17)=st
C
C
   61	format(3A2,1X,24(A4,1X))
C
	write(*,79) stat
	if ((st.eq.kw).or.(st.eq.ss)) GOTO 51        !ADD 2014
      if ((st.eq.qn).and.(mn.eq.10)) GOTO 51		 !ADD 2014
	if ((st.eq.eb).and.(mn.eq.7)) GOTO 51		 !ADD 2014

C
		    OPEN(11,FILE=INFILEF,action='READ')
	OPEN(12,FILE=INFILEH,action='READ')
	DO i=1,4
	read(11,160,end=773) txt            !NEW
	if (i.eq.1) then
	stat=txt(1:13)
	endif
	read(12,160,end=774) txt			 !NEW
	ENDDO
 773	close(unit=11)						!NEW
 774  close(unit=12)						!NEW
      GOTO 51								!NEW
C
	write(*,76) AYEAR,AMN,glat,glon,stat,st
	write(13,76) AYEAR,AMN,glat,glon,stat,st
	do i=0,23
	iresf(i)=0
	iresh(i)=0
	enddo
  200	FORMAT(3A2,2X,24(I3,A2))
C
   62 FORMAT(1X,4A2,1X,F6.1)
	ixf=0
	iyf=0
	ixh=0
	iyh=0
      DO 777 idy=1,ndnr
	call blet2(idy,ADY)
  601	READ(11,200,end=771) BYR,BMN,BDY,(iresf(k),Cf(k),k=0,23)
  	READ(12,200,end=772) BYR,BMN,BDY,(iresh(k),Ch(k),k=0,23)
	if ((AYR.eq.BYR).and.(AMN.eq.BMN).and.(ADY.eq.BDY)) then
	goto 302
	                else
	goto 601
	        endif
C
  302	CONTINUE
      DO iut=0,23
      resf(iut)=iresf(iut)/10.
	resh(iut)=float(iresh(iut))
	call blet2(iut,AUT)
	fht='t'		! TEC
	call subinterp1(fht,glat,glon,byr,bmn,bdy,aut,res)
      x(iut+1)=res/10.
      ENDDO
C
  	DO mut=0,23,2
	fm=resf(mut)
	call blet2(mut,AUT)
C Obs foF2:
	ixf=ixf+1
      xfof2(ixf)=fm
C	write(*,62) BYR,BMN,BDY,AUT,fm
C
	hm=resh(mut)
C Obs hmF2:
	ixh=ixh+1
      xhmf2(ixh)=hm
c	write(*,62) BYR,BMN,BDY,AUT,hm
C
	fht='f'		! foF2
	call subinterp1(fht,glat,glon,byr,bmn,bdy,aut,res)
	resf(mut)=res/10.
C
	fht='h'		! hmF2
	call subinterp1(fht,glat,glon,byr,bmn,bdy,aut,res)
	resh(mut)=res
	ENDDO
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C polyfit foF2 and hmF2
C
      fht='f'        ! foF2
  305	      do k=1,24
	fof2=resf(k-1)
 	ynm(k)=1.24*fof2**2  ! NmF2
	y(k)=ynm(k)  !First NmF2
	enddo
       call subpoly3c(x,y,p1f,p2f,p3f,rmsef)
	fht='h'
	      do k=1,24
	y(k)=resh(k-1) ! hmF2
	enddo
      call subpoly3c(x,y,p1h,p2h,p3h,rmseh)
C Output profile:
	DO jj=0,23,2
	call blet2(jj,AUT)
	xx=x(jj+1)  !TEC
	yproxf=xx*(p1f*xx+p2f)+p3f  ! Polyval
	yproxh=xx*(p1h*xx+p2h)+p3h  ! Polyval
	  if(yproxf.ge.4.) then
	yproxf=sqrt(yproxf/1.24)
	else
	yproxf=resf(jj)
	endif
      	if ((yproxh.lt.175.).or.(yproxh.gt.500.)) then
     	    yproxh=resh(jj)
	    endif
C
C Proxy foF2:
	iyf=iyf+1
      yfof2(iyf)=yproxf
C	write(*,62) BYR,BMN,BDY,AUT,yproxf
C
C Proxy hmF2:
	iyh=iyh+1
      yhmf2(iyh)=yproxh
c	write(*,62) BYR,BMN,BDY,AUT,yproxh
C
	ENDDO
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C 200	format(3A2,24(1X,I4))
  777 CONTINUE  ! next day
  771 close(unit=11)
  772 close(unit=12)
c ***************************************************************
C
 111	format(A6,1X,24(I4,1X))
	GOTO 33
C
    2 CONTINUE
       write(*,*) 'INPUT FILE IS NOT IN YOUR DIRECTORY '
 	pause ' '
	goto 34

   33	continue
C
C  Add Calculation of corr-coef and RMSE
C (1) between xfof2 & yfof2 (DO ii=1,ndnr
	ncor=ndnr*12
      call subrmscor(ncor,xfof2,yfof2,corf,rmsf)
C (2)	between xhmf2 & yhmf2 (DO ii=1,ndnr
      call subrmscor(ncor,xhmf2,yhmf2,corh,rmsh)
C write(13,...) !results <corfoF2,rmsefoF2, corhmF2,rmsehmF2> write to output file
	write(*,133) st,AYEAR,AMN,corf,rmsf,corh,rmsh
	write(13,133) st,AYEAR,AMN,corf,rmsf,corh,rmsh
   51	CONTINUE             ! next station
      close(unit=13)
	GOTO 201		   ! next month
  133	FORMAT(1X,A2,1X,A4,A2,4(1X,F7.3))
      pause ' ' 
   34   STOP
       END

C***********************************************************
C=========================================
      subroutine subpoly3c(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

C+++++++++++++++++++++++++++++++++++++++++++++++++++++++