character*180 fleo1,fdir(6),foname,fmonth(20),fname,fout,fpre	
	character*200 chavg,tit
	character*40 stnname(40)
	CHARACTER*6		MMYYYY
	real		ouarr(50,8,6)
	Integer Iarr(20),istay,iendy
	Integer daym(12) /31,28,31,30,31,30,31,31,30,31,30,31/
	real    fmaxd(126,221),fdem(126,221),firatio(126,221)
	integer imaxd(126,221),ircod(126,221),iumd(126,221),itang(126,221)
	real		xdata(126,221)

	real	AXj,AFj,Sfj,Sxj,rfxj,MSEj,MSEcj,MSSSj,rmse,rrmse
	real	em_t,ea_t,ep_t
c	real	Xj(7320),Fj(7320) !7320 for 20 year running !18300 for 50 yrs
	real	Xj(18300),Fj(18300)

	real	sumwith(20,50)

	istay=1960
	iendy=2000


	max_nr=126
	max_nc=221

	fpre	='Discharge_'

	fdir(1)='Allchange_HI'
	fdir(2)='noclimate_HI'
	fdir(3)='novegetation_HI'
	fdir(4)='noirrarea_HI'
	fdir(5)='Nochange_HI'
	fdir(6)='Noch_rept60_HI'
	num =6

	fmonth(1)='tangnaihai'
	fmonth(2)='lanzhou'
	fmonth(3)='qingtongxia'
	fmonth(4)='shizuishan'
	fmonth(5)='toudaoguai'
	fmonth(6)='longmen'
	fmonth(7)='sanmenxia'
	fmonth(8)='huayuankou'
	fmonth(9)='lijin'
	
	do i=1,num
	foname= trim(fpre)//trim(fdir(i))//'.txt'
	open(10,file=foname,status='unknown')
c		write(10,*) trim(fdir(i))
		sumwith=0.
		do j=1,9
		fname='../'//trim(fdir(i))//'/HYDRO_OUTPUT/'
     $		//trim(fmonth(j))//'_y.txt'
		open(1,file=fname,status='old')
		read(1,*)
			do iy=istay,iendy
				k=iy-istay+1
					read(1,*) i_y,Qy,Qy_obv
					sumwith(j,k)=Qy
					sumwith(j+10,k)=Qy_obv
					Xj(k)=Qy_obv
					Fj(k)=Qy
			enddo
			Ni=iendy-istay+1
			k=iendy-istay+1
			CALL MSSS_Murphy
     $(Ni,Xj,Fj,AFj,AXj,Sfj,Sxj,rfxj,MSEj,MSEcj,MSSSj,em_t,ea_t,ep_t)
			sumwith(j,k+1)=(AFj-AXj)/AXj*100
			sumwith(j,k+2)=sqrt(MSEj)
			sumwith(j,k+3)=sqrt(MSEj)/AXj
			sumwith(j,k+4)=MSSSj
		enddo
		
c		do j=9,2,-1
c			do k=1,iendy-istay+1
c				sumwith(j,k)=sumwith(j,k)-sumwith(j-1,k)
c			enddo
c		enddo

		write(10,'(A4,1x,20(A10,1x))') "YEAR",(trim(fmonth(j)),"OBV",j=1,9)
		do iy=istay,iendy
			k=iy-istay+1
			write(10,'(I4,1x,20(F10.3,1x))') 
     $			iy,(sumwith(j,k),sumwith(j+10,k),j=1,9)
		enddo
c		k=iendy-istay+1
c		do iy=1,4
c			write(10,'(I4,1x,20(F10.3,1x))') 
c     $			iy,(sumwith(j,k+iy),sumwith(j+10,k+iy),j=1,9)
c		enddo

	close(10)
	enddo


	end

c	The least-squares line uses a straight line Y = a + bX 
c	to approximate the given set of data (x1,y1), (x2, y2), ...
c	b: slope, a: intercept	
c	AVy: averaged value of Y
c	Sign: significance 
c	The statistical significance of the trends is evaluated using the 
c	Student¡¯s t-test with the following significance parameter:
c	 t = r \sqrt{(n-2)/(1-r^2)}
c	where r is the correlation coefficient between $t$ and $y$.
	Subroutine HPINT(N,X,Y,a,b,r2,AVy,Sign)
	DIMENSION X(*), Y(*)
	REAL	X, Y,slope,y_intercept,AVy,Sign
	INTEGER N
	REAL	a,b,r2,bbru,bbrd,bbr

	a=0.0
	b=0.0
	r2=0.0
	AVy=0.0
	Sign=0.0
	SUMx=0.0
	SUMy=0.0
	SUMxy=0.0
	SUMxx=0.0
	Sumyy=0.0

	do i = 1, n
		SUMx = SUMx + x(i)
		SUMy = SUMy + y(i)
		SUMxy = SUMxy + x(i)*y(i)
		SUMxx = SUMxx + x(i)*x(i)
		Sumyy = SUMyy + y(i)*y(i)
	end do
	
	AVy = SUMy / N

	slope = ( SUMx*SUMy - n*SUMxy ) / ( SUMx*SUMx - n*SUMxx )
	y_intercept = ( SUMy - slope*SUMx ) / n

