character*180 fleo1,fdir(6),foname,fmonth(20),fname,foname2
	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)

c summary for level 2 sub-basins
	integer basin2(126,221) !basin2 id file
	integer sub_id(100),nsubid !sub basin id
	integer sub_n(100) !valid data number
	real sub_ser(50,100,6)  !valid data time series (yrs,basins,case) 
	real sub_sum(100)  !sum of valiad data
	real trend(126,221),retrend(126,221) !results
	real subtr(100),subrtr(100) !results

	istay=1960
	iendy=2000

	nsy=3

	max_nr=126
	max_nc=221

	foname	='_series.txt'

	fdir(1)='Allchange_HI'
	fdir(2)='noclimate_HI'
	fdir(3)='novegetation_HI'
	fdir(4)='noirrarea_HI'
	fdir(5)='Nochange_HI'
c	fdir(6)='obv_year.txt'

	fmonth(1)='tgwir'
	fmonth(2)='wthdr_h.asc'
	fmonth(3)='dSiB.asc'
	fmonth(4)='ETmass_h.asc'
	fmonth(5)='fshort.asc'
	fmonth(6)='GWdep_h.asc'
	fmonth(7)='gwsoil_h.asc'
	fmonth(8)='raet.asc'
	fmonth(9)='raht.asc'
	fmonth(10)='retnfw.asc'
	fmonth(11)='tc.asc'
	fmonth(12)='td.asc'
	fmonth(13)='tg.asc'
	fmonth(14)='tprec_h.asc'
	fmonth(15)='wreq.asc'
	fmonth(16)='wthd.asc'
	fmonth(17)='assimn.asc'	
	fmonth(18)='roff_h.asc'



	fname='irrcode.asc'
	CALL readfile_int (fname,nr,nc,x0,y0,s,
     $	nodata,ircod,max_nr,max_nc)
	fname='iumd.asc'
	CALL readfile_int (fname,nr,nc,x0,y0,s,
     $	nodata,iumd,max_nr,max_nc)
	fname='tangnaihai_loca.asc'
	CALL readfile_int (fname,nr,nc,x0,y0,s,
     $	nodata,itang,max_nr,max_nc)
	fname='wsdem.asc'
	CALL readfile_float (fname,nr,nc,x0,y0,s,
     $	fnodata,fdem,max_nr,max_nc)
	fname='irriratio.asc'
	CALL readfile_float (fname,nr,nc,x0,y0,s,
     $	fnodata,firatio,max_nr,max_nc)
	fname='basin2.asc'
	CALL readfile_int (fname,nr,nc,x0,y0,s,
     $	nodata,basin2,max_nr,max_nc)
	fname='basin2_list'
	open(1,file=fname,status='old')
	do i=1,100
	   read(1,*,end=20) sub_id(i)
	enddo
 20	close(1)
	nsubid = i-1
	sub_n=0

	do j=1,18
	sub_ser=0.  !valid data time series (yrs,basins,case) 
	subtr=0.
	subrtr=0.		!results
	foname2='HI'//trim(fmonth(j))//trim(foname)
	open(10,file=foname2,status='unknown')
	ouarr =0.
	do i=1,5
	do iy=istay,iendy
	do im=1,12
	Call ConvMMYYYY(im,iy,MMYYYY)
	fleo1='../'//trim(fdir(i))//'/HYDRO_OUTPUT/'
     $	//MMYYYY//'_'//trim(fmonth(j))
	CALL readfile_float (fleo1,nr,nc,x0,y0,s,
     $	fnodata,xdata,max_nr,max_nc)
	if(trim(fmonth(j)).eq.'assimn.asc') then
	do inr=1,max_nr
	  do inc=1,max_nc
		if (xdata(inr,inc).ne.-9999.) then
			xdata(inr,inc)=xdata(inr,inc)*(10.**8)
		else
			xdata(inr,inc)=-9999.
		endif
		enddo
	enddo
	endif

	ftotv = 0.0
	ntotv = 0
	firrv = 0.0
	nirrv = 0
	fhi   = 0.0
	nhi   = 0
	fhetao=0.0
	nhetao=0

	fup=0.0
	nup=0
	fmid=0.0
	nmid=0
	fdown=0.0
	ndown=0
	ftang=0.0
	ntang=0


	do inr=1,max_nr
	  do inc=1,max_nc
	     if (fdem(inr,inc).gt.-1.) then
				if (xdata(inr,inc).ne.-9999.) then
		ftotv = ftotv + xdata(inr,inc)
		ntotv = ntotv + 1
				endif
	     endif
	     if (ircod(inr,inc).gt.-1) then
				if (xdata(inr,inc).ne.-9999.) then
		firrv = firrv + xdata(inr,inc)
		nirrv = nirrv + 1
				endif
	     endif
	     if (ircod(inr,inc).eq.55) then
				if (xdata(inr,inc).ne.-9999.) then
		fhetao = fhetao + xdata(inr,inc)
		nhetao = nhetao + 1
				endif
	     endif
	     if  (fdem(inr,inc).gt.-1.) then
	     if (firatio(inr,inc).gt.30.) then
				if (xdata(inr,inc).ne.-9999.) then
		fhi = fhi + xdata(inr,inc)
		nhi = nhi + 1
				endif
	     endif
	     endif

	     if (fdem(inr,inc).gt.-1.) then
			 if (iumd(inr,inc).eq.1) then
				if (xdata(inr,inc).ne.-9999.) then
		fup = fup + xdata(inr,inc)
		nup = nup + 1
				endif
			 endif
	     endif
	     if (fdem(inr,inc).gt.-1.) then
			 if (iumd(inr,inc).eq.2) then
				if (xdata(inr,inc).ne.-9999.) then
		fmid = fmid + xdata(inr,inc)
		nmid = nmid + 1
				endif
			 endif
	     endif
	     if (fdem(inr,inc).gt.-1.) then
			 if (iumd(inr,inc).eq.3) then
				if (xdata(inr,inc).ne.-9999.) then
		fdown = fdown + xdata(inr,inc)
		ndown = ndown + 1
				endif
			 endif
	     endif
	     if (fdem(inr,inc).gt.-1.) then
			 if (itang(inr,inc).gt.0) then
				if (xdata(inr,inc).ne.-9999.) then
		ftang = ftang + xdata(inr,inc)
		ntang = ntang + 1
				endif
			 endif
	     endif

	  enddo
	enddo
	ftotv = ftotv / ntotv
	firrv = firrv / nirrv
	fhi   = fhi / nhi
	fhetao= fhetao / nhetao
	fup		= fup / nup
	fmid	= fmid / nmid 
	fdown	= fdown / ndown
	ftang	= ftang / ntang

	ii =iy -istay +1
	ouarr(ii,1,i)=ouarr(ii,1,i)+ftotv
	ouarr(ii,2,i)=ouarr(ii,2,i)+firrv
	ouarr(ii,3,i)=ouarr(ii,3,i)+fhi
	ouarr(ii,4,i)=ouarr(ii,4,i)+fhetao
	ouarr(ii,5,i)=ouarr(ii,5,i)+fup
	ouarr(ii,6,i)=ouarr(ii,6,i)+fmid
	ouarr(ii,7,i)=ouarr(ii,7,i)+fdown
	ouarr(ii,8,i)=ouarr(ii,8,i)+ftang

	sub_n=0
	sub_sum=0.
	do inr=1,max_nr
	  do inc=1,max_nc
	     if (fdem(inr,inc).gt.-1.) then
		if (xdata(inr,inc).ne.-9999.) then
		   ii=basin2(inr,inc)
		   if (ii.gt.0) then
		   sub_n(ii)=sub_n(ii)+1
		   sub_sum(ii)=sub_sum(ii)+xdata(inr,inc)
		   endif
		endif
	     endif
	  enddo
	enddo
	ii =iy -istay +1
	do jj=1,100
	   if(sub_n(jj).ne.0) then
	      sub_sum(jj)=sub_sum(jj)/sub_n(jj)
	      sub_ser(ii,jj,i)=sub_ser(ii,jj,i)+sub_sum(jj)
	   endif
	enddo
	
	enddo
	enddo	
	enddo

	ouarr=ouarr/12.
