C........................................................Apr 2014
C Add output ht=h05t
C
C Extract (1) for noon (2) for midnight
C Add output hmF2obs (1,2,3)
C
C NRMSEF=RMSE/avf*100%;   NRMSET=RMSE/avt*100%;
C .......................................................Mar 2014
C RMSE (1) between foF2-obs and fcF2-fitTEC; (3)hmF2obs and hcF2-fitTEC, 
C RMSE (2) between foF2-obs and fcF2-CCIR, (4) hmF2obs and hcF2-CCIR  at given sites
c RMSE (5) between TEC-obs and TEC-fitTEC; (6)TECobs and TEC-CCIR, (7) TECobs - TEC-fit(foF2&hmF2)
C outputA:  (1),(2),(3),(4),(5),(6),(7),cnthrs
c 
C outputB:UT (if foF2obs>0.),xhi,W,tau1=fitTEC,tau2=CCIR,tau3=fitfohm,Hsc1,Hsc2,Hsc3
C
C 2 days in Sep. 2002, 2 days in Oct. 2003; 31 stations
C
C
C Input files stYRMNDY.out [ fcF2-fitTEC, fcF2-CCIR]
C Input files stYRMNDY.br stYRMNg.txt (fcF2// from TECgps)
C & styrf.br(from IDCE) 	
C
C <<< to replace only line(s) for given day(s)
C 
C Output: annual file stYRf.br for Conver & /dat8/10/styrmn.txt
C
      PROGRAM DERMSHF7
C---------------------------------------------------------------
C data for 24 UT=0,1,2,...23hUT daily
C
C
C*********TEST
	CHARACTER*128 OUTFILE2,infilest,infile1,infile2,infile3
cold	+,OUTFILE1
	character*2 ayr,amn,ST,sts(31),cyr,cmn,cdy
	+,ADY,stinp
	INTEGER*2 IDY,KWIND(0:23)
	dimension glats(31),glons(31)
	+,CC(7,0:23),XAV(0:23)
	character*1 fm
	CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
	CHARACTER*4 STAT,STATS(31)
	dimension im(12),fobs(0:23),tobs(0:23)	!.ins;.inp files
	+,hsc1(0:23),hsc2(0:23),hsc3(0:23),tau1(0:23),tau2(0:23),tau3(0:23)
	+,t1(0:23),t2(0:23),t3(0:23),xhi(0:23)
	+,f1(0:23),f2(0:23),f3(0:23)
	+,h1(0:23),h2(0:23),h3(0:23)
	+,ht1(0:23),ht2(0:23),ht3(0:23)
      COMMON /BLM/ STAT,GLAT,GLON
			COMMON /BL1/CC,XAV
	DATA  IM/31,28,31,30,31,30,31,31,30,31,30,31/
	DATA STS
	+/'at','ch','co','ds','eg','ga','gs','jr','kk','ld'
     +,'mg','mh','mm','mo','nq','ok','pa','ro','rv','sd'
     +,'sf','so','ss','tk','tq','tr','tz','vt','wk','wp'
     +,'yg'/
c YEAR MMDD UThr LThr XHI  SSN  COV  Kpm    L  Glati Glong Mlati Mlong MoDip  hmF2  foF2    NmF2   Nes      hpl    Nepl   ECbot  ECtop   ECpl   TEC    TAU  h05b  h05t   Hsc
c    1    2    3    4   5    6    7    8    9     10    11    12    13    14    15    16      17    18       19      20      21     22     23    24     25    26    27    28
c         x         x         x         x         x         x         x         x         x         x         x         x         x         x         x         x         x
c        10        20        30        40        50        60        70        80        90       100       110       120       130       140       150       160       170
c 2016  601   .0  1.6 115.  36.  91. 2.7  1.56  38.4  23.6  36.7 103.6  47.1 331.5  5.01  311057.  24179. 20200.    296.   1.87   9.29   3.29  14.45  464.5 274.2 481.4  265.0

C
	IRUN=1	 ! 1 for foF2, 2 for hmF2
c
	      CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT
C
C read stations coordinates
	infilest='inputs1.inp'
  76	format(1X,I2,1X,F5.1,1X,F5.1,1X,A4,1X,A2)
   77	format(30X,F5.1,1X,F5.1,40X,A4,9X,A2)
		    OPEN(109,FILE=INFILEST,action='READ')
	do i=1,31
	read(109,77,end=78) glat,glon,stat,stinp
	if (stinp.eq.sts(i)) then
	write(*,76) i,glat,glon,stat,stinp
	glats(i)=glat
	glons(i)=glon
	stats(i)=stat
	 endif
	enddo
   78	close(unit=109)

C
  201	      WRITE(*,*)' ENTER yr,mn :'
		READ(*,15) AYR,AMN
   15	FORMAT(A2)
  202	      WRITE(*,*)' ENTER DY :'
		READ(*,15) ADY
	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)
		read(ady,*) xdy
 	idy=int(xdy)


