      program comparmedfht
C.......................................................Apr 2020
C
C Compare model IRI-Plas-TEC(foF2,hmF2 from GMF2), with observation
C monthly median data / model
C
      character*120 infile1,infile2,outfile,txt
		character*2 AYR,st
	character*4 AYEAR
      dimension CC(31,24),xmed(24),IM(12)
     +,tobs(288),tmod(288)
	DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
	common /BL1/CC,XMED
C
C+      infile1='stYEAR_1.out' ! CCIR	model
	infile1='stYEARfh_2.out' ! GMF2	model
C
	infile2='d:\web\graf\YEAR\styr.med' ! TEC-jpl data obs
C
C+	AYEAR='2014' 
      AYEAR='2019' 
C+      AYR='14'
      AYR='19'
	read(ayr,*) ryr
	iyr=int(ryr)
	if(int(ryr/4.)*4.eq.iyr) THEN 
 	im(2)=29
	                          ELSE
	im(2)=28
	ENDIF
c
   31 format(2(1X,A25))
	infile1(3:6)=AYEAR
C
      infile2(13:16)=AYEAR
	infile2(20:21)=AYR
C
      outfile='comparmedfht.res'
	OPEN(10,file=outfile,access='APPEND')
C
  1      write(*,*) ' Enter st 00 for stop:'
	read(*,*) st
	if (st.eq.'00') goto 77
       infile1(1:2)=st
	infile2(18:19)=st

	write(*,31) infile1,infile2
	write(10,31) infile1,infile2

	OPEN(11,file=infile1,action='READ')
	OPEN(12,file=infile2,action='READ')
	do k=1,7
	read(11,*) txt   ! title line
	enddo

C
C Read Model results for month:
       tav=0.
       DO jj=1,288
	read(11,33,end=51) tmod(jj)
	read(12,34,end=52) tobs(jj)
	tav=tav+tobs(jj)
	enddo
	tav=tav/288.
C
C RMSE  and PE
       dev=0.
	zet=0.
	DO jj=1,288
	dev=dev+(tobs(jj)-tmod(jj))**2
	zet=zet+(tobs(jj)-tav)**2
	ENDDO
	RMSE=sqrt(dev/288.)
	PE=(1.-sqrt(dev/zet))*100. ! %
C
        write(*,35) AYEAR,RMSE,PE
	write(10,35) AYEAR,RMSE,PE
   70 CONTINUE ! 
   51   close(unit=11)
   52   close(unit=12)
	goto 1
c++++++++++++++++++++++++++++++++++++++++++++++++++++++
   77	close(unit=10)
C
   33 format(144X,F6.2) ! TEC from IRI-Plas.out file
   34 format(26X,F5.1) ! TEC-jpl obs
   35 format(1X,A4,1X,F7.2,1X,F8.3) ! tec RMSE PE
       pause ' '
	STOP
	END