      program fcnrmsd4a
C Total rms for 8 years
C.......................................................Aug. 2015
C.........CBK....f-maps
C.......................................................Aug. 2015
C Limit EQs by EIA area [-20 +20], [90 190]
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.  
C+ Vsigma Q days ranked for 27-prec days
C
C+ Calculate total average RMS, NRMS
C
	dimension fx(71,0:72),fy(71,0:72),IM(12)
	character*64 eqlist,infile1,infile2,outfile,txt
	character*4 AYEAR
	character*2	AYR,AMN,ADY,AUT,BYR,BMN,BDY,QYR,QMN,QDY
	character*1 sq
	CHARACTER*6 YYMMDDMAX,YYMMDDMIN
	REAL ut1
	DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/ 
C
   31	format(A33)
   33	format(73(1X,F4.2))
   35	format(1X,A4,2A2,1X,F5.2,1X,3(1X,F6.2),2(1X,F5.2),1X,F5.1)
   36	format(A1)
C
	ivcnt=0
	iyear0=1998
C
   98 write(*,*) ' Enter s or q for EQ-storm or EQ-quiet: '
	read(*,36) sq

C	outfile='EQnrmsdvx.txt'	  ! 
      outfile='EQfcaver4.txt'	  !
c	outfile(9:9)=sq
	OPEN(15,file=outfile,access='APPEND')
	if (sq.eq.'s') then
	write(*,*) ' EQ storm'
	write(15,*) ' EQ storm'
	else
	write(*,*) ' EQ quiet'
	write(15,*) ' EQ quiet'
	endif
      	write(*,*) ' Vs rank Qdays' ! Vs  TEMP
		write(15,*) ' Vs rank Qdays'  !Vs TEMP
C?		write(*,*) ' Mag Qdays'  ! TEMP
C?	write(15,*) ' Mag Qdays'  !	TEMP
C
C  100	eqlist='exteq6s.txt'
  100	eqlist='EQnrmsdvx.txt'	  !
	eqlist(9:9)=sq
C for Mag-Q-days:
C?	eqlist(8:8)='m' !TEMP
C
	open(13,file=eqlist,action='READ')
	read(13,31) txt !extra line
	read(13,31) txt !extra line
Crem  103	eqcnt=0
      eqcnt=0.
	faveq=0.
	stdx2eq=0.
	rmseq=0.
	rmsneq=0.
	devxeq=0.
	devxyeq=0.
	aveq=0.
103	iyear0=iyear0+1
  101	read(13,34,end=27,err=26) AYEAR,BMN,BDY,UT1
 	+,glat,glon,depth,eqmag,hlt,xhi,fav,stdx2,rms,rmsn,YYMMDDMIN,AUT
	if (AYEAR.eq.'2007') goto 27
C++ only EIA area:
	if ((glat.lt.-20.).or.(glat.gt.20.)) goto 101
	if ((glon.lt.90.).or.(glon.gt.190.)) goto 101
CREM+	if (xhi.gt.90.) goto 101  ! only day
CREM+	if (xhi.le.90.) goto 101  ! only night
C Season
C EQN
CREM+	IF((bmn.eq.'03').or.(bmn.eq.'04').or.(bmn.eq.'09').or.
CREM+	+(bmn.eq.'10')) then
CREM+      goto 105 !EQN
CREM+	else
CREM+	goto 101
CREM+	endif

C 	
!WNT
CREM+	IF (glat.ge.0) then 
CREM+	 IF ((bmn.eq.'01').or.(bmn.eq.'02').or.(bmn.eq.'11')
CREM+	+ .or.(bmn.eq.'12')) then    !North WNT
CREM+	goto 105
CREM+	 else
CREM+	goto 101
CREM+	 endif
CREM+	                else
CREM+	  IF ((bmn.eq.'05').or.(bmn.eq.'06').or.(bmn.eq.'07')
CREM+	+  .or.(bmn.eq.'08')) then	  !South WNT 
CREM+	 goto 105 !
CREM+	else
CREM+	goto 101
CREM+	endif
CREM+	ENDIF

C	 SMR
	IF (glat.lt.0) then 
	 IF ((bmn.eq.'01').or.(bmn.eq.'02').or.(bmn.eq.'11')
	+ .or.(bmn.eq.'12')) then    !South SMR
	goto 105
	 else
	goto 101
	 endif
	                else
	  IF ((bmn.eq.'05').or.(bmn.eq.'06').or.(bmn.eq.'07')
     +  .or.(bmn.eq.'08')) then	  !North SMT 
	 goto 105 !
	else
	goto 101
	endif
	ENDIF
C++ EIA:
  105	BYR=AYEAR(3:4)
	read(BMN,*) rmn
	imn=int(rmn)

	read(AYEAR,*) ryear
	iyyyy=int(ryear)
	if (iyyyy.eq.iyear0) then 
	goto 102
	else
	backspace(13)
	goto 27
	endif
C	if (BYR.eq.'16') goto 27   ! end of EQs
	if (AYEAR.eq.'2007') goto 27   ! end of EQs
