program main

!*****************************************************************************80
!! MAIN is the main program for KMEANS_TEST.
!  Discussion:
!    KMEANS_TEST tests the KMEANS library.
!  Licensing:
!    This code is distributed under the GNU LGPL license.
!  Modified:
!    06 October 2009
!  Author:
!    John Burkardt
!
  Use NetCDF
  implicit none
  Real, Parameter :: Rmiss=-99.0
  Integer, Parameter :: Imiss=-99
  Real, Allocatable :: X(:,:), Lon(:), Lat(:), Reg(:,:), xtime(:)
  Real, Allocatable :: Mask_Land(:), Mask(:,:), xlon(:),xlat(:)
  Integer :: Landpoint, Nst, ist, Nlat, Nlon, ilat, ilon, Ntime,i,j, irow, icol, N_Clus

  Real (kind=8), Allocatable :: point(:,:)   ! (M,N)
  Real, Allocatable :: C(:,:), XX(:,:)
  Integer,Allocatable :: Z(:)

  Integer (kind=4) :: M, N, iclu, cluster_num  ! Dims, Points, Num of clusters
  Integer (kind=4) :: Min_Obs, Max_Clu
  Logical :: status

  Integer, Allocatable :: Clus(:,:)      ! (N,N_Clus)
  Real,    Allocatable :: Clus_var(:,:)  ! (N_Clus,N_Clus)
  Real,    Allocatable :: Mean_var(:)   ! (N_Clus)
  Real,    Allocatable :: Total_Ene(:)   ! (N_Clus)
  Integer, Allocatable :: Clus_popu(:,:) ! (N_Clus,N_Clus)
  Integer(kind=4), Allocatable :: Popu(:) ! (N_Clus,N_Clus)
  Integer(kind=4), Allocatable :: Clus_out(:) ! (N)      - B
  Real (kind=8), Allocatable :: Vari(:)   ! (N_Clus)     - DEV
  Real (kind=8), Allocatable :: Center(:,:) ! (N_Clus,M) - D

  !Call Proc_land_mask_AWC  ! Return Mask, LandMask, Lon(Nst),Lat(Nst)
  Call Read_data           ! Return: X(Landpoint,Ntime), xlon(Nlon), xlat(Nlat) 
  print*,"Read data Done..."
  print*,"Num of points, Num of Var: ", Landpoint,Ntime

  M=Ntime; N=Landpoint
  cluster_num=5
  Min_Obs=5; Max_Clu=cluster_num

  Allocate (point(N,M),Popu(cluster_num),Clus_out(N),Vari(cluster_num),Center(cluster_num,M)) ! , XX(M,N))
  !Do irow=1,N
  !   point(1:M,irow) = X(irow,1:M)
  !   XX(1:M,irow) = X(irow,1:M)
  !Enddo

  !cluster_num=10
  !Call clustr ( point, d, Vari, Clus_out, Popu, N, M, cluster_num, 1, N )

  point(:,:)=X(:,:)
  N_Clus=cluster_num
  Allocate (Clus(Landpoint,N_Clus),Clus_var(N_Clus,N_Clus),Clus_popu(N_Clus,N_Clus))
  Allocate (Total_Ene(N_Clus),Mean_var(N_Clus))
  Clus=Imiss; Clus_popu=Imiss; Clus_var=Rmiss; Mean_var=Rmiss

  Do j=1,N_Clus
     iclu=j
     print*,"Case ",j

     Call clustr ( point, Center, Vari, Clus_out, Popu, N, M, iclu, Min_Obs, Max_Clu )
     print*,"Clusters: ",Clus_out
     Clus(1:Landpoint,j)=Clus_out(1:N)  ! Clusters

     Clus_popu(1:j,j)=Popu(1:iclu)  ! Population of each cluster
     Mean_var(j)=sum(Vari(1:iclu))/iclu
     Clus_var(1:j,j)=Vari(1:iclu)

     !Clus_popu(1:N_Clus,j)=Popu(1:cluster_num)  ! Population of each cluster
     !Mean_var(j)=sum(Vari(1:cluster_num))/cluster_num
     !Clus_var(1:N_Clus,j)=Vari(1:cluster_num)
  Enddo

 ! For writing to file:
  Allocate (Reg(Nlon,Nlat))
  Reg(:,:)=Rmiss   ! Regions
  ist=0
  Do ilat=1,Nlat
     Do ilon=1,Nlon
        if (Mask(ilon,ilat)/=Rmiss) then
           ist=ist+1
           Reg(ilon,ilat)=Clus(ist,N_Clus)
        endif
     Enddo
  Enddo
  Call Write_Clus_to_NC

 ! Clusters:
  print*,"Clusters"
  print"(3A7,20(3x,A1,I2.2))","STT","Lon","Lat",(("C",i),i=1,N_Clus)
  Do i=1,Landpoint
     print"(I7,2F7.3,20I6)",i,Lon(i),Lat(i),(Clus(i,1:N_Clus))
  Enddo

 !Populations:
  print*,"Populations"
  print"(A5,20(6x,A1,I2.2))","Gr",(("C",i),i=1,N_Clus)
  Do i=1,N_Clus
     print"(I5,20I9)", i, (Clus_popu(i,1:N_Clus))
  Enddo

 ! Variances:
  print*,"Variances"
  print"(A5,20(9x,A1,I2.2))","Gr",(("C",i),i=1,N_Clus)
  Do i=1,N_Clus
     print"(I5,20F12.1)", i, (Clus_var(i,1:N_Clus))
  Enddo

 ! Total variances:
  print*,"Mean variances"
  print"(A3,2A15)","Gr","Mean_Var" !,"Tot_Ene"
  Do i=1,N_Clus
     print"(I2,20F15.2)", i, Mean_var(i) !, Total_Ene(i)
  Enddo

  pause

  Deallocate (X,Lon,Lat,xtime)
  Deallocate (point, Center,Vari, Clus_out, Popu)
  Deallocate (Clus,Clus_var,Clus_popu,Total_Ene,Mean_var)

  stop 0
  !-----------------------------------------------------------------------------
  CONTAINS
      !--------------------------------------------------------------------
      Subroutine Write_Clus_to_NC
        Character (Len=200) :: Fname="Clus_pca_pdsi.nc"

        Call wr_ncdf_2d (Fname,xlon,xlat,nlon,nlat,Reg,"reg","Drought regions","-","lon","lat", Rmiss)
      End Subroutine Write_Clus_to_NC
      !--------------------------------------------------------------------
      Subroutine Read_data
        Character (Len=200) :: Fname
        Real, Allocatable :: PDSI(:,:,:), tmp(:)
        Integer :: ist, ilon, ilat, ipoint

        !Fname="/home/tanpv/Progs/Dr_SEA/Results/BL_APH.nc"
        !Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,ntime,PDSI,"ori_pdsi")

        Fname="/home/tanpv/Progs/Dr_SEA/EOF_New/pca_APH.nc"
        Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,ntime,PDSI,"pdsi")

        Nst=Nlon*Nlat
        Allocate (Mask(Nlon,Nlat))
        Allocate (lon(Nst), lat(Nst), tmp(Ntime))
        ist=0; Landpoint=0; Mask(:,:)=Rmiss
        Do ilat=1, Nlat
           Do ilon=1, Nlon
              ist=ist+1;
              Lon(ist)=xlon(ilon); Lat(ist)=xlat(ilat)
              tmp(1:Ntime)=PDSI(ilon,ilat,1:Ntime)
              if (.not.Ar1_Empty(tmp,Ntime)) then
                 Landpoint=Landpoint+1; Mask(ilon,ilat)=1
              endif
           Enddo
        Enddo
        Allocate (X(Landpoint,Ntime))
        ipoint=0
        Do ilat=1, Nlat
           Do ilon=1, Nlon
              !tmp(1:Ntime)=PDSI(ilon,ilat,1:Ntime)
              if (Mask(ilon,ilat)/=Rmiss) then
              !if (.not.Ar1_Empty(tmp,Ntime)) then
                 ipoint=ipoint+1
                 X(ipoint,1:Ntime)=PDSI(ilon,ilat,1:Ntime)
              endif
           Enddo
        Enddo
      End Subroutine Read_data
   !------------------------------------------------------------------------------
   Subroutine Proc_land_mask_AWC
   Integer :: ist, ilon, ilat
   Real, Allocatable :: xtmp(:,:),xlon(:),xlat(:)
   Character (Len=200) :: MaskFname="/work/users/hmo/tanpv/CORDEX_SEACLID/Land_mask_sea.nc"

     Call rd_ncdf_2d (trim(MaskFname),xlon,xlat,Nlon,Nlat,Mask, "land")

     Nst=Nlon*Nlat
     Allocate (Mask_Land(Nst))
     Allocate (lon(Nst), lat(Nst))
     ist=0; Landpoint=0
     Do ilat=1, Nlat
        Do ilon=1, Nlon
           ist=ist+1; Mask_Land(ist)=Mask(ilon,ilat);
           Lon(ist)=xlon(ilon); Lat(ist)=xlat(ilat)
           if (Mask_Land(ist)/=Rmiss) Landpoint=Landpoint+1
        Enddo
     Enddo
     Deallocate (xlon,xlat)
   End Subroutine Proc_land_mask_AWC
   !------------------------------------------------------------------------------
   Logical Function Ar1_Empty(X,N)
     Integer :: N, i
     Real :: X(N)
     Logical :: Empty
     Empty=.false.
     Do i=1,N
        if (X(i)==Rmiss) then
           Empty=.true.; exit
        endif
     Enddo
     Ar1_Empty=Empty
   End Function Ar1_Empty
   !------------------------------------------------------------------------------
  Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
  Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
!!!end
!-----------------------------------------------------------------------------
 !Include "kmean_lib.f90"
 Include "kmean_asa058.f90"
!-----------------------------------------------------------------------------
!------------------ End Prog. ------------------
END