!====================================================================================================!
Program Rape_Data
Implicit None
Integer :: ii,i,Nmon,i5,NDay, iyr,Yr1, Yr2
Integer :: Num_Var, Year, Nst
Character (Len=3), Allocatable  :: Var(:)
Character (Len=21), Allocatable :: StaName(:)
Character (Len=80) :: filename,fntem,fntem2
Integer, Allocatable :: Reg(:)                ! Climatological Sub-regions:
Real, Allocatable :: Lon(:), Lat(:), H(:) 
Real, Allocatable :: Dat(:,:,:)      !  Array stores data of one year for one var of one station (Nday,Nmon,Yr1:Yr2)
real, Allocatable :: Dat4d(:,:,:,:,:)  ! Array stores data of one year for one var of all station (isc,Nday,Nmon,Yr1:Yr2)
real, Allocatable :: mdhtmp(:,:,:,:)
Character (Len=1), Allocatable :: Data_Map(:,:,:) !  Data info
Integer :: ivar, ist, iday, imon, Flag
!====================================================================================================!
integer, parameter :: nsta=64, styr=1981, enyr=2010
real, dimension (nsta) :: oslat, oslon, osheight
character(len=12) :: osname(nsta)
real :: q33,q66,TB
real, dimension(800,20) ::  percent33,percent66
real, dimension(800,20,12) ::  mdhpercent33,mdhpercent66
integer :: np,isc,ik,j,l
real, dimension (500000) :: pdata
real, dimension (31) :: mdh4sub
Character (len=21) ::  sname(800)
Character(len=50) :: dirin, dirout
data dirin /"./"/
data dirout /"./output/"/
!====================================================================================================!
! Ghi chu cua MDH:
! Input: tai ./ gom:
!   + Obs_Daily_Bin_70.dat:        so lieu dang binary cua cac yeu to tai cac tram
!   + Station_info_observed64.txt: thong tin cac tram
! Output: tai ./output/ gom:
!   + Station_info.txt    : thong tin ve cac tram
!   + TENTRAM_TENBIEN.txt : (rat nhieu) tram voi cac bien tuong ung 
!   + TRUNGBINH_Tx.txt
!   + TRUNGBINH_Tm.txt
!====================================================================================================!
Open(1,file=trim(dirin)//"Obs_Daily_Bin_70.dat",form="unformatted",status="Old")
Read(1) Num_Var       ! Integer (Num of Var)
Read(1) Nst           ! Integer (Num of Station)
Read(1) Yr1, Yr2      ! Integer (Start and End year of the period)
Read(1) Nmon, NDay    ! Integer (Num of Mon in the year and num of day in the mon)
!-------------------------------------------------------------------------------------
Allocate (Dat(Yr1:Yr2,Nmon,Nday))
Allocate (Dat4d(nst,Num_Var,Yr1:Yr2,Nmon,Nday))
Allocate (mdhtmp(nst,Num_Var,Yr1:Yr2,Nmon))
Allocate (Data_Map(Num_Var,Nst,Yr1:Yr2))
Allocate (Reg(Nst), Var(Num_Var), StaName(Nst))
Allocate (Lon(Nst), Lat(Nst), H(Nst))
!-------------------------------------------------------------------------------------
Data_Map(:,:,:) = "-"
StaName(:) = ""; Reg(:) = -99; Lon(:) = -99; Lat(:) = -99
print*,"Num_var, Nst=", Num_Var, Nst
print*, Yr1, Yr2, Nmon, NDay
open(25,file=trim(dirout)//'Station_info.txt',status='unknown')
write(25,*) 'Name                   Lon      Lat    Eliv Region'
Do ivar=1, Num_Var      ! (=8)
   read(1) Var(ivar)
   Do ist=1, Nst
      read(1) Flag
      if (Flag == 1) then  ! data Var(ivar) are avialabe at this station
         read(1) StaName(ist)
         read(1) Lon(ist), Lat(ist), H(ist)
         read(1) Reg(ist)
         if(ivar.eq.2) write(25,'(A21,2F8.3,F8.1,I4)') StaName(ist),Lon(ist), Lat(ist), H(ist),Reg(ist)
         read(1) Dat
         Dat4d(ist,ivar,:,:,:)=Dat(:,:,:)
         Call Map (ivar, ist)
      endif
   Enddo
Enddo
Close (25)
!====================================================================================================!
! Write data to files
!====================================================================================================!
isc=0
Do ist=1, Nst
  fntem=StaName(ist)
  call getc(fntem,ii)
  if(ii.gt.1) isc=isc+1
  Do ivar=1, Num_Var
    fntem=StaName(ist)
    call getc(fntem,ii)
    If(ii.gt.1) then
      fntem2=fntem(1:ii) // '_' // Var(ivar)
      call getc(fntem2,i5)
      filename=trim(dirout)// fntem2(1:i5) // '.txt'
      call getc(filename,ii)
      write(*,*) ist, StaName(ist),filename(1:ii)
      open(90,file=filename(1:ii),status='unknown')
      write (90,*) 'R: Daily rainfall amount'
      write (90,*) 'Tx: Daily Max Temp'
      write (90,*) 'Tm: Daily Min Temp'
      write (90,*) 'T2m: Daily 2m Temp'
      write (90,*) 'BH: Daily Potential Evaporation'
      write (90,*) 'Um: Daily Min Relative Humidity'
      write (90,*) 'U13: Relative Humidity at 13h'
      write (90,*) 'Vx: Daily Max Wind Speed'
      write (90,*) 'Variable  Station_Name Lat  Lon Hight Rgion_ID'
      write (90,'(a3,x,A21,2F8.3,F8.1,I4)') Var(ivar),StaName(ist),Lon(ist),Lat(ist),H(ist),Reg(ist)
      q33=0.33
      q66=0.66
!      np=0
!      pdata=0.0
      Do iyr=Yr1,Yr2
        write (90,'(i4)') iyr
! Them doan tinh trung binh thang de thay cach tinh percentile 
        Do imon=1,Nmon
          l=0
          mdh4sub=flag
          Do iday=1,Nday
            if(Dat4d(ist,ivar,iyr,imon,iday).gt.-9.0.and.iyr.ge.styr.and.iyr.le.enyr) then
              l=l+1
              mdh4sub(l)=Dat4d(ist,ivar,iyr,imon,iday) 
            endif 
          Enddo  
          mdhtmp(ist,ivar,iyr,imon)=TB(mdh4sub,l,flag)
        Enddo
! Doan nay giu nguyen de ghi file
        Do iday=1, Nday
          Do imon=1,Nmon
            if(Dat4d(ist,ivar,iyr,imon,iday).gt.-9.0.and.iyr.ge.styr.and.iyr.le.enyr) then
              np=np+1
              pdata(np)=Dat4d(ist,ivar,iyr,imon,iday)
            endif
          Enddo
          write (90,'(I3,12F7.1)') iday, (Dat4d(ist,ivar,iyr,imon,iday), imon=1,Nmon)
        Enddo
      sname(isc)=StaName(ist)
      Enddo
! Sua cach tinh percentile
      Do imon=1,Nmon 
        np=0
        pdata=0.0
        Do iyr=styr,enyr
          if (mdhtmp(ist,ivar,iyr,imon).gt.-9.0) then
               np=np+1
               pdata(np)=mdhtmp(ist,ivar,iyr,imon)
          endif
        Enddo
        if(Var(ivar)(1:1).eq.'R')then
          mdhpercent33(isc,ivar,imon)=percentile_1(pdata,np,q33,.true.)
          mdhpercent66(isc,ivar,imon)=percentile_1(pdata,np,q66,.true.)
        else
          mdhpercent33(isc,ivar,imon)=percentile_1(pdata,np,q33,.false.)
          mdhpercent66(isc,ivar,imon)=percentile_1(pdata,np,q66,.false.)
        endif
      Enddo
! ========================= !
      close(90)
    Endif
  Enddo
Enddo
!====================================================================================================!
open(22,file=trim(dirout)//"observed_q33q66.txt",status='unknown')
write(22,*) 'PLEASE NOTE: TO COMPARE WITH MODEL, the observed q33,q66 is compute from 1981-2010 only'
do ivar=1,Num_Var
  write(22,*) 'q33 q66 for ', Var(ivar)
  do imon=1,Nmon
      write(22,*) 'Mon :', imon
      do i=1,isc
        write(22,'(a15,2F10.2)') sname(i)(1:15),mdhpercent33(i,ivar,imon),mdhpercent66(i,ivar,imon)
      enddo
  enddo
enddo
close (22)
!====================================================================================================!
! Get only used 64 station percentiles 
!=====================================================================!
! Read station informations
!=====================================================================!
open(3,file='Station_info_observed64.txt',status='old')
read(3,*)
do i=1,nsta
  read(3,*,end=55) osname(i), oslat(i),oslon(i), osheight(i)
!  write(*,*) sname(i), slat(i),slon(i), sheight(i)
enddo
close (3)
55 continue
!=====================================================================!
open(22,file=trim(dirout)//"Q33Q66_T2m.txt",status='unknown')
write(22,*) 'PLEASE NOTE: TO COMPARE WITH MODEL, the observed q33,q66 is compute from 1981-2010 only'
do ivar=1,Num_Var
  if (trim(Var(ivar)).eq.'T2m') then
    write(22,*) 'q33 q66 for ', Var(ivar)
    do imon=1,Nmon
    write(22,*) 'Mon:  ', imon
    do j=1,nsta
      do i=1,isc
        if (trim(sname(i)).eq.trim(osname(j))) then
          write(22,'(a15,2F10.2)') sname(i)(1:15),mdhpercent33(i,ivar,imon),mdhpercent66(i,ivar,imon)
        endif
      enddo
    enddo
    enddo
  endif
enddo
close (22)
!====================================================================================================!
!Call Write_Data_Map
Deallocate(Dat)
Deallocate(Dat4d)
CONTAINS
!---------------------------------------------------------------------
Subroutine Write_Data_Map
Integer :: iyr, iyt, ist
Character (Len=30) :: Fmt
Character (Len=40) :: Var_desc
Open (4,file="Data_Map.txt")
Do iyt=1, Num_Var
   Call Var_Description (Var(iyt), Var_desc)
   write (4,"(/A100)") "Available Data of "//trim(Var_desc)//" ("//trim(Var(iyt))//")"
   Write (Fmt,"(A,I2.2,A)") "(A15,A3,",Yr2-Yr1+1,"(1x,I2.2))"
   write (4,Fmt) "Station Name   ","Reg",(mod(iyr,100), iyr=Yr1, Yr2)
   Do ist=1, Nst
      if (Empty_Sta(iyt, ist)) Cycle
      Write (Fmt,"(A,I2.2,A)") "(A15,I2,1x",Yr2-Yr1+1,"(1x,A1,1x))"
      write (4,Fmt) StaName(ist), Reg(ist), (Data_Map(iyt,ist,iyr), iyr=Yr1, Yr2)
   Enddo
Enddo
Close (4)
End Subroutine Write_Data_Map
!---------------------------------------------------------------------
Logical Function Empty_Sta(iyt, ist)
Integer :: iyr, iyt, ist
Logical :: Empty
Empty = .true.
Do iyr=Yr1, Yr2
   if (Data_Map(iyt, ist, iyr) == "x") then
      Empty = .false.
      exit
   endif
Enddo
Empty_Sta = Empty
End Function Empty_Sta
!---------------------------------------------------------------------
Subroutine write_data_one_station
  Integer :: iyr, imon, iday
  Do iyr=Yr1,Yr2
     write (26,*) iyr
     Do iday=1, Nday
      write (26,'(I3,12F7.1)') iday, (Dat(iyr,imon,iday), imon=1,Nmon)
     Enddo
  Enddo
End Subroutine write_data_one_station
!---------------------------------------------------------------------
Subroutine Map (iyt, ist)
Integer :: iday, imon, iyr, iyt, ist
Logical :: Empty
Do iyr=Yr1, Yr2
 Empty = .true.
 Do imon=1,Nmon
   Do iday=1,Nday
      if (Dat(iyr,imon,iday) /= -99.) then
         Empty = .false.
         exit
      endif
   Enddo
 Enddo
 if (.not.Empty) Data_Map(iyt, ist, iyr) = "x"
Enddo 
End Subroutine Map
!---------------------------------------------------------------------
Subroutine Var_Description(Varname, Descr)
Character (Len=*) :: Varname, Descr
Select Case (trim(Varname))
  Case ("R")
    Descr = "Daily rainfall amount"
  Case ("Tx")
    Descr = "Daily Maximum Temperature"
  Case ("Tm")
    Descr = "Daily Minimum Temperature"
  Case ("T2m")
    Descr = "Daily 2m mean Temperature"
  Case ("BH")
    Descr = "Daily Potential Evaporation"
  Case ("Um")
    Descr = "Daily Minimum Relative Humidity"
  Case ("U13")
    Descr = "Relative Humidity at 13h"
  Case ("Vx")
    Descr = "Daily Maximum Wind Speed"
  Case ("Vtb")
    Descr = "Daily mean Wind Speed"
End Select
End Subroutine Var_Description
!---------------------------------------------------------------------
        subroutine getc(charvar,icount)
        character*80 :: charvar
        integer :: icount
        icount=0
        do i=1,80
        if(charvar(i:i).ne.' ') then
         icount=icount+1
        endif
        enddo
        return
        end subroutine getc
!---------------------------------------------------------------------
real function percentile_1(X1,N1,p,Rain)
  Integer :: N, i, j, N1
  Real :: X1(N1), X(n1), tmp, p, per  ! (p value mustbe in [0-1])
  Logical :: Rain
  if (Rain) then
     n=0
     Do i=1,N1; if (X1(i)>=1.) then; n=n+1; X(n)=X1(i); endif; Enddo
  else
     n=N1; X(1:N) = X1(1:N1)
  endif
  Do i = 1, n-1
     Do j = i+1, n
        if (x(j) < x(i)) then
           tmp = x(j)
           x(j) = x(i)
           x(i) = tmp
        endif
     Enddo
  Enddo
  tmp = p*(n+1)
  if (tmp < 1) then
     per = x(1)
  elseif (tmp > n) then
     per = x(n)
  else
     i = int(tmp)
     j = i+1
     per = x(i) + (x(j)-x(i))*(tmp-real(i))
  endif
  percentile_1 = per
End Function Percentile_1
!====================================================================================================!
END
	!=================================================================================!
	! Subroutine va Function di kem
	!=================================================================================!
	Function TB(X,N,flg)
	Integer :: N,i,flg
	Integer :: count
	Real*4 X(N), temp
	temp = 0.
	count = 0
	Do i = 1,N
	    If (X(i) /= flg) then
		count = count + 1
		temp = temp + X(i)
	    Endif
	Enddo
	TB = temp/real(count)
	Return
	End Function
