! Chuong trinh tinh PCs, EOFs theo cach moi (2019-10-10)
      Use NetCDF
      Implicit None
      Real, Parameter :: Miss=-99.0, Lat_min=-11.0  ! Cat phan Autralia
      Integer, Parameter :: Nidx=24

      Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:)
      Real*8, Allocatable :: EOFs(:,:), PCs(:,:)
      Real*8, Allocatable   :: Eigenvalues(:), ECumulate(:)
      !---------------------------------------------------------------------------------

      Integer :: Yr1, Yr2, Nmon, Nyear ! Thoi ky
      Integer :: Nt, Nx   ! Kich thuoc thuc cua mang X
      Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K
      Integer :: switch1, standard1, PC_Stard1, Nmode1    ! Clim
      Integer :: switch2, standard2, PC_Stard2, Nmode2    ! sst
      Character (Len=250) :: Fname, Fname1, MaskFname, pca_out, eof_out
      Character (Len=250) :: Fnclim, pca_clim, eof_clim
      Character (Len=250) :: FnSST, pca_sst, eof_sst, pca_prec, eof_prec
      Real :: Slon1, Slon2, Slat1, Slat2, Miss_sst           ! sst domain
      Integer :: Yr_sta

      Integer :: i,j,L
      ! Mo ta thong tin SL dau vao:
      !---------------------------------------------------------------------------------
      Real, allocatable :: xlon(:), xlat(:), xtime(:)
      Integer :: Nlon, Nlat, Ntime, ilon, ilat
      Logical :: Status

      Real, allocatable :: PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:)
      Real, allocatable :: ClimIDX(:,:)
      Real, allocatable ::  Mask(:,:), Mask_all(:,:), Mask_dr(:,:)
      Real, allocatable :: Dr_IDX(:,:,:,:)
      !---------------------------------------------------------------------------------
      ! SST:
      Real, allocatable :: lon_sst(:), lat_sst(:)
      Real, allocatable :: sst(:,:), Mask_sst(:,:)   ! (Nt,Nx)
      Integer :: Nlon_sst, Nlat_sst, SeaPoint, Mov_Ave_Step
      !---------------------------------------------------------------------------------
      ! Composite maps of SST:
      Real, allocatable :: sst_dry(:,:), sst_wet(:,:), sst_nor(:,:)
      Real, allocatable :: sst_cli(:,:)   ! (Nlon_sst, Nlat_sst): SST with Areal mean of SPEI<-0.5 and >-0.5
      Real, allocatable :: sst_3d(:,:,:)                ! (Nlon_sst, Nlat_sst,Ntime)
      Real, allocatable :: Areal_mean_Dr(:),sst_1d(:)                    ! (Ntime)
      Real, allocatable :: sst_N(:), sst_S(:), sst_G(:)
      Real, allocatable :: Cor_Dr_SSTG(:,:), Pval_Dr_SSTG(:,:)
      !---------------------------------------------------------------------------------
      ! Rainfall:
      Real, allocatable :: Prec(:,:), lon_prec(:), lat_prec(:), Mask_prec(:,:)
      Integer :: Ntime_prec, Nlon_prec, Nlat_prec, PrecPoint

      Integer :: Nst, Nst1, iidx
      Integer, Parameter :: NDr_idx=1 ! PDSI, SPEI12, SPI12
      Character (Len=5), Parameter :: IDX_nm(NDr_idx)= &
                                      (/"spe12"/) !,"spi12"/)
      Character (Len=255) :: Filename
      !---------------------------------------------------------------------------------
      ! Reading data:
      !Call rd_ncdf_2d (trim(MaskFname),xlon,xlat,Nlon,Nlat,Mask, "land")
      !Do ilat=1,Nlat
      !   if (xlat(ilat) < Lat_min) Mask(:,ilat) = Miss
      !Enddo
      !print*,trim(Fname)
      !Call Read_Drough_IDX (Fname)

      Yr1=1960; Yr2=2019; Nmon=12
      Nyear = Yr2-Yr1+1            ! So nam co SL
      Ntime = Nyear*Nmon 
      Allocate (sst_N(Ntime), sst_S(Ntime), sst_G(Ntime), sst_1d(Ntime))
      !---------------------------------------------------------------------------------
      ! Composite SST maps:
      !---------------------------------------------------------------------------------
      ! SST info:
      Yr_sta     = 1891
      Miss_sst   = 1.e+20         ! Missing values
      Slon1      = 190
      Slon2      = 210
      Slat1      =  30
      Slat2      =  40
      !Slat2      =  50

      Fname="/work/projects/G09/SEACLID/SST/sst.mon.mean_COBE.nc"
      Call Read_sst_data (Fname)
      sst_N(1:Ntime) = sst_1d(1:Ntime)
      Slon1      = 190
      Slon2      = 210
      Slat1      = -5
      Slat2      =  5
      Call Read_sst_data (Fname)
      sst_S(1:Ntime) = sst_1d(1:Ntime)
      sst_G(1:Ntime) = sst_S(1:Ntime) - sst_N(1:Ntime)

      Fname= &
     "/work/projects/G09/SEACLID/GPCC_CRU_SEA/SPEI_SPI_GPCC_CRU_59_19.nc"
      Call Read_Drough_IDX (Fname)

      Fname="/home/tanpv/Progs/2021/NC/Cor_Dr_SSTG.nc"
      print*,"Calculating Corr Dr and SSTG..."
      Call Corre_Dr_SSTG (Fname)
      CONTAINS
      !================================================================================= 
      Subroutine Anomaly (X,Ntime)
      Integer :: Ntime
      Real :: X(Ntime), Y(Ntime), TB(Nmon)
      Integer :: ic, it, imon, imon1
      Real :: tmp

      if (mod(Ntime,Nmon)/=0) print*,"Something wrong !"
      Do imon=1,Nmon
         tmp=0; ic=0
         Do it=1,Ntime
            imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon
            if (imon1==imon) then
               if (X(it)/=Miss) then
                  ic=ic+1; tmp=tmp+X(it)
               endif
            endif
         Enddo
         if (ic>0) then; TB(imon)=tmp/ic; else; TB(imon)=Miss; endif
      Enddo
      Do it=1,Ntime
         imon=mod(it,Nmon); if (imon==0) imon=Nmon
         Y(it)=X(it)-TB(imon)
      Enddo
      X(:) = Y(:)
      End Subroutine Anomaly
      !--------------------------------------------------------------------
      Subroutine Correlation(X,Y,N,Corr,prob,Miss)
        Integer :: N, N1, i, j
        Real :: X(n), Y(n), Miss, Xtb, Ytb, Sx, Sy, Rxy, tmp, t
        Real :: Prob, Corr , T_Dist

        N1=0; Xtb=0; Ytb=0; Sx=0; Sy=0; Rxy=0
        Do i=1,N
           if (X(i)/=Miss.AND.Y(i)/=Miss) then
              n1=n1+1
              Xtb=Xtb+X(i); Sx=Sx+X(i)*X(i)
              Ytb=Ytb+Y(i); Sy=Sy+Y(i)*Y(i)
              Rxy=Rxy+X(i)*Y(i)
           endif
        Enddo
        if (N1>2) then
           Xtb=Xtb/Real(N1); Ytb=Ytb/Real(N1)
           Sx=Sx/Real(N1)-Xtb*Xtb; Sy=Sy/Real(N1)-Ytb*Ytb
           tmp=Rxy/Real(N1)
           Rxy=(tmp-Xtb*Ytb)/Sqrt(Sx*Sy)
           t = ABS(Rxy/ ( Sqrt(1-Rxy**2) / Sqrt(Real(N1-2)) ) )
           Prob = T_Dist (t, Real(N1-2))
        else
           Rxy=Miss; prob=Miss
        endif
        Corr=Rxy

      End Subroutine Correlation
      !---------------------------------------------------------------------------------
      Subroutine Read_Drough_IDX (Fname)
         Integer :: ist, ilon, ilat, it, it1, it2, Nt, imon, imon1, iyr
         Integer :: nlat1, nlat0! , nlon, nlat
         Character (Len=*) :: Fname
         Real :: tmp, xmean(12), ic(12)
         Real, Allocatable :: X0(:), X1(:,:,:), X2(:,:,:), X3(:,:,:)
         !Real, Allocatable :: xlon(:), xlat(:)

         Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPEI12, "spei12")
         !Deallocate (xlon, xlat, xtime)
         Where (IsNaN(SPEI12)); SPEI12=Miss; Endwhere

         Ntime = Nt-Nmon     ! Skip first year due to SPI and SPEI were missed
         it1 = Nmon+1; it2 = Nt

         Allocate (X1(nlon, nlat, Ntime),X2(nlon, nlat, Ntime),X3(nlon, nlat, Ntime))

         X2(:,:,1:Ntime)=SPEI12(:,:,it1:it2) !Nmon+1:Nt)

         Deallocate (SPEI12)
         Allocate (SPEI12(nlon,nlat,Ntime))

         SPEI12=X2;

         Allocate (Dr_IDX(nlon,nlat,Ntime,NDr_idx))

         Dr_IDX(:,:,:,:) = Miss

         Deallocate (X1, X2, X3)
         Allocate (Mask(nlon,nlat))
         Mask = Miss
         Do ilat=1, Nlat
            Do ilon=1, Nlon
               !if (Mask(ilon,ilat)>0) then
                  if (.not.Empty_Ar1(SPEI12(ilon,ilat,1:Ntime), Ntime)) then
                     Mask(ilon,ilat)=1
                     Dr_IDX(ilon,ilat,:,1)=SPEI12(ilon,ilat,:)
                  else
                     Mask(ilon,ilat)=Miss
                     SPEI12(ilon,ilat,:)=Miss
                  endif
               !endif
            Enddo
         Enddo
         Nst=ist   ! Number of gridpoints
         print*,"Reading data ... Done !"
      End Subroutine Read_Drough_IDX
      !--------------------------------------------------------------------------------------
      Subroutine Corre_Dr_SSTG (Fname)
        ! Lon_sst(Nlon_sst), Lat_sst(Nlat_sst)
        Character (Len=*) :: Fname
        Integer :: ilon,ilat
        Real, Allocatable :: X(:), Y(:)
        Real :: core, prob

        Allocate (Cor_Dr_SSTG(Nlon,Nlat), Pval_Dr_SSTG(Nlon,Nlat))
        Allocate (X(Ntime), Y(Ntime))

        Cor_Dr_SSTG=Miss; Pval_Dr_SSTG=Miss
        Y(1:Ntime) = sst_G(1:Ntime)
        print*,"Checking Nlon,Nlat,Ntime: ",Nlon,Nlat,Ntime; ! pause
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              if (Mask(ilon,ilat)/=Miss) then
                 X(1:Ntime)=SPEI12(ilon,ilat,1:Ntime)
                 Call Correlation(X,Y,Ntime,Core,prob,Miss)
                 Cor_Dr_SSTG (ilon,ilat) = Core
                 Pval_Dr_SSTG(ilon,ilat) = prob
              endif
           Enddo
        Enddo
        Call wr_ncdf_2d (Fname,xlon,xlat, Nlon,Nlat, Cor_Dr_SSTG, &
                     "cor","Correlation between SPEI and SSTG", &
                     "1", "lon","lat", Miss)
        Call wr_ncdf_2d_apd (Fname,Pval_Dr_SSTG, "pval", Nlon,Nlat, &
                 "p-values", "1", "lon", "lat", Miss)
      End Subroutine Corre_Dr_SSTG 
      !---------------------------------------------------------------------------------
      Subroutine Read_sst_data (Fname)
        Character (Len=*) :: Fname
        !Integer, Parameter :: Yr_sta=1891   ! Nam bat dau chuoi SL SST
        Real, Allocatable :: PDSI(:,:,:), tmp(:), XX(:,:,:),YY(:,:,:)
        Real, Allocatable :: X1d(:), Y1d(:), tlon(:), tlat(:), xtime(:)
        Real, Allocatable :: X2d(:,:)
        Integer :: ist, i, j, ilon, ilat, ipoint, nt, it, nlon, nlat
        Integer :: it1, it2, ilon1, ilon2,ilat1,ilat2, imon, M, Md2
        !Character (Len=255) :: Fn

        print*, "Reading SST ..." 

        Call rd_ncdf_3d (Fname,tlon,tlat,xtime,nlon,nlat,nt,XX,"sst") !trim(varn))
        Where (XX==Miss_sst) XX=Miss
        Where (IsNaN(XX)); XX=Miss; Endwhere
        Where (XX==Miss_sst); XX=Miss; Endwhere
        i=1; Do while (tlon(i)<Slon1); i=i+1; Enddo; ilon1=i
        i=1; Do while (tlon(i)<Slon2); i=i+1; Enddo; ilon2=i-1
        j=1; Do while (tlat(j)>Slat1); j=j+1; Enddo; ilat1=j
        j=1; Do while (tlat(j)>Slat2); j=j+1; Enddo; ilat2=j-1
        it1=(Yr1-Yr_sta)*Nmon+1 ! + Nmon         ! Bo nam dau tien
        it2=it1+Ntime-1
        print*,"SST: Yr_sta,Yr1,it1,it2,Ntime: ", Yr_sta,Yr1,it1,it2,Ntime
        ! Lay them M/2 phan tu dau va M/2 phan tu cuoi cho Moving Average:
        !M = Mov_Ave_Step
        !Md2=M/2
        !it1=it1-Md2; it2=it2+Md2
        Nt = Ntime                ! Luu gia tri Ntime
        Ntime = it2-it1+1         ! Kich thuoc chuoi moi dung trong Moving Average

        Nlon_sst=ilon2-ilon1+1;
        Nlat_sst=ilat1-ilat2+1       ! Chu y: ilat1>ilat2 vi Latitude giam tu North to South

        Allocate (YY(Nlon_sst,Nlat_sst,Ntime),Lon_sst(Nlon_sst), Lat_sst(Nlat_sst))
        Allocate (sst_3d(Nlon_sst,Nlat_sst,Nt))         ! Luu lai de su dung trong composite maps, Nt = do dai chuoi chua them Md2
        Allocate (Mask_sst(Nlon_sst,Nlat_sst), X1d(Ntime))
        Mask_sst(:,:)=Miss
        !X1d(1:Ntime) = xtime(it1:it2)
        i=0
        Do ilon=1,Nlon
           if (ilon>=ilon1.AND.ilon<=ilon2) then
              i=i+1; Lon_sst(i)=tlon(ilon)
              j=Nlat_sst+1
              Do ilat=1,Nlat
                 if (ilat>=ilat2.AND.ilat<=ilat1) then
                    j=j-1; Lat_sst(j)=tlat(ilat)
                    YY(i,j,1:Ntime) = XX(ilon,ilat,it1:it2)
                    ! Anomaly:
                    X1d(1:Ntime)=YY(i,j,1:Ntime)
                    Call Anomaly (X1d,Ntime)
                    YY(i,j,1:Ntime)=X1d(1:Ntime)
                 endif
              Enddo
           endif
        Enddo
        sst_3d(:,:,1:Nt) = YY(:,:,1:Nt) ! YY(:,:,Md2+1:Ntime-Md2)    ! Luu lai
        print*,"Checking Nt and Ntime: ", Nt,Ntime; if (Nt/=Ntime) stop
        Allocate (X2d(Nlon_sst,Nlat_sst))
        Do it=1,Ntime 
           X2d(:,:)=sst_3d(:,:,it)
           sst_1d(it)=Areal_Mean (X2d,Nlon_sst,Nlat_sst)
           !print*,sst_1d(it)
        Enddo 
        Deallocate (YY, XX, X1d,tlat,tlon,xtime)
        Deallocate (Mask_sst,Lon_sst,Lat_sst,sst_3d,X2d)
      End Subroutine Read_sst_data
      !---------------------------------------------------------------------------------
      Real Function Areal_Mean (X,Nlon,Nlat)
        Integer :: Nlon,Nlat, ilon,ilat
        Real :: X(Nlon,Nlat),lon1,lon2,lat1,lat2
        Real :: ic, tmp

        ic=0; tmp=0
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              if (X(ilon,ilat)/=Miss) then
                 ic=ic+1; tmp=tmp+X(ilon,ilat)
              endif
           Enddo
        Enddo
        Areal_Mean=tmp/ic
      End Function Areal_Mean
      !---------------------------------------------------------------------------------
      Subroutine Moving_Ave (X,N,m,Y)
         Integer :: N,m, i,j,k   ! m must be an odd number
         Real :: X(N), Y(N-m+1), tmp

         Do j=1,N-m+1
            tmp=0
            Do i=j,j+m-1
               tmp=tmp+X(i)
            Enddo
            Y(j)=tmp/Real(m)
         Enddo
      End Subroutine Moving_Ave
      !---------------------------------------------------------------------------------
      Logical Function Empty_Ar1(X,N)
         Integer :: N,i
         Real :: X(N)
         Logical :: Empty
         Empty=.False.
         Do i=1,N
            if (X(i)==Miss.OR.IsNaN(X(i)).OR.ABS(X(i))>1000) then
               Empty=.True.; Exit
            endif
         Enddo
         Empty_Ar1=Empty
      End Function Empty_Ar1

      !---------------------------------------------------------------------------------

      !----------------------------------------------------------------------------------!
      Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
!      Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90"

      END
      !----------------------------------------------------------------------------------!
      Include "/home/tanpv/Mylibs/Student_Prob.f90"