C start foF2, hmF2 or TEC
cold	outfile1='dermshf3.txt'
	outfile2='dermshf7.txt'
cold	OPEN(110,FILE=outfile1,access='APPEND')
	OPEN(111,FILE=outfile2,access='APPEND')
	fm='f' !Temp
	infile1='d:\web\graf\2002\styrf.br'	! foF2 -obs.
	infile1(15:16)=AYR
	infile1(20:21)=AYR
	infile2='styrmndy.inp'	! TEC -obs.
	infile2(3:4)=AYR
	infile2(5:6)=AMN
	infile2(7:8)=ADY
C
	infile3='styrmndy.out'
	infile3(3:4)=AYR
	infile3(5:6)=AMN
	infile3(7:8)=ADY

   75	format(3A2,2X,24(F3.1,2X))

  199	nds=1
	sumsf1=0.
	sumsf2=0.
	sumst1=0.
	sumst2=0.
	sumst3=0.
	cnttot=0.
	avftot=0.
	avttot=0.
C    
	DO 88 ist=1,31
	st=sts(ist)
	write(*,18) st,ayr,amn,ady
cold	write(110,18) st,ayr,amn,ady
	write(111,18) st,ayr,amn,ady
	infile1(18:19)=st
	infile2(1:2)=st
	infile3(1:2)=st
C Read foF2obs.
	OPEN(109,FILE=infile1,action='READ')
	do k=0,23
	xhi(k)=0.
	fobs(k)=0.
	tobs(k)=0.
	hsc1(k)=0.
	hsc2(k)=0.
	hsc3(k)=0.
	tau1(k)=0.
	tau2(k)=0.
	tau3(k)=0.
	f1(k)=0.
	f2(k)=0.
	f3(k)=0.
	h1(k)=0.
	h2(k)=0.
	h3(k)=0.
	ht1(k)=0.
	ht2(k)=0.
	ht3(k)=0.
	enddo
  198	read(109,75,end=197) CYR,CMN,CDY,(fobs(k),k=0,23)
	do k=0,23
	cc(7,k)=fobs(k)
	enddo
	if ((cyr.eq.ayr).and.(cmn.eq.amn).and.(cdy.eq.ady)) then
	goto 197
	else
	do k=0,23
	do j=2,7
	cc(j-1,k)=cc(j,k)
	enddo
	cc(7,k)=0.
	enddo
	goto 198
	endif
  197	close(unit=109)
	call rmedian7(7)
C
  170	format(74X,F6.2) 
	OPEN(109,FILE=infile2,action='READ')
	do k=0,23
	read(109,170) tobs(k)	   ! TECobs
	enddo
   	close(unit=109)

C Read foF2-tecfit, foF2-ccir, foF2obs
	OPEN(108,FILE=infile3,action='READ')
	read(108,79) txt  !extra line
	do k=0,23
	read(108,178) xhi(k),xmlat,xmlon,xmodip,h1(k),f1(k),t1(k)
	+,tau1(k),ht1(k),hsc1(k)