CREM	if (xhi.gt.90.) then  !exclude nighttime
CREM	if (xhi.lt.90.) then  !exclude daytime
CC	IF ((imn.eq.3).or.(imn.eq.4).or.(imn.eq.9).or.(imn.eq.10)) THEN !EQN
CCC	IF ((glat.ge.0).and.((imn.lt.3).or.(imn.gt.10))) THEN !North WNT
CCC	IF ((glat.lt.0).and.((imn.lt.3).or.(imn.gt.10))) THEN !South SMR
CCC	if (eqcnt.eq.0) then
CREM	write(*,*) 'Day'
CREM		write(15,*) 'Day'
Crem		write(*,*) 'Night'
Crem	write(15,*) 'Night'
Crem	endif
CREM      goto 101
CREM      endif
C
CC			write(*,*) 'Equinox' 
CC	write(15,*) 'Equinox' 
CCC	write(*,*) 'Winter'
CCC	write(15,*) 'Winter' 
CCC	write(*,*) 'Summer' 
CCC	write(15,*) 'Summer' 
CCC	endif
CCC	goto 102
CCC	                 ENDIF
CCC	IF ((glat.lt.0).and.((imn.gt.4).and.(imn.lt.9))) THEN !South WNT
CCC	IF ((glat.ge.0).and.((imn.gt.4).and.(imn.lt.9))) THEN !North SMR
CCC	if (eqcnt.eq.0) then
CCC	write(*,*) 'Winter'
CCC	write(15,*) 'Winter' 
CCC	write(*,*) 'Summer' 
CCC	write(15,*) 'Summer' 
CCC	endif
CCC	goto 102
CCC                                        ELSE
CCC	goto 101
CCC	                                  ENDIF
  102    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(BDY,*) rdy
	idy=int(rdy)
C
C Calculate rank-27 Qday:
	call subrankf(iyr,imn,idy,YYMMDDMAX,jmax,YYMMDDMIN,jmin)
	QYR=YYMMDDMIN(1:2)
	QMN=YYMMDDMIN(3:4)
	QDY=YYMMDDMIN(5:6)

	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)
C
CREM	infile1='d:\web\dfs\YR\MN\fcMNDDUT.jYR' ! Q-day
	infile1='d:\Mapf\20YR\MN\fcMNDDUT.fYR' ! Q-day
C
	infile1(11:12)=QYR
	if ((QYR.eq.'98').or.(QYR.eq.'99')) then
	infile1(9:10)='19'
	                     else
	infile1(9:10)='20'
	endif
	infile1(14:15)=QMN
	infile1(19:20)=QMN
	infile1(27:28)=QYR
C
	infile1(21:22)=QDY
	infile1(23:24)=AUT
	OPEN(11,FILE=infile1,action='READ')
	infile2=infile1
	infile2(9:12)=AYEAR
	infile2(14:15)=AMN
	infile2(19:20)=AMN
	infile2(27:28)=AYR
	infile2(21:22)=ADY

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

	qlat=-62.5
C	    do ilat=1,71  !read  for q-day
	    do ilat=12,60  !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=0.
	if (cnt.gt.0.) stdx2=devx/cnt
		stdx=0.
	if (stdx2.gt.0.) stdx=sqrt(stdx2)
	rms=0.0
	if ((devxy.gt.0.0).and.(cnt.gt.0.)) rms=sqrt(devxy/cnt)
	rmsn=0.0
	if (stdx.gt.0.) rmsn=rms/stdx
C Output results:
      aveq=aveq+fxav
	devxeq=devxeq+devx
	devxyeq=devxyeq+devxy
	eqcnt=eqcnt+1

c      write(*,34) AYEAR,BMN,BDY,UT1,glat,glon,depth,eqmag,hlt,xhi
c	+,fav,stdx2,rms,rmsn,YYMMDDMIN,AUT
c	write(15,34) AYEAR,BMN,BDY,UT1,glat,glon,depth,eqmag,hlt,xhi
c	+,fav,stdx2,rms,rmsn,YYMMDDMIN,AUT
	 ivcnt=ivcnt+1
	GOTO 101 ! Next EQ

   27	continue
c>      aveq=0.
		if (iyear0.lt.2007) goto 103
      if (eqcnt.gt.0.) aveq=aveq/(eqcnt*24.)
	stdx2=0.
 	 if (eqcnt.gt.0.) stdx2=devxeq/(eqcnt*24.)
	stdx=0.
	if(stdx2.gt.0.) stdx=sqrt(stdx2)
	rms=0.
	if (eqcnt.gt.0.) rms=sqrt(devxyeq/(eqcnt*24.))
	rmsn=0.
	if (stdx.gt.0.) rmsn=rms/stdx
c	write(*,*) iyear0  !Temp
c		write(15,*) iyear0  !Temp
	write(*,39) iyear0,aveq,stdx,rms,rmsn,eqcnt 
	write(15,39) iyear0,aveq,stdx,rms,rmsn,eqcnt

c	if (iyear0.lt.2007) goto 103
	write(15,*) ivcnt
      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),1X,A6,1X,A2)
   39	format(1X,I4,4(1X,F7.3),1X,F6.0)
	goto 400
   26	write(*,*) 'Error in input file=',eqlist
  400	continue ! 
	write(*,*) ivcnt
	pause ' '
	stop
	END
C
