      program fcnrmsd
C.......................................................July 2015
C Normalized NRMSD Deviation of fcF2 maps at EQ-quiet and EQ-storm
C Moving window: = 24 grids around EQ epicenter
C Select area +-10 deg: 4 lats>glat,+4lats<glat; 2longs>glong+2 longs <glong
C
C Regarding fc-map and std for Geom-qdays and Wp-qdays.  
	dimension fx(71,0:72),fy(71,0:72),IM(12)
	character*64 qdays,eqlist,infile1,infile2,outfile
	character*33 txt
	character*4 AYEAR,CYEAR(17)
	character*2	AYR,AMN,ADY,AUT,CMN(12),QDY,BYR,BMN,BDY
	+,qmag(17,12),qion(17,12)
	character*1 sq,mi
	real ut1
	DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/ 
C
	qdays='qdays.txt'	 ! Q-days
	OPEN(10,file=qdays,action='read')	

	read(10,31) txt ! extra line
	write(*,31) txt

	DO i=1,17
	do k=1,12
	read(10,32,end=111) CYEAR(i),CMN(k),qmag(i,k),qion(i,k)
	enddo
	ENDDO
  111	close(unit=10)

   31	format(A33)
   32	format(A4,A2,2(1X,A2))

   33	format(73(1X,F4.1))
   35	format(3X,3A2,1X,F5.2,1X,3(1X,F6.2),2(1X,F5.2),1X,F5.1)
   36	format(A1)
c   37	format(A13,5X,I3)
C
C
   98 write(*,*) ' Enter s or q for EQ-storm or EQ-quiet: '
	read(*,36) sq
C   	sq='s' ! storm EQ	 !TEMP
   99	write(*,*) ' Enter m or i for Q-mag or Q-ion days: '
	read(*,36) mi
C   	mi='m' ! Mag-Qday	 !TEMP

	outfile='EQnrmsdxx.txt'	  ! 
	outfile(8:8)=mi
	outfile(9:9)=sq
	OPEN(15,file=outfile)
	if (sq.eq.'s') then
	write(*,*) ' EQ storm'
	write(15,*) ' EQ storm'
	else
	write(*,*) ' EQ quiet'
	write(15,*) ' EQ quiet'
	endif

	if (mi.eq.'m') then
	write(*,*) ' Mag Qdays'
	write(15,*) ' Mag Qdays'
	else
	write(*,*) ' Iono Qdays'
	write(15,*) ' Iono Qdays'
	endif
C
  100	eqlist='exteq6s.txt'
	eqlist(7:7)=sq
	open(13,file=eqlist,action='READ')
	eqcnt=0
  101	read(13,35,end=27,err=26) BYR,BMN,BDY,UT1
 	+,glat,glon,depth,eqmag,hlt,xhi
	if (BYR.eq.'16') goto 27   ! end of EQs
C
	read(BYR,*) ryr
	iyr=int(ryr)
		if(int(ryr/4.)*4.eq.iyr) THEN 
		 im(2)=29
	ndnyr=366
	                          ELSE
	IM(2)=28
	ndnyr=365
	ENDIF
C
	read(BMN,*) rmn
	imn=int(rmn)
	read(BDY,*) rdy
	idy=int(rdy)
	IF (ut1.lt.0.) THEN
	kut1=0
	                ELSE
	ut2=ut1+0.005
	kut1=int(ut2)+1
C
	if (kut1.eq.24) then
	kut1=0
	idy=idy+1
	   if (idy.gt.im(imn)) then
	   idy=1
	   imn=imn+1
      	if (imn.gt.12) then
      	imn=1
      	iyr=iyr+1
      	endif
	   endif
      endif
	ENDIF
C
	call blet2(iyr,AYR)
	call blet2(imn,AMN)
	call blet2(idy,ADY)
	call blet2(kut1,AUT)
	if (iyr.eq.99) then 
	kyr=1
	else
	kyr=iyr+2
	endif
C
	if (mi.eq.'m') then
	QDY=qmag(kyr,imn)
	else
	QDY=qion(kyr,imn)
	endif
C
	if (QDY.eq.ADY) GOTO 101 ! Avoid rms calculation
	AYEAR=CYEAR(kyr)
C
	infile1='d:\web\dfs\YR\MN\fcMNDDUT.jYR' ! Q-day
C
	infile1(12:13)=AYR
	infile1(15:16)=AMN
	infile1(20:21)=AMN
	infile1(28:29)=AYR