CTEMP	      if (xmlat.lt.60.) goto 83   ! Select MLAT Zone 3
CTEMP	      if (xmlat.gt.40.) goto 83   ! Select MLAT Zone 1
CTEMP 	      if ((xmlat.gt.60.).or.(xmlat.lt.40.)) goto 83   ! Select MLAT Zone 2
	enddo
	read(108,79) txt  !extra line
	do k=0,23
	read(108,178) xhi(k),xmlat,xmlon,xmodip,h2(k),f2(k),t2(k)
	+,tau2(k),ht2(k),hsc2(k)
	enddo
	read(108,79) txt  !extra line
	do k=0,23
	read(108,178,end=81) xhi(k),xmlat,xmlon,xmodip,h3(k),f3(k),t3(k)
	+,tau3(k),ht3(k),hsc3(k)
	enddo
   81	continue

C
C start cycle on stations
   79 format(A128)
  178 format(21X,F4.0,32X,3(1X,F5.1),1X,F5.1,1X,F5.2,54X,F6.2
     +,1X,F6.1,7X,F5.1,1X,F6.1)

c
	glat=glats(ist)
	glon=glons(ist)
	sumf1=0.
	sumf2=0.
	sumt1=0.
	sumt2=0.
	sumt3=0.
	cnt=0.
	rmsf1=0.
	rmsf2=0.
	rmst1=0.
	rmst2=0.
	rmst3=0.
	avf=0.
	avt=0.
	do 122 k=0,23
	XHIMAX=0.         ! midnight
	XHIMIN=100.
	do jj=0,23
	if (xhi(jj).gt.xhimax) then
	jmax=jj
	xhimax=xhi(jj)
	endif
	if (xhi(jj).lt.xhimin) then
	jmin=jj
	xhimin=xhi(jj)
	endif
	enddo
CTEMP      if (xhi(k).lt.90.) goto 122   ! Select night
CTEMP	if (k.ne.jmax) goto 122			 !midnight CTEMP
	 if (xhi(k).gt.90.) goto 122   ! Select day
CTEMP  if (k.ne.jmin) goto 122	  !noon	 CTEMP
	if (fobs(k).gt.0.) then
	avf=avf+fobs(k)
	avt=avt+tobs(k)
	sumf1=sumf1+(fobs(k)-f1(k))**2
	sumf2=sumf2+(fobs(k)-f2(k))**2
	sumt1=sumt1+(tobs(k)-t1(k))**2
	sumt2=sumt2+(tobs(k)-t2(k))**2
	sumt3=sumt3+(tobs(k)-t3(k))**2
	 cnt=cnt+1.
	xpar=fobs(k)**2
	xmed=xav(k)**2
	delog=0.
	call windex(xpar,xmed,delog,iw)
	KWIND(k)=iw
cold	write(*,19) k,xhi(k),iw,fobs(k),tobs(k)
cold	+,tau1(k),hsc1(k),tau2(k),hsc2(k),tau3(k),hsc3(k)
cold	write(111,19)  k,xhi(k),iw,fobs(k),tobs(k)
cold	+,tau1(k),hsc1(k),tau2(k),hsc2(k),tau3(k),hsc3(k)
	rsc1=(h1(k)+hsc1(k))/h1(k)
	rsc2=(h2(k)+hsc2(k))/h2(k)
	rsc3=(h3(k)+hsc3(k))/h3(k)
c+	write(*,119) k,xhi(k),iw,f1(k),h1(k),rsc1
c+	+,f2(k),h2(k),rsc2,f3(k),h3(k),rsc3
c+	write(111,119) k,xhi(k),iw,f1(k),h1(k),rsc1
c+	+,f2(k),h2(k),rsc2,f3(k),h3(k),rsc3
	write(*,119) k,xhi(k),iw,tobs(k),f1(k),h1(k),ht1(k),t1(k),rsc1
 	+,f2(k),h2(k),ht2(k),t2(k),rsc2,f3(k),h3(k),ht3(k),t3(k),rsc3
	write(111,119) k,xhi(k),iw,tobs(k),f1(k),h1(k),ht1(k),t1(k),rsc1
 	+,f2(k),h2(k),ht2(k),t2(k),rsc2,f3(k),h3(k),ht3(k),t3(k),rsc3

	endif
  122	continue