c	SSE = 0.0
c	SST = Sumyy - SUMy**2 / n
c	do i = 1, n
c		SSE = SSE + (y_intercept+slope*X(i)-Y(i))**2
c	enddo

c	bbru = n * SUMxy - SUMx * SUMy
c	bbrd = sqrt ( (n*SUMxx-SUMx*SUMx) * (n*Sumyy - SUMy*SUMy)  )
	bbr  = (n * SUMxy - SUMx * SUMy)/
     $ sqrt ( (n*SUMxx-SUMx*SUMx) * (n*Sumyy - SUMy*SUMy)  )

	r2= bbr*bbr
	Sign=bbr*sqrt( (n-2)/(1-bbr*bbr))*100.
cccccccccccccccccc
c	AVy=Sign
cccccccccccccccccc
	b = slope
	a = y_intercept
	END

c Verfication scores, refer to Murphy, 1988
c recommendate by WMO
c also refer to Arora 1999
c xij,fij (i=1, 2, .. Ni) denote time series of observations and forecasts
c Axj, AFj: averages
c Sfj, Sxj: sqrt ( sample variances)
c rfxj:	  product moment correlations of the forecasts and observations
c MSEj:	  mean squared error of the forcasts
c MSEcj:	  mean squared error of climatology forecasts (Murphy, 1988)
c MSSSj:	  Verfication scores
c em_t,ea_t,ep_t: error associated with overall mean, amplitude, and phase
c 
	subroutine MSSS_Murphy
     $(Ni,Xj,Fj,AFj,AXj,Sfj,Sxj,rfxj,MSEj,MSEcj,MSSSj,em_t,ea_t,ep_t)
	integer	Ni,Na
	DIMENSION Xj(*), Fj(*), Xa(Ni), Fa(Ni)
	real	AXj,AFj,Sfj,Sxj,rfxj,MSEj,MSEcj,MSSSj
	real	em_t,ea_t,ep_t

	integer i,j
	real	sumx,sumf,sumxx,sumff,sumfx,sumfx2
	real	em2,ea2,ep2,etot2
	
	Na=0
	do i=1,Ni
		if (abs(Xj(i)+9999.).gt.0.1.and.abs(Fj(i)+9999.).gt.0.1) then
			Na=Na+1
			Xa(Na)=Xj(i)
			Fa(Na)=Fj(i)
		endif
	enddo

	if (Na.lt.3) then
		AFj	=-9999.
		AXj	=-9999.
		Sfj	=-9999.
		Sxj	=-9999.
		rfxj	=-9999.
		MSEj	=-9999.
		MSEcj	=-9999.
		MSSSj	=-9999.
		em_t	=-9999.
		ea_t	=-9999.
		ep_t	=-9999.
	else
		sumx=0.
		sumf=0.
		do i=1,Na
			sumx=sumx+Xa(i)
			sumy=sumy+Fa(i)
		enddo
		AXj = sumx / Na
		AFj = sumy / Na

		sumxx=0.
		sumff=0.
		do i=1,Na
			sumxx=sumxx+ ( Xa(i)-AXj )*( Xa(i)-AXj )
			sumff=sumff+ ( Fa(i)-AFj )*( Fa(i)-AFj )
		enddo
		Sxj = sqrt( sumxx / (Na - 1) )
		Sfj = sqrt( sumff / (Na - 1) )

		sumfx=0.
		do i=1,Na
			sumfx=sumfx+ ( Xa(i)-AXj )*( Fa(i)-AFj )
		enddo
		rfxj = sumfx / Na / (Sfj*Sxj)

		sumfx2=0.
		do i=1,Na
			sumfx2=sumfx2+ ( Fa(i)-Xa(i) )*( Fa(i)-Xa(i) )
		enddo
		MSEj = sumfx2 / Na

		MSEcj = sumxx / Na
		MSSSj = 1 - MSEj / MSEcj
		
		em2 = (AXj-AFj)*(AXj-AFj)
		ea2 = ( sqrt(sumff/Na)-sqrt(sumxx/Na) )**2.
		ep2 = 2.*sqrt(sumff/Na)*sqrt(sumxx/Na)
     $		*(1.-sumfx/Na/(sqrt(sumff/Na)*sqrt(sumxx/Na))  )
		etot2=sumfx2/ Na
		em_t	= em2 / etot2 * 100.
		ea_t	= ea2 / etot2 * 100.
		ep_t	= ep2 / etot2 * 100.	
	endif

	end !subroutine MSSS_Murphy


      subroutine strlen(str, l1,l2)
      character str*120
      integer i,l1,l2,k
      k=0
      do i = 1, 200
        if(k.eq.0 .and. str(i:i).NE.' ') then
          l1=i
          k=1
       elseif(k.eq.1.and. str(i:i).EQ.' '.or.ICHAR(str(i:i)).eq.0) then
          l2 = i-1
          return
        endif
      end do
      l2 = i
      return
      end

