program statsbreaker
c	This is a program that should extract the forecasts in linear array to calculate
c	a 24 hour running statistical count for mean, standard deviation 
c	for each station entered.
c
c
      real forecast(4416), ztime(184,24), hrmean(24), zsum, stdev(24)
      real zdev, zseq(184, 24), a1(184), temp , q1(24), q2(24), q3(24)
      CHARACTER aname*50
      CHARACTER  aname4*67
      CHARACTER  stnm3*5
	integer icount, leastindex, szz, ja, ka , ib, ic, id
      integer stnm

c	Initializing the forecast array with -9999 flags
	DO i=1,4416,1
	forecast(i)=-9999
	ENDDO
c	read inthe station WMO number
c	WRITE(6,*) 'Input stnm   '
c	READ(5,'(i5)') stnm
      OPEN(UNIT=13, file='goodstns.txt')
      DO ia=1,34
      READ(13,'(i5)') stnm
      write (stnm3,'(i5)') stnm


c	now I concatonate the files to make a path and file name
      aname='/snowball2/s0/apaget/pblmodel/1DPBL/results/'//stnm3(1:5)
	aname4=aname(1:49)//'/forecast'//stnm3(1:5)//'.txt'
c	now I open and read in the ascii forecast
      OPEN(UNIT=16, FILE=aname4, STATUS='old')
	DO i=1,4416
	READ(16, *) forecast(i)
	ENDDO
	CLOSE(16)
c	write everything to one array
	DO i=1,184
		DO j=1,24
			ztime(i,j)=forecast((i-1)*24+j)
		ENDDO
	ENDDO

c	just a test print to make sure the file is reading right
c	DO j=1,24
c	print*, ztime(1,j)
c	ENDDO

c	now I calculate the mean
	DO j=1,24
		icount=0
		zsum=0.
		DO i=1,184
			IF(ztime(i,j).GT.0) THEN
					zsum=zsum+ztime(i,j)
					icount=icount+1
			ENDIF
		ENDDO
		IF(icount.ne.0) hrmean(j)= zsum/icount
c	print*, icount
	ENDDO
c 	now I calculate the stdev
	DO j=1,24
		icount=0
		zdev=0.
		DO i=1,184
			IF(ztime(i,j).GT.0) THEN
					zdev=zdev+(ztime(i,j)-hrmean(j))**2
					icount=icount+1
			ENDIF
		ENDDO
		IF(icount.ne.0) stdev(j)= SQRT(zdev/icount)
c	print*, icount
	ENDDO

c	initialize the zsequence array with flags
	DO j=1,24,1
		DO i=1,184
			zseq(i,j)=-9999
		ENDDO
	ENDDO
c 	now I sequence the data from least to greatest
	DO j=1,24,1
c		j=1
		icount=1
		DO i=1,184
			IF(ztime(i,j).GT.0.0) THEN
				a1(icount)=ztime(i,j)
				icount=icount+1
			ENDIF
		ENDDO	
		szz=icount-1
		ja=0
		ka=0
		leastindex=0
		temp=0
		DO ja=1,szz
			leastindex=ja
			DO ka=ja+1,szz
				IF(a1(ka).lt.a1(leastindex)) THEN
					leastindex=ka
				ENDIF
			ENDDO
c	swap values in the array
			temp=a1(ja)
			a1(ja)=a1(leastindex)
			a1(leastindex)=temp
		ENDDO
		DO i=1,184,1
			zseq(i,j)=a1(i)
		ENDDO
	ENDDO

c 	calculate the quartiles
c 	I got this information from http://www.vias.org/tmdatanaleng/cc_quartile.html
	DO j=1,24,1
		icount=0
		DO i=1,184,1
			IF(zseq(i,j).GT.30.0) THEN
				icount=icount+1
			ENDIF
		ENDDO
c 	print*, icount
c	this will calculate the 1st quartile of obs.
		ic=ANINT(0.25*(icount+1))
		q1(j)=zseq(ic,j)
c		print*, 'ic is   ', ic
c	this will calculate the 2nd quartile if odd number of obs.
		IF(MOD(icount,2.).EQ.1) THEN
			ib=(icount+1)/2
			q2(j)=zseq(ib,j)
c		print*, 'ib is   ', ib
		ENDIF
c	this will calculate the 2nd quartile if even number of obs.
		IF(MOD(icount,2.).EQ.0) THEN
			ib=(icount)/2
			q2(j)=(zseq(ib,j)+zseq(ib+1,j))/2
c		print*, 'ib is   ', ib
		ENDIF
c	this will calculate the 3rd quartile of obs
		id=ANINT(0.75*icount)
		q3(j)=zseq(id,j)
c		print*, 'id is   ', id
		IF(ib.le.ic.or.ib.ge.id.or.ic.ge.id) print*, 'flag'
	ENDDO

c	just a test print to make sure the mean and stdev calculator is right
c	print*, stnm3
c	print*, 'hour   ', 'mean height    ', 'STDEV    ', 'Q1       ','Q2        ', 'Q3       ', 'Mean-Median'
c	DO j=1,24
c	print*, j-1, hrmean(j), stdev(j), q1(j), q2(j), q3(j), hrmean(j)-q2(j)
c	ENDDO

	OPEN(UNIT=23, FILE=stnm3//'.txt', STATUS='new')
	WRITE(23, '(i5)') stnm
	DO j=1,24
	WRITE(23, 403) j-1, hrmean(j), stdev(j), q1(j), q2(j), q3(j), hrmean(j)-q2(j)
  403 FORMAT(i2,3x,6F9.3)
	ENDDO

      ENDDO
c	DO i=1,184
c		IF(zseq(i,24).EQ.-9999) print*, i,zseq(i,24)
c		print*, zseq(46,24),  zseq(92,24), zseq(93,24), zseq(138,24)
c	ENDDO
      end