      program comparmedhm
C.......................................................Mar 2020
C
C Compare Shubin model hmF2 with observation
C monthly median data / model
C
      character*120 infile1,infile2,outfile,outmed,txt
		character*2 AYR,st,AMN
	character*4 AYEAR
	character*6 YYMMDD
      dimension CC(31,24),xmed(24),IM(12),ymod(24)
     +,yrmed(288),yrssn(288)
	DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
	common /BL1/CC,XMED
C
      infile1='testmohm2014.out'
		OPEN(13,file=infile1,action='READ')
	read(13,*) txt ! extra line
C
	infile2='d:\web\dat6\YR\moyrmnh.txt' ! hmR2 data
	st='mo' 
	AYEAR='2014' 
      AYR='14'
	read(ayr,*) ryr
	iyr=int(ryr)
	if(int(ryr/4.)*4.eq.iyr) THEN 
 	im(2)=29
	                          ELSE
	im(2)=28
	ENDIF

      outfile='comparmedhm.res'
	OPEN(10,file=outfile,access='APPEND')
	outmed='stYRmedhm.txt'
	outmed(1:2)=st
	outmed(3:4)=AYR
C
      GOTO 90  ! Annual RMSE
      	OPEN(11,file=outmed,access='APPEND')
	write(*,32) AYEAR,st,' hmF2med'
	write(11,32) AYEAR,st,' hmF2med'
C
   30 FORMAT(A6,2X,24(F3.0,2X)) ! hmF2
   31 FORMAT(1X,2A2,1X,I2,2(2X,F6.1)) ! hmF2
   32 FORMAT(1X,A4,1X,A2,A8) ! hmF2med
   33 format(79X,F6.1) ! hmF2
        infile2(13:14)=AYR
	infile2(18:19)=AYR
C
	DO 70 imn=1,12
	call blet2(imn,AMN)
      infile2(20:21)=AMN
		OPEN(12,file=infile2,action='READ')
      do i=1,4
	read(12,*) txt ! extra line
	enddo
	ndmn=im(imn)
	DO kk=1,ndmn
	read(12,30,end=50) YYMMDD,(CC(kk,k),k=1,24) 
	ENDDO
   50	close(unit=12)
C
       call submedmon(ndmn,24)
      xav=0.
        do jj=1,24
	xav=xav+xmed(jj)
	 enddo
	xav=xav/24.
C
C Read Model results for month:
       DO jj=1,24
	read(13,33,end=51) ymod(jj)
	enddo
C
C RMSE  and PE
       dev=0.
	zet=0.
	DO jj=1,24
	dev=dev+(xmed(jj)-ymod(jj))**2
	zet=zet+(xmed(jj)-xav)**2
	ENDDO
	RMSE=sqrt(dev/24.)
	PE=(1.-sqrt(dev/zet))*100. ! %
C
        do iut=0,23
	write(*,31) AYR,AMN,iut,xmed(iut+1),ymod(iut+1)
      write(11,31) AYR,AMN,iut,xmed(iut+1),ymod(iut+1)
	 enddo
        write(*,34) AYR,AMN,dev,zet,RMSE,PE
	write(10,34) AYR,AMN,dev,zet,RMSE,PE
   70 CONTINUE ! next month
   51   close(unit=13)
      close(unit=11)
	close(unit=10)
C	Annual RMSE , PE
   90	infile1=outmed
      OPEN(10,file=outfile,access='APPEND')
		OPEN(12,file=infile1,action='READ')
	read(12,*) txt
	write(*,*) txt
	write(10,*) txt
C
	yrav=0.
	do kk=1,288
	read(12,31) AYR,AMN,iut,yrmed(kk),yrssn(kk)
	yrav=yrav+yrmed(kk)
	enddo
	yrav=yrav/288.
C
       dev=0.
	zet=0.
	DO jj=1,288
	dev=dev+(yrmed(jj)-yrssn(jj))**2
	zet=zet+(yrmed(jj)-yrav)**2
	ENDDO
	RMSE=sqrt(dev/288.)
	PE=(1.-sqrt(dev/zet))*100. ! %
        write(*,35) AYR,AMN,dev,zet,RMSE,PE
	write(10,35) AYR,AMN,dev,zet,RMSE,PE
	close(unit=10)
		close(unit=12)
C
   34 format(1X,2A2,2(1X,F8.2),1X,F7.2,1X,F8.3) ! foF2
   35 format(1X,2A2,2(1X,F10.2),1X,F7.2,1X,F8.3) ! foF2
       pause ' '
	STOP
	END