implicit none
	
integer, parameter :: Nyear=30, startyear=1981, Nstation=64, &
                      ift2m=100,iftmin=101,iftmax=102,ifprep=103,of=104    ! File id

real, parameter :: Undef = -99.0

! Each station can be calculated seperately, no need to include the station index
real, dimension(Nyear,12,31) :: T2m, Tmax, Tmin, Prep
real, dimension(12,2000)  ::  Tmin_tot, Tmax_tot, Prep_tot  ! The whole monthly data
integer, dimension(12) :: N_total

real, dimension(12) :: T2m_cli,Tmin_cli,Tmax_cli,  &    !(1)
                       Rtot_Cli	, &						!(2)
					   Rx1day_Cli , Rx5day_Cli ,   &	!(3)
					   NH35_Cli,			       &	!(4)
					   NC15_Cli, NC13_Cli,         & 	!(5)
					   NR50_Cli,				   &    !(6)
					   TX90p,TN10p,R95p,R99p,      &    !(7)
					   NTX90p_Cli, NTN10p_Cli,     &    !(8)
					   NR95p_Cli, NR99p_Cli				!(9)
					   	   
real, dimension(12) :: threshold


integer :: ista, im, iy,id, year
real :: dummy
character(LEN=20) :: stname


open(ift2m,file="data/T2m_OUT.txt",status='old')
open(iftmin,file="data/T2min_OUT.txt",status='old')
open(iftmax,file="data/T2max_OUT.txt",status='old')
open(ifprep,file="data/Tpr_OUT.txt",status='old')

open(of,file="data/output.txt",status='unknown')


!Loop for each station
Do ista=1,Nstation
	 !Read T2m
	write(*,*)ista,'/',Nstation
	call read_a_station(ift2m,T2m)
	call read_a_station(iftmin,Tmin)
	call read_a_station(iftmax,Tmax)
	call read_a_station(ifPrep,Prep)


	write(of,*)stname
	write(of,'(A10,12i10)') 'Thang',1,2,3,4,5,6,7,8,9,10,11,12
	write(of,'(A10,12A10)') '------','------','------','------','------','------','------','------','------','------','------','------','------'

	!(1)			   
	call calc_cli(T2m,T2m_cli)
	call calc_cli(Tmin,Tmin_cli)
	call calc_cli(Tmax,Tmax_cli)
    write(of,'(A10,12f10.3)')'T2m_Cli',T2m_cli(:)
    write(of,'(A10,12f10.3)')'Tmin_Cli',Tmin_cli(:)
    write(of,'(A10,12f10.3)')'Tmax_Cli',Tmax_cli(:)	 
	!(2)
	call calc_cli_tot(Prep,Rtot_Cli)
    write(of,'(A10,12f10.3)')'Rtot_Cli',Rtot_Cli(:)
	!(3)
	call calc_cli_xday(Prep,Rx1day_Cli,1)
	call calc_cli_xday(Prep,Rx5day_Cli,5)
    write(of,'(A10,12f10.3)')'Rx1day_Cli',Rx1day_Cli(:)
    write(of,'(A10,12f10.3)')'Rx5day_Cli',Rx5day_Cli(:)
	!(4)
	threshold=35.
	call calc_ge_cli(Tmax,threshold,NH35_Cli,0)
    write(of,'(A10,12f10.3)')'NH35_Cli',NH35_Cli(:)
	!(5)
	threshold=15.
	call calc_ge_cli(Tmin,threshold,NC15_Cli,1)
	threshold=13.	
	call calc_ge_cli(Tmin,threshold,NC13_Cli,1)
		
    write(of,'(A10,12f10.3)')'NC15_Cli',NC15_Cli(:)
	write(of,'(A10,12f10.3)')'NC13_Cli',NC13_Cli(:)

	!(6)
	threshold=50.
 	call calc_ge_cli(Prep,threshold,NR50_Cli,0)
    write(of,'(A10,12f10.3)')'NR50_Cli',NR50_Cli(:)

    !(7)
	call total_series(Tmax,Tmax_tot,N_total)
	call calc_mon_percentile(Tmax_tot,N_total,90.,TX90p);
	write(of,'(A10,12f10.3)')'TX90p',TX90p(:)

 	call total_series(Tmin,Tmin_tot,N_total)
	call calc_mon_percentile(Tmin_tot,N_total,10.,TN10p);
	write(of,'(A10,12f10.3)')'TN10p',TN10p(:)

 	call total_series(Prep,Prep_tot,N_total)
	call calc_mon_percentile(Prep_tot,N_total,95.,R95p);
	call calc_mon_percentile(Prep_tot,N_total,99.,R99p);

	write(of,'(A10,12f10.3)')'R95p',R95p(:)
  	write(of,'(A10,12f10.3)')'R99p',R99p(:)

	!(8)
 	call calc_ge_cli(Tmax,TX90p,NTX90p_Cli,0)
 	call calc_ge_cli(Tmin,TN10p,NTN10p_Cli,1)
	write(of,'(A10,12f10.3)')'NTX90p_Cli',NTX90p_Cli(:)
  	write(of,'(A10,12f10.3)')'NTN10p_Cli',NTN10p_Cli(:)
	
	!(9)
 	call calc_ge_cli(Prep,R95p,NR95p_Cli,0)
 	call calc_ge_cli(Prep,R99p,NR99p_Cli,0)
	write(of,'(A10,12f10.3)')'NR95p_Cli',NR95p_Cli(:)
  	write(of,'(A10,12f10.3)')'NR99p_Cli',NR99p_Cli(:)


Enddo




Stop

!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
Contains