C
	infile1(22:23)=QDY
	infile1(24:25)=AUT
	OPEN(11,FILE=infile1,action='READ')
	infile2=infile1
	infile2(22:23)=ADY

	OPEN(12,FILE=infile2,action='READ') ! EQ day 
	    do ilat=1,71  !read  for EQ-day
	read(12,33,end=233,err=26) (fy(ilat,k),k=0,72)
	enddo
  233	close(unit=12)  !

	qlat=-90.0
	    do ilat=1,71  !read  for q-day
	qlat=qlat+2.5
	read(11,33,end=133,err=26) (fx(ilat,k),k=0,72)
	if ((glat.ge.qlat).and.(glat.lt.qlat+2.5)) iqlat=ilat
	enddo
  133	close(unit=11)  !
C 
	qlon=-185.
	xlon=glon
	if (glon.gt.180.) xlon=glon-360.
	do ilon=0,71
	qlon=qlon+5.
	if (qlon.ge.180.) qlon=qlon-360.

	if ((xlon.ge.qlon).and.(xlon.lt.qlon+5.0)) iqlon=ilon
	enddo
	iqlon_1=iqlon-1
	if (iqlon_1.lt.0) iqlon_1=iqlon_1+70
	iqlon1=iqlon+1
	if (iqlon1.gt.72) iqlon1=iqlon1-72
	iqlon2=iqlon+2
	if (iqlon2.gt.72) iqlon2=iqlon2-72

	fxav=0.
	cnt=0.
      	do ilat=iqlat-3,iqlat+4
	if ((ilat.ge.iqlat-1).and.(ilat.lt.iqlat+3)) then
c	do ilon=iqlon-1,iqlon+2
      fxav=fxav+fx(ilat,iqlon_1)+fx(ilat,iqlon)+fx(ilat,iqlon1)
     $+fx(ilat,iqlon2)
	cnt=cnt+4.
c	enddo
	               else
c	do ilon=iqlon,iqlon+1
      fxav=fxav+fx(ilat,iqlon)+fx(ilat,iqlon1)
	cnt=cnt+2.
c	enddo
	endif
		enddo
	fav=fxav/cnt
C  std
	devx=0.
	devxy=0.
      	do ilat=ilat-3,ilat+4
	    if ((ilat.ge.iqlat-1).and.(ilat.lt.iqlat+3)) then
c	do ilon=iqlon-1,iqlon+2
	devx=devx+(fav-fx(ilat,iqlon_1))**2.+(fav-fx(ilat,iqlon))**2.
	&+(fav-fx(ilat,iqlon1))**2.+(fav-fx(ilat,iqlon2))**2.
	devxy=devxy+(fx(ilat,iqlon_1)-fy(ilat,iqlon_1))**2.
	&+(fx(ilat,iqlon)-fy(ilat,iqlon))**2.
	&+(fx(ilat,iqlon1)-fy(ilat,iqlon1))**2.
	&+(fx(ilat,iqlon2)-fy(ilat,iqlon2))**2.
C	enddo
                                else
C	do ilon=iqlon,iqlon+1
	devx=devx+(fav-fx(ilat,iqlon))**2.+(fav-fx(ilat,iqlon1))**2.
	devxy=devxy+(fx(ilat,iqlon)-fy(ilat,iqlon))**2.
	&+(fx(ilat,iqlon1)-fy(ilat,iqlon1))**2.
C	enddo
                   endif
	enddo
	stdx2=devx/cnt
		stdx=0.
	if (stdx2.gt.0.) stdx=sqrt(stdx2)
	rms=0.0
	if (devxy.gt.0.0) rms=sqrt(devxy/cnt)
	rmsn=0.0
	if (stdx.gt.0.) rmsn=rms/stdx
C Output results:
      write(*,34) AYEAR,BMN,BDY,UT1,glat,glon,depth,eqmag,hlt,xhi
	+,fav,stdx2,rms,rmsn
	write(15,34) AYEAR,BMN,BDY,UT1,glat,glon,depth,eqmag,hlt,xhi
	+,fav,stdx2,rms,rmsn
	GOTO 101 ! Next EQ

   27	close(unit=15)
C
   34	format(1X,A4,2A2,1X,F5.2,1X,3(1X,F6.2),2(1X,F5.2),1X,F5.1,1X,F5.2
     +,3(1X,F7.3))
	goto 400
   26	write(*,*) 'Error in input file=',eqlist
  400	continue ! 
	pause ' '
	stop
	END
C