c	This subroutine convert Month/Year to MMYYYY format CHARACTER
	subroutine ConvMMYYYY(MM,YYYY,MMYYYY)
	implicit none
	integer MM,YYYY
	CHARACTER*6	MMYYYY

	if (MM.lt.1.or.MM.gt.12) then
		PRINT *, 'ERROR INPUT MONTH, MONTH SHOULD BETWEEN 1-12'
		STOP
	else if (MM.lt.10) then
		MMYYYY(1:1)='0'
		Write (MMYYYY(2:2),'(i1)') MM
	else 
		Write (MMYYYY(1:2),'(i2)') MM
	endif

	if (YYYY.lt.1000.or.YYYY.gt.9999) then
		PRINT *, 'ERROR INPUT YEAR, YEAR SHOULD BETWEEN 1000-9999'
		STOP
	else 
		Write (MMYYYY(3:6),'(i4)') YYYY
	endif

	end


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c	Read Write ascii file, the file format is ArcInfo ascii file format 
c	Input: 
c		name:			filename of the read/write file
c		nr:				nrows
c		nc:				ncols
c		nnr,nnc:		the size of input/output array
c		x0:				xllcorner	
c		y0:				yllcorner
c 		s:				cellsize
c		znodata:		NODATA_value
c		array(/write):	the array will be written into the file			
c	Output(/read):
c		array:			the file data will be read into the array
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	subroutine readfile_int(name,nr,nc,x0,y0,s,znodata,array,nnr,nnc)
	implicit none
      integer znodata,nc,nr
	integer nnr,nnc
	character*180  name
	integer array(nnr,nnc)
	character*14 ncols,nrows
      character*14 xllcorner,yllcorner
      character*14 cellsize
      character*14 nodata
      real x0,y0,s
	integer i,j

	open(1, file= name ,status='old')
      read(1,'(A, i)') ncols,nc
      read(1,'(A, i)') nrows,nr
      read(1,'(A, f)') xllcorner,x0
      read(1,'(A, f)') yllcorner,y0
      read(1,'(A, f15.0)') cellsize,s
      read(1,'(A, i)') nodata,znodata    
      do i=1,nr        
      	read(1,*,err=10) (array(i,j), j=1,nc)
		goto 11
10		array(i,j)=-9999.
11		continue
      end do          
      close (1)
	end

	subroutine writefile_int(name,nr,nc,x0,y0,s,znodata,array,nnr,nnc)
	implicit none
      integer znodata,nc,nr
	integer nnr,nnc
	character*180  name
	integer array(nnr,nnc)
      real*4 x0,y0,s
	integer i,j
	
      open(1,file=name,status='unknown')
      write(1,'(A, i)') 'ncols',nc
      write(1,'(A, i)') 'nrows',nr
      write(1,'(A, f15.3)') 'xllcorner',x0
      write(1,'(A, f15.3)') 'yllcorner',y0
      write(1,'(A, f)') 'cellsize',s
      write(1,'(A, i15)') 'NODATA_value',znodata
      do i=1,nr
      	do j=1,nc
      		write (1,'(i,$)') array(i,j)
      	end do
		write (1,*)
      end do
      close(1)
	end

	subroutine readfile_float(name,nr,nc,x0,y0,s,znodata,array,nnr,nnc)
	implicit none
      integer nc,nr
	integer nnr,nnc
	character*180  name
	real array(nnr,nnc),znodata
	character*5 ncols,nrows
      character*9 xllcorner,yllcorner
      character*14 cellsize
      character*12 nodata
      real x0,y0,s
	integer i,j

c	print *, name
	open(1, file= name ,status='old')
      read(1,'(A, i)') ncols,nc
      read(1,'(A, i)') nrows,nr
      read(1,'(A, f)') xllcorner,x0
      read(1,'(A, f)') yllcorner,y0
      read(1,'(A, f15.0)') cellsize,s
      read(1,'(A, f15.0)') nodata,znodata  
      do i=1,nr 
      	read(1,*) (array(i,j), j=1,nc)
      end do          
      close (1)
	end

	subroutine writefile_float(name,nr,nc,x0,y0,s,znodata,array,nnr,nnc)
	implicit none
      integer nc,nr
	integer nnr,nnc
	character*180  name
	real array(nnr,nnc), znodata
      real x0,y0,s
	integer i,j

      open(1,file=name,status='unknown')
      write(1,'(A, i)') 'ncols',nc
      write(1,'(A, i)') 'nrows',nr
      write(1,'(A, f15.3)') 'xllcorner',x0
      write(1,'(A, f15.3)') 'yllcorner',y0
      write(1,'(A, f)') 'cellsize',s
      write(1,'(A, f)') 'NODATA_value',znodata
      do i=1,nr
      	do j=1,nc
      		write (1,'(f,$)') array(i,j)
      	end do
		write (1,*)
      end do
      close(1)
	end