subroutine read_a_station(ifile,field)
	real, dimension(Nyear,12,31) :: field
	integer :: year, month,ifile

	read(ifile,*) stname
	
	do iy=1,Nyear
	   read(ifile,*)	year
	   do id=1,31
		  !  1   10.838    9.973   14.739   17.146   21.078   25.004   21.766   24.403   20.480   23.164   16.613    5.960
	      read(ifile,'(i3,12f9.3)')month,field(iy,:,id)
	   enddo
	enddo
end subroutine

!--calc climate ave for one station
subroutine calc_cli(F,F_cli)
	real, dimension(Nyear,12,31) :: F
	real, dimension(12) :: F_cli, total
	total=0
	F_cli=0
	do im=1,12
	   do iy=1,Nyear
	     do	id=1,31
		    if (F(iy,im,id).ne.Undef) then
			    F_cli(im)=F_cli(im)+F(iy,im,id)
				total(im)=total(im)+1
			endif
		 enddo
	   enddo
	   F_cli(im)=F_cli(im)/total(im)
	enddo

end subroutine

!--calc monthly climate total for one station
subroutine calc_cli_tot(F,F_cli)
	real, dimension(Nyear,12,31) :: F
	real, dimension(12) :: F_cli
	F_cli=0
	do im=1,12
	   do iy=1,Nyear
	     do	id=1,31
		    if (F(iy,im,id).ne.Undef) then
			    F_cli(im)=F_cli(im)+F(iy,im,id)
			endif
		 enddo
	   enddo
	   F_cli(im)=F_cli(im)/Nyear
	enddo
end subroutine

!--calc Monthly maximum consecutive x-day  (e.g. x=1,x=5)
subroutine calc_cli_xday(F,F_cli,x)
    integer :: x
	real, dimension(Nyear,12,31) :: F
	real, dimension(Nyear,12) :: Fmaxxday
	real, dimension(12) :: F_cli
	real :: Fxday
	F_cli=0
	Fmaxxday = 0
	!First calc max x-day values for each month,year
	do im=1,12
	   do iy=1,Nyear
		 Fmaxxday(iy,im) = 0
	     do	id=1,31-x
		    call usum(F(iy,im,id:id+x),x,Fxday)
			if (Fxday.gt.Fmaxxday(iy,im)) Fmaxxday(iy,im) = Fxday
		 enddo
	   enddo
	enddo

	!Then, calculate the average value
	do im=1,12
	   F_cli(im) = sum(Fmaxxday(:,im))/Nyear
	enddo
end subroutine

!--Ngreater
!--Calc. the climatological of number of day that is greater than or equal a threhold
subroutine calc_ge_cli(F,thre,N_cli,flag)
    real, dimension(12) :: thre	  ! threshold of each month
	real, dimension(Nyear,12,31) :: F
	integer, dimension(Nyear,12) :: N
	real, dimension(12) :: N_cli
	integer :: flag                 ! If flag =1, use less than instead of greater than

	!First, calculate the number of day that is greater than a threhold
	do im=1,12
	   do iy=1,Nyear
		 N(iy,im) = 0
	     do	id=1,31
		    if (F(iy,im,id).ne.Undef) then
			    if (flag.eq.0) then
			       if (F(iy,im,id).ge.thre(im)) N(iy,im) =  N(iy,im)+1
				else
				   if (F(iy,im,id).le.thre(im)) N(iy,im) =  N(iy,im)+1
				endif
			endif
		 enddo
	   enddo
	enddo
	
	do im=1,12
	   N_cli(im) = 1.*sum(N(:,im)) / Nyear
	enddo

end subroutine


!-Caclate monhtly percentile
subroutine calc_mon_percentile(F_total,N_total,k,F_per)
    real :: k
	real, dimension(12,2000) :: F_total	
	integer, dimension(12) :: N_total
	real, dimension(12) :: F_per

	real, dimension(2000) :: Atemp

	do im=1,12
		Atemp(:)=F_total(im,:)
		call percentile(Atemp,N_total(im),k,F_per(im))
	enddo
     
end subroutine


!--- Concat all data to get a long monthly array for all year
subroutine total_series(F,F_total,N_total)
	real, dimension(Nyear,12,31) :: F
	real, dimension(12,2000) :: F_total    ! maximum ~65 years
	integer, dimension(12) :: N_total
	F_total=0
	N_total=0
	do im=1,12
	   do iy=1,Nyear
	      do id=1,31
		     if (F(iy,im,id).ne.Undef) then
			    N_total(im)=N_total(im)+1
			    F_total(im,N_total(im))	= F(iy,im,id)
			 endif
		  enddo
	   enddo
	enddo

end subroutine

!Calculate percentile a SORTED data
!e.g. k is from 0-100, note
subroutine percentile(F,n,k,per)
    integer :: n, ip
    real :: k, per
	real, dimension(n) :: F
	call bubble(F,n)
	ip=NINT(k*n/100.)
	per = F(ip)
end subroutine


!--calculate sum of an array (don't cout Undef values)
subroutine usum(F,n,ave)
   integer :: n,i,total
   real, dimension(n) :: F
   real :: ave
   ave=0
   do i=1,n
      if (F(i).ne.Undef) then
		ave=ave+F(i)
	  endif
   enddo
end subroutine 
!return p,q in ascending order
Subroutine Order(p,q)
real :: p, q, temp

  if (p>q) then
    temp=p
    p=q
    q=temp
  end if
  return
end	 Subroutine

!Buuble sorting of integer array A
Subroutine Bubble(A, n)
integer :: n,i,j
real,dimension(n) :: A
  do i=1, n
    do j=n, i+1, -1
      call Order(A(j-1), A(j))
    end do
  end do
  return
end	  Subroutine




End