c+  119	format(1X,I2,1X,F4.0,1X,I2,3(1X,F6.2,1X,F5.1,1X,F7.3))
  119	format(1X,I2,1X,F4.0,1X,I2,1X,F6.2,3(1X,F6.2,1X,F5.1,1X,F5.1
     +,1X,F6.2,1X,F7.3))
	if (cnt.gt.0.) then
	 fav=avf/cnt
	tav=avt/cnt
	rmsf1=sqrt(sumf1/cnt)/fav*100.
	rmsf2=sqrt(sumf2/cnt)/fav*100.
	rmst1=sqrt(sumt1/cnt)/tav*100.
	rmst2=sqrt(sumt2/cnt)/tav*100.
	rmst3=sqrt(sumt3/cnt)/tav*100.
	endif
cold	write(*,17) glat,glon,xmlat,xmlon,xmodip,rmsf1,rmsf2
cold	+,rmst1,rmst2,rmst3,cnt
cold	write(110,17) glat,glon,xmlat,xmlon,xmodip,rmsf1,rmsf2
cold	+,rmst1,rmst2,rmst3,cnt
c
	sumsf1=sumsf1+sumf1
	sumsf2=sumsf2+sumf2
 	sumst1=sumst1+sumt1
	sumst2=sumst2+sumt2
	sumst3=sumst3+sumt3
	cntot=cntot+cnt
	avftot=avftot+avf
	avttot=avttot+avt
   83	close(unit=108)

   88	continue
c	
	 rec=0.
	ret=0.
	if (cntot.gt.0.) then
	faver=avftot/cntot
	taver=avttot/cntot
	ref1=sqrt(sumsf1/cntot)/faver*100.
	ref2=sqrt(sumsf2/cntot)/faver*100.
	ret1=sqrt(sumst1/cntot)/taver*100.
	ret2=sqrt(sumst2/cntot)/taver*100.
	ret3=sqrt(sumst3/cntot)/taver*100.
	endif
cold	write(*,16) ref1,ref2,ret1,ret2,ret3,cntot
cold	write(110,16) ref1,ref2,ret1,ret2,ret3,cntot
cold	write(110,115)
	write(111,115)
C
c ***************************************************************
C
  115	format('%--------------------------------------------')
   16	format('%',5(1X,F5.1),2X,F4.0)
   17 format(1X,5(1X,F5.1),5(1X,F7.2),1X,F4.0)
   18	format('%',4(1X,A2))
   19	format(1X,I2,1X,F4.0,1X,I2,2(1X,F6.2),6(1X,F6.1))

C
C++
C output 
C
C
cold	close(unit=110)
	close(unit=111)
  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
      write(*,*) outfile2,infile1,infile2
      pause ' ' 
   34   STOP
       END

C***********************************************************
      subroutine blet2(ilet,alet2)
c     nn=1,2,3,4 

	CHARACTER*1 IN(0:9)
	INTEGER*2 ilet
	DATA IN/'0','1','2','3','4','5','6','7','8','9'/
	CHARACTER*2 alet2
CLjuba	a2='00'
	alet2='00'             !NEW+7
	j10=ilet/10
	j1=ilet-j10*10
    5	do i=0,9
	if (j10.eq.i) then
	alet2(1:1)=IN(i)
	endif
	enddo
   7	do i=0,9 
     	if (j1.eq.i) then
 	alet2(2:2)=IN(i)
	endif
	enddo 
	return
	end
C
      integer function ndoy(iyr,imn,idy)
C To define day-of-year
C
	       DIMENSION IM(12)
	integer*2 iyr,imn,idy
	DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
	z1=iyr/4.0
      jz1=int(z1)*4
      IF(jz1.EQ.iyr) THEN
               IM(2)=29
        ELSE
                IM(2)=28
	       ENDIF
C Day-of-year LDA1
		mosum=0
      if(imn.gt.1) then
         do 1234 i=1,imn-1
1234    mosum=mosum+im(i)
         endif
      ndoy=mosum+idy	
	END FUNCTION
C ===============================================================