c	write(10,'(6(A10,1x))') 'iy','im','ftotv','firrv','fhi','fhetao'

	do k=1,8
	do i=1,5
		do iy=istay,iendy
			ii=iy -istay +1
			Xj(ii)=ii
			Fj(ii)=ouarr(ii,k,i)
		enddo
		N=ii
		call HPINT(N,Xj,Fj,a,b,r2,AVy,Sign)
		ouarr(ii+1,k,i)=a
		ouarr(ii+2,k,i)=b
		ouarr(ii+3,k,i)=r2
		ouarr(ii+4,k,i)=AVy
		ouarr(ii+5,k,i)=Sign
	enddo
	enddo

	sub_ser=sub_ser/12.
	do i=1,5
	subtr=0.
	subrtr=0.
	trend=0.
	retrend=0.

	do k=1,100 !loop of sub basins
	   do iy=istay+nsy,iendy
	      ii=iy -(istay+nsy) +1
	      Xj(ii)=ii
	      if (i.eq.1) Fj(ii)=sub_ser(ii,k,i)
				if (i.ne.1) Fj(ii)=sub_ser(ii,k,1)-sub_ser(ii,k,i)
	   enddo
	   N=iendy-(istay+nsy) +1
	   call HPINT(N,Xj,Fj,a,b,r2,AVy,Sign)
	   ftr=b*N
	   fretr=b*N/AVy
	   do jj=1,nsubid
	      if(sub_id(jj).eq.k) then
		 subtr(jj)=ftr
		 subrtr(jj)=fretr
	      endif
	   enddo
	   do inr=1,max_nr
	      do inc=1,max_nc
		 if(basin2(inr,inc).eq.k)then
		 trend(inr,inc)=ftr
		 retrend(inr,inc)=fretr
		 endif
	      enddo
	   enddo
	enddo

	fname='./HIlv2rS1/'//trim(fdir(i))//trim(fmonth(j))//'_trS1.txt'
	open(71,file=fname,status='unknown')	
	do jj=1,nsubid
	   write(71,*) sub_id(jj),subtr(jj),subrtr(jj)
	enddo
	close(71)
	fname='./HIlv2rS1/'//trim(fdir(i))//trim(fmonth(j))//'_trendS1.asc'
	call writefile_float (fname,max_nr,max_nc,x0,y0,s,
     $	fnodata,trend,max_nr,max_nc)
	fname='./HIlv2rS1/'//trim(fdir(i))//trim(fmonth(j))//'_retrendS1.asc'
	call writefile_float (fname,max_nr,max_nc,x0,y0,s,
     $	fnodata,retrend,max_nr,max_nc)
	
	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