! Maximum Covariance Analysis of Drought Inde:wqx

      Use NetCDF
      Real, Parameter :: Miss=-99.
      Integer, Parameter :: Nmon=12, Yr1=1980, Yr2=2019, Nidx=24
      Integer, Parameter :: Nmode=6, Nreg=6
      !Real, Parameter :: Slon1=70, Slon2=280, Slat1=-20, Slat2=40
      Real, Parameter :: Slon1=110, Slon2=290, Slat1=-20, Slat2=20

      Real, Allocatable :: Lon(:), Lat(:), xtime(:)
      Real, Allocatable :: Lon_sst(:), Lat_sst(:), Mask_sst(:,:)
      Real, Allocatable :: Mask(:,:), xlon(:), xlat(:), ylon(:), ylat(:)
      Integer :: Landpoint, Nst, ist, Nlat, Nlon, ilat, ilon, Ntime,i,j
      Integer :: SeaPoint, Nsst, Nlon_sst, Nlat_sst
      Integer :: Yr_sta
      Real :: Miss_sst
      !---------------------------------------------------------------- 
      Real, Allocatable :: slp3d(:,:,:), us3d(:,:,:), u853d(:,:,:), u203d(:,:,:)
      Integer, parameter :: Nvar=3
      Real, Allocatable :: slpuv(:,:)
      Integer :: Nlon_slp, Nlat_slp, Npoint 
      Real, Allocatable :: lon_slp(:), lat_slp(:)
      !---------------------------------------------------------------- 

      ! For MCA:
      Real, Allocatable :: Xbar(:), Sx(:), Ybar(:), Sy(:)
      Real, Allocatable :: X(:,:), Y(:,:), Xc(:,:), Yc(:,:)
      Real, Allocatable :: C(:,:), U(:,:), V(:,:), WW(:,:), W(:)
      Real, Allocatable :: A(:,:), B(:,:)   ! Timeseries of X & Y
      Real, Allocatable :: Uk(:,:), Vk(:,:), WWk(:,:), Wk(:)
      Real, Allocatable :: Ak(:,:), Bk(:,:), Corre(:), Probs(:)
      Real, Allocatable :: Accum(:), Percent(:)         ! Accumulated percentage of Singular values
      Real, Allocatable :: XX(:), YY(:), Dr_IDX(:,:,:)  ! (Nlon,Nlat,Ntime)
      Real, Allocatable :: U3D(:,:,:),V3D(:,:,:)        ! Uk,Vk in 3D array
      Integer :: Nt, Nx, Ny, isys 
      Real, Allocatable :: Homo(:,:,:),Hete(:,:,:)      ! Homo and Hetero maps (Nlon,Nlat,Nmode)
      Real, Allocatable :: p_Ho(:,:,:),p_He(:,:,:)      ! Pval

      Integer :: Mov_Ave_Step                    ! Buoc truot trong Moving average                
      Real :: Sumx,Corr,prob

      Real, Allocatable :: xx1d(:), yy1d(:)
      Logical :: Calc_mtx, Normalize, SVD_MatLab
      !--------------------------------------------------------------------

      Character (Len=255) :: Fname, Fname1, Fnu,Fnv,Fnab,Dir_Out
      !Character (Len=15)  :: Dat_src="Ens" ! "APH" ! "Ens"
      Character (Len=15)  :: Dat_src="CRU" ! "Ens"
      Character (Len=15)  :: varn="spei_12"
      !--------------------------------------------------------------------
      ! LAPACK:
      Integer :: M, N, LDA, LDU, LDVT, LWMAX, INFO, LWORK
      DOUBLE PRECISION, Allocatable :: AA(:,:), UU(:,:), VT(:,:), SS(:), WORK(:)
      EXTERNAL ::        DGESVD
      !--------------------------------------------------------------------

      Mov_Ave_Step = 13           ! Smooth with moving average

      Ntime=(Yr2-Yr1+1)*Nmon

      Fname="/work/projects/G09/SEACLID/GPCC_CRU_SEA/"// &
            "SPEI_SPI_GPCC_CRU_59_19_New.nc"
      Call Read_SPEI (Fname)
      print*,"Reading dought index ... done ! Ntime,Landpoint=",Ntime,Landpoint
 
      Call Read_SLP_UV
      print*,"Reading SLP,U,V ... done ! Ntime, Npoint=",Ntime,Npoint

      Nt=Ntime; Nx=Landpoint; Ny=Npoint ! SeaPoint

      Allocate (Xc(Nx,Nt),Yc(Ny,Nt),Xbar(Nx),Ybar(Ny), Sx(Nx), Sy(Ny))
      !Allocate (C(Nx,Ny),U(Nx,Ny), W(Ny), V(Ny,Ny), A(Nt,Nx), B(Nt,Ny))
      Allocate (C(Nx,Ny),U(Nx,Ny), W(Ny), V(Ny,Ny), A(Nx,Nt), B(Ny,Nt))
      Allocate (WW(Ny,Ny))
      Allocate (Uk(Nx,Nmode), Vk(Ny,Nmode), WWk(Nmode,Nmode), Wk(Nmode))
      !Allocate (Ak(Nt,Nmode), Bk(Nt,Nmode), Corre(Nmode),Probs(Nmode))
      Allocate (Ak(Nmode,Nt), Bk(Nmode,Nt), Corre(Nmode),Probs(Nmode))
      Allocate (Accum(Nmode),Percent(Nmode))
      !Allocate (U3D(Nlon,Nlat,Nmode),V3D(Nlon_sst,Nlat_sst,Nmode))
      Allocate (U3D(Nlon,Nlat,Nmode),V3D(Nlon_slp,Nvar,Nmode))
      Allocate (XX(Nt),YY(Nt))
      !--------------------------------------------------------------------

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

      print*,"Centering data ..."
   
      !Normalize = .false.          ! Normalize data (true/false)
      Normalize = .True.
      !Calc_mtx=.false.              ! Calculate C-matrix (true/false)
      Calc_mtx=.true.              ! Calculate C-matrix (true/false)
      SVD_MatLab = .True.          ! Use SVD procedure from MatLab (true) or R (false)
      !--------------------------------------------------------------------

      Call Centered_Data (Normalize)            ! Su dung de tinh Ak, Bk

      if (Calc_mtx) then            ! Chi tinh Ma tran Covariance

         ! Covariance matrix:
         print*,"Calculating Covariance matrix..."
         !C = MatMul(Xc,Transpose(Yc))/Real(Nt-1)
         Call Calc_Cova
         print"(20F8.3)",C(Nx-20:Nx,Ny-20:Ny)

         ! Write C to nc file:
         Allocate (xx1d(Nx), yy1d(Ny))
         Do i=1,Nx; xx1d(i)=i; Enddo
         Do i=1,Ny; yy1d(i)=i; Enddo
         Call wr_ncdf_2d ("C_Matrix.nc",xx1d,yy1d,Nx,Ny,C,"cor", &
              "-", "-", "x", "y", Miss)
         Deallocate (xx1d,yy1d)

         ! Calculating SVD:
         if (SVD_MatLab) then         ! MatLab
            isys = system ("./run_mlab.sh ")
            Fnu="/home/tanpv/Progs/2021/NC/U.nc"
            Fnv="/home/tanpv/Progs/2021/NC/V.nc"
            Fnab="/home/tanpv/Progs/2021/NC/S.nc"
            Call Read_UVW_MatLab (Fnu,Fnv,Fnab)
         else         ! R
            isys = system ("./run_r.sh ")
            Fnu="/home/tanpv/Progs/2021/NC/u.nc"
            Fnv="/home/tanpv/Progs/2021/NC/v.nc"
            Fnab="/home/tanpv/Progs/2021/NC/w.nc"
            Call Read_UVW (Fnu,Fnv,Fnab)
         endif
         Call Calc_3D_AB

         !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      else                          ! Xu ky KQ sau khi thuc hien SVD bang MatLab hoac R

         if (SVD_MatLab) then         ! Use MatLab run
            Fnu="/home/tanpv/Progs/2021/NC/U.nc"
            Fnv="/home/tanpv/Progs/2021/NC/V.nc"
            Fnab="/home/tanpv/Progs/2021/NC/S.nc"
            Call Read_UVW_MatLab (Fnu,Fnv,Fnab)
         else                          ! Use R run
            Fnu="/home/tanpv/Progs/2021/NC/u.nc"
            Fnv="/home/tanpv/Progs/2021/NC/v.nc"
            Fnab="/home/tanpv/Progs/2021/NC/w.nc"
            Call Read_UVW (Fnu,Fnv,Fnab)
         endif
         Call Calc_3D_AB

      endif

      Dir_Out="/home/tanpv/Progs/2021/NC/"
      !--------------------------------------------------------------------
      ! Write to NC files:
      Fnu=trim(Dir_Out)//"U_Dr_Multi.nc"
      Fnv=trim(Dir_Out)//"V_Dr_Multi.nc"
      Fnab=trim(Dir_Out)//"AB_Dr_Multi.nc"
      Call Write_to_NC (Fnu,Fnv,Fnab)

      ! Calculate and write Homogenerous and Heterogenerous correlation maps: 
      !Call Homo_And_Hetero_Map
      !Fnu = trim(Dir_Out)//"Ho_He_SST.nc"
      !Call Write_Homo_Hete_Maps_to_NC (Fnu)
      !--------------------------------------------------------
      print*,"Wk", Wk
      print*,"Corr=", Corre
      print*,"Prob=",Probs
      print*,"Perc=",Percent
      print*,"Accum",Accum

      CONTAINS
      !----------------------------------------------------------------------------------------------
      Subroutine Calc_MYI
         ! Multi Year Index

      End Subroutine Calc_MYT
      !----------------------------------------------------------------------------------------------
      Subroutine Calc_Cova
         Integer :: i,j,k
         Real :: tmp

         Do j=1,Nx
            Do k=1,Ny
               tmp=0
               Do i=1,Nt
                  tmp=tmp+Xc(j,i)*Yc(k,i)
               Enddo
               C(j,k)=tmp/Real(Nt-1)
               !if (IsNaN(C(j,k))) then
               !   print*,Xc(j,:); pause
               !   print*,Yc(k,:); pause
               !   print*,j,k,C(j,k)
               !endif
            Enddo
         Enddo
      End Subroutine Calc_Cova
      !----------------------------------------------------------------------------------------------
      Subroutine Get_UVW (Fname,X,Nx,Ny,varn)
      Character (Len=*) :: Fname,varn
      Integer :: N1,N2,ncid, dimid, varid
      Real, Allocatable :: X(:,:)
      Logical :: status

      status = nf90_open(fname, nf90_NoWrite, ncid)
      status = nf90_inq_dimid(ncid,"x",dimid)
      status = nf90_inquire_dimension(ncid,dimid,len=Nx)
      status = nf90_inq_dimid(ncid,"y",dimid)
      status = nf90_inquire_dimension(ncid,dimid,len=Ny)

      Allocate (X(Nx,Ny))
      status = nf90_inq_varid(ncid,trim(varn),varid)
      status = nf90_get_var(ncid,varid,X)
      status =  nf90_close(ncid)
      End Subroutine Get_UVW

      !--------------------------------------------------------------------
      Subroutine Read_UVW_MatLab (Fu,Fv,Fw)
        Character (Len=*) :: Fu,Fv,Fw
        Real, Allocatable :: UU(:,:)
        Real, Allocatable :: xlo(:), xla(:)
        Integer :: nlon,nlat,i,j

        Call rd_ncdf_2d (Fw,xlo,xla,nlon,nlat,UU,"S")
        print*, nlon,nlat,nx,ny

        !Do i=1,Ny; W(i) = SQRT(UU(i,i)); Enddo
        Do i=1,Ny; W(i) = UU(i,i); Enddo

        Deallocate (xlo,xla,UU)

        Call rd_ncdf_2d (Fu,xlo,xla,nlon,nlat,UU,"U")
        print*, nlon,nlat,nx,ny
        Do i=1,Nx
           U(i,1:Ny) = UU(i,1:Ny) !*W(i)
        Enddo
        Deallocate (xlo,xla,UU)
        Call rd_ncdf_2d (Fv,xlo,xla,nlon,nlat,UU,"V")
        print*, nlon,nlat,nx,ny

        Do i=1,Ny
           V(i,1:Ny) = UU(i,1:Ny) !*W(i)
        Enddo
        Deallocate (xlo,xla,UU)
        !Call rd_ncdf_2d (Fw,xlo,xla,nlon,nlat,UU,"S")
        !print*, nlon,nlat,nx,ny

        !Do i=1,Ny; W(i) = UU(i,i); Enddo
        !Deallocate (xlo,xla,UU)
      End Subroutine Read_UVW_MatLab
      !--------------------------------------------------------------------
      Subroutine Read_UVW (Fu,Fv,Fw)
        Character (Len=*) :: Fu,Fv,Fw
        Real, Allocatable :: UU(:,:)
        Real, Allocatable :: xlo(:), xla(:)
        Integer :: nlon,nlat,i,j

        Call rd_ncdf_2d (Fu,xlo,xla,nlon,nlat,UU,"var")
        print*, nlon,nlat,nx,ny
        Do i=1,Nx
           U(i,1:Ny) = UU(i,1:Ny)
        Enddo
        Deallocate (xlo,xla,UU)
        Call rd_ncdf_2d (Fv,xlo,xla,nlon,nlat,UU,"var")
        print*, nlon,nlat,nx,ny

        Do i=1,Ny
           V(i,1:Ny) = UU(1:Ny,i)
        Enddo
        Deallocate (xlo,xla,UU)
        Call rd_ncdf_2d (Fw,xlo,xla,nlon,nlat,UU,"var")
        print*, nlon,nlat,nx,ny

        W(1:Ny) = UU(1:Ny,1)
        Deallocate (xlo,xla,UU)
      End Subroutine Read_UVW
      !--------------------------------------------------------------------
      Subroutine Calc_3D_AB
         Real :: Sumx
         Integer Nmin 
         Real, Allocatable :: X1d(:), Y1d(:)
         ! Sign inversion of U and V:
         !U = -U; V = -V

         print*,"Nx,Ny,Nmode: ",Nx,Ny,Nmode
         Wk(1:Nmode) = W(1:Nmode)                  ! Singular values
         Do i=1,Nmode
            Uk(:,i)=U(:,i); Vk(:,i)=V(:,i)
         Enddo

         ! Singular matrix:
         WW=0;  Do j=1,Ny; WW(j,j)=W(j); Enddo      ! Full 
         WWk=0; Do j=1,Nmode; WWk(j,j)=W(j); Enddo  ! Nmode


         ! Timeseries of Nmode
         Ak = MatMul(Transpose(Uk),Xc)
         Bk = MatMul(Transpose(Vk),Yc)

         ! Correlations:
         Allocate (X1d(Nt), Y1d(Nt))
         Do i=1,Nmode
            X1d(1:Nt)=Ak(i,1:Nt); Y1d(1:Nt)=Bk(i,1:Nt)
            Call Correlation (X1d,Y1d,Nt,Corr,Prob, Miss)
            Corre(i)=Corr; Probs(i)=prob
         Enddo
         Deallocate (X1d,Y1d)
         ! Accum (%) of singular values:
         Nmin=Min(Nx,Ny)
         !print*,"(30F6.1)",W(1:Nmin)
         Sumx=Sum(W(1:Nmin)); Accum(:)=0
         Do i=1,Nmode
            Percent(i)=W(i)/Sumx*100
            Accum(i)=Sum(Percent(1:i))
         Enddo
         print*,"Var: ",Percent(:)
         print*,"Cum: ",Accum(:)
         ! Convert Uk to U3D:
         U3D=Miss
         Do i=1,Nmode
            j=0
            Do ilat=1,Nlat
               Do ilon=1,Nlon
                  if (Mask(ilon,ilat)==1) then
                     j=j+1
                     U3D(ilon,ilat,i) = Uk(j,i)
                  endif
               Enddo
            Enddo
         Enddo

         ! Convert Vk to V3D:
         V3D=Miss
         Do i=1,Nmode
            j=0
            Do ilat=1,Nvar !Nlat_sst
               Do ilon=1,Nlon_slp ! Nlon_sst
                  !if (Mask_sst(ilon,ilat)==1) then
                     j=j+1
                     V3D(ilon,ilat,i) = Vk(j,i)
                  !endif
               Enddo
            Enddo
         Enddo
  
      End Subroutine Calc_3D_AB
      !--------------------------------------------------------------------
      Subroutine Read_SLP_UV
        Real, Parameter :: Miss_slp=-9.96921e+36 ! U: -9.96921e+36
        Integer, Parameter :: Yr_sta=1948, Yr_sta_olr=1975
        Character (Len=255) :: Fname
        Real, Allocatable :: X3d(:,:,:), X4d(:,:,:,:)
        Real, Allocatable :: X1d(:), Y1d(:), tlon(:), tlat(:), tlev(:), xtime(:)
        Integer :: i,j, ilon,ilat, nlon,nlat,nlev, nt
        Integer :: it1, it2, ilon1, ilon2,ilat1,ilat2, imon, M, Md2
        Real :: ic, tmp

        !Fname="/work/projects/G09/SEACLID/NCEP/slp.mon.mean.nc"
        !print*, "Reading SLP ..."
        !Call rd_ncdf_3d (Fname,tlon,tlat,xtime,nlon,nlat,nt,X3d,"slp")

        Fname="/work/projects/G09/SEACLID/NCEP/olr.mon.mean_74_21.nc"
        Call rd_ncdf_3d (Fname,tlon,tlat,xtime,nlon,nlat,nt,X3d,"olr") !trim(varn))

        Where (X3d==Miss_slp) X3d=Miss
        Where (IsNaN(X3d)); X3d=Miss; Endwhere

        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

        !it1=(Yr1-Yr_sta)*Nmon+1; it2=it1+Ntime-1
        it1=(Yr1-Yr_sta_olr)*Nmon+8; it2=it1+Ntime-1        ! Since 6-1974

        Nlon_slp=nlon; Nlat_slp=ilat1-ilat2+1
        print*,"Check: Nlon SLP, Slat1,Slat2: ",Nlon,tlat(ilat2),tlat(ilat1)

        Allocate (slp3d(Nlon_slp,Nlat_slp,Ntime))

        Do ilon=1,Nlon_slp
           j=Nlat_slp+1
           Do ilat=1,Nlat
              if (ilat>=ilat2.AND.ilat<=ilat1) then
                 j=j-1
                 slp3d(ilon,j,1:Ntime) = X3d(ilon,ilat,it1:it2)
                 !print*,slp3d(ilon,j,1:Ntime),nt,Ntime
              endif
           Enddo
        Enddo
        !print*,slp3d, Nlon_slp,Nlat_slp,j,Ntime; pause
        !
        Fname ="/work/projects/G09/SEACLID/NCEP/uwnd_srf.mon.mean.nc"
        Deallocate (tlon,tlat,xtime,X3d)
        Call rd_ncdf_3d (Fname,tlon,tlat,xtime,nlon,nlat,nt,X3d,"uwnd") 
        Where (X3d==Miss_slp) X3d=Miss
        Where (IsNaN(X3d)); X3d=Miss; Endwhere

        print*,"Check: Nlon u_srf: ",Nlon

        it1=(Yr1-Yr_sta)*Nmon+1; it2=it1+Ntime-1    ! Diff to OLR
        Allocate (us3d(Nlon_slp,Nlat_slp,Ntime))

        Do ilon=1,Nlon_slp
           j=Nlat_slp+1
           Do ilat=1,Nlat
              if (ilat>=ilat2.AND.ilat<=ilat1) then
                 j=j-1
                 us3d(ilon,j,1:Ntime) = X3d(ilon,ilat,it1:it2)
                 !print*,us3d(ilon,j,1:Ntime), ilon,j,Ntime; ! pause
              endif
           Enddo
        Enddo
        !pause
        !
        Fname="/work/projects/G09/SEACLID/NCEP/uwnd.mon.mean.nc"
        !level = 1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, 50, 
        !    30, 20, 10 ;
        Deallocate (tlon,tlat,xtime,X3d)
        Call rd_ncdf_4d (Fname,tlon,tlat,tlev,xtime,nlon,nlat,nlev,nt,X4d,"uwnd")
        Where (X4d==Miss_slp) X4d=Miss
        Where (IsNaN(X4d)); X4d=Miss; Endwhere
        print*,"Check: Nlon u: ",Nlon

        Allocate (u853d(Nlon_slp,Nlat_slp,Ntime),u203d(Nlon_slp,Nlat_slp,Ntime))

        Do ilon=1,Nlon_slp
           j=Nlat_slp+1
           Do ilat=1,Nlat
              if (ilat>=ilat2.AND.ilat<=ilat1) then
                 j=j-1; 
                 u853d(ilon,j,1:Ntime) = X4d(ilon,ilat,3,it1:it2)
                 u203d(ilon,j,1:Ntime) = X4d(ilon,ilat,10,it1:it2)
                 !print*,u203d(ilon,j,1:Ntime), ilon,j,Ntime
              endif
           Enddo
        Enddo
        !pause
        print*,"it1,it2,ilat1,ilat2,Nlat_slp=",it1,it2,ilat1,ilat2,Nlat_slp; ! pause
        if (Nlon_slp/=Nlon) then; print*,"check??????"; pause; endif
        ! Combine:
        Deallocate (tlon,tlat,tlev,xtime,X4d)
        !!!!!!!!
        Allocate (X4d(Nlon_slp,Nlat_slp,Ntime,Nvar))
        X4d(1:Nlon_slp,1:Nlat_slp,1:Ntime,1)=slp3d(1:Nlon_slp,1:Nlat_slp,1:Ntime)
        X4d(1:Nlon_slp,1:Nlat_slp,1:Ntime,2)=u853d(1:Nlon_slp,1:Nlat_slp,1:Ntime)
        X4d(1:Nlon_slp,1:Nlat_slp,1:Ntime,3)=u203d(1:Nlon_slp,1:Nlat_slp,1:Ntime)
        if (Nvar==4) &
        X4d(1:Nlon_slp,1:Nlat_slp,1:Ntime,4)=us3d(1:Nlon_slp,1:Nlat_slp,1:Ntime)

        Npoint=Nvar*Nlon_slp

        Allocate (slpuv(Npoint,Ntime),X1d(Nlat_slp))
        slpuv=Miss
        i=0
        Do ivar=1,Nvar
           !print*,X4d(:,:,:,ivar); pause
           Do ilon=1,Nlon_slp
              i=i+1
              Do it=1,Ntime
                 ic=0; tmp=0
                 Do ilat=1,Nlat_slp
                    if (X4d(ilon,ilat,it,ivar)/=Miss) then
                       ic=ic+1; tmp=tmp+X4d(ilon,ilat,it,ivar)
                    endif
                 Enddo
                 if (ic>0) slpuv(i,it)=tmp/ic
              Enddo
              !print*,slpuv(i,:),i
           Enddo
        Enddo

        ! Anomaly:
        Allocate (Y(Npoint,Ntime))
        !print*,"Check: Nlon, Npoint: ",Nlon, Npoint
        !Deallocate (X1d); Allocate (X1d(Ntime))
        !Do i=1,Npoint
        !   X1d(1:Ntime)=slpuv(i,1:Ntime)
        !   if (.not.Ar1_Empty(X1d(1:Ntime),Ntime, Miss)) then
        !      Call Anomaly (X1d,Ntime)
        !      slpuv(i,1:Ntime)=X1d(1:Ntime)
        !      !print*,slpuv(i,1:Ntime)
        !   endif
        !Enddo
        Y(:,:)=slpuv(:,:)
      End Subroutine Read_SLP_UV
      !--------------------------------------------------------------------
      
      !--------------------------------------------------------------------
      Subroutine Homo_And_Hetero_Map
        ! Calculate Homogenerous correlation map = r[X(t), a_k (t)]
        ! and     Heterogenerous correlation map = r[X(t), b_k (t)]

        Integer :: ilon,ilat,it,i,j
        Real :: ctmp, ptmp
        Real, Allocatable :: Ho(:,:),He(:,:),pHo(:,:),pHe(:,:),X1d(:), A1d(:), B1d(:)

        Allocate (X1d(Nt),A1d(Nt),B1d(Nt),Ho(Nx,Nmode),He(Nx,Nmode))
        Allocate (pHo(Nx,Nmode),pHe(Nx,Nmode))

        Do i=1,Nmode
           A1d(1:Nt)=Ak(i,1:Nt); B1d(1:Nt)=Bk(i,1:Nt)
           Do j=1,Nx
              X1d(1:Nt)=X(j,1:Nt)
              Call Correlation(X1d,A1d,Nt,ctmp,ptmp,Miss)
              Ho(j,i)=ctmp; pHo(j,i)=ptmp
              Call Correlation(X1d,B1d,Nt,ctmp,ptmp,Miss)
              He(j,i)=ctmp; pHe(j,i)=ptmp
           Enddo
        Enddo
        ! Convert to 3D:
        Allocate (Homo(Nlon,Nlat,Nmode),p_Ho(Nlon,Nlat,Nmode))
        Allocate (Hete(Nlon,Nlat,Nmode),p_He(Nlon,Nlat,Nmode))
         Homo=Miss; Hete=Miss; p_Ho=Miss; p_He=Miss
         Do i=1,Nmode
            j=0
            Do ilat=1,Nlat
               Do ilon=1,Nlon
                  if (Mask(ilon,ilat)==1) then
                     j=j+1
                     Homo(ilon,ilat,i)=Ho(j,i); p_Ho(ilon,ilat,i)=pHo(j,i)
                     Hete(ilon,ilat,i)=He(j,i); p_He(ilon,ilat,i)=pHe(j,i)
                  endif
               Enddo
            Enddo
         Enddo

      End Subroutine Homo_And_Hetero_Map
      !--------------------------------------------------------------------
      Subroutine Write_Homo_Hete_Maps_to_NC (Fnu)
        Character (Len=*) :: Fnu
        Real, Allocatable :: xmod(:)
        Character (Len=150) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm
        Integer :: i,j,k

        Allocate (xmod(Nmode))
        Do i=1,Nmode; xmod(i)=i; Enddo
        t_unit="months since 1962-01-01 "
        var_ln="Homogenerous correlation map r[Dr(t), a_k (t)]"
        Call wr_ncdf_3d (Fnu,xlon,xlat,xmod,nlon,nlat,Nmode,Homo,"ho", &
             t_unit, var_ln, "-", "lon", "lat", "mod", Miss)
        var_ln="p-values of Homo"
        Call wr_ncdf_3d_apd (Fnu,p_Ho,"pho",nlon,nlat,Nmode, var_ln, &
             "-", "lon", "lat", "mod", Miss)
        var_ln="Heterogenerous correlation map r[Dr(t), b_k (t)]"
        Call wr_ncdf_3d_apd (Fnu,Hete,"he",nlon,nlat,Nmode, var_ln, &
             "-", "lon", "lat", "mod", Miss)
        var_ln="p-values of Hete"
        Call wr_ncdf_3d_apd (Fnu,p_He,"phe",nlon,nlat,Nmode, var_ln, &
             "-", "lon", "lat", "mod", Miss)
      End Subroutine Write_Homo_Hete_Maps_to_NC
      !--------------------------------------------------------------------
      Subroutine Write_to_NC (Fnu,Fnv,Fnab)
        Character (Len=*) :: Fnu, Fnv, Fnab  ! Filename of U, V and A,B
        Real, Allocatable :: xmod(:), X2D(:,:), xtime(:)
        Real, Allocatable :: xvar(:), xlon1(:), X3d(:,:,:)
        Character (Len=150) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm
        Integer :: i,j,k

        print*,"Check in Write_to_NC:"
        print*,"Nx,Ny,Nt,Nvar,nlon,nlat,Nlon_slp,Nmode:"
        print*,Nx,Ny,Nt,Nvar,nlon,nlat,Nlon_slp,Nmode
        Allocate (xmod(Nmode),X2D(Nmode,Ny),xtime(Nt))

        Do i=1,Nmode; xmod(i)=i; Enddo
        Do i=1,Nt; xtime(i)=i; Enddo

        Allocate (xvar(Nvar),xlon1(Nlon_slp))
        Do i=1,Nvar; xvar(i)=i; Enddo
        Do i=1,Nlon_slp; xlon1(i)=(i-1)*2.5; Enddo

        ! U matrix:
        t_unit="months since 1962-01-01 "
        var_ln="MCA for drought index and SST - U-matrix"
        Call wr_ncdf_3d (Fnu,xlon,xlat,xmod,nlon,nlat,Nmode,U3D,"u", &
             t_unit, var_ln, "-", "lon", "lat", "mod", Miss)

        ! V matrix:
        !Allocate (X3d(Nlon_slp,Nvar,Nmode))
        !Do i=1,Nvar
        !   Do j=1,Nlon_slp
        !      X3d(j,i,1:Nmode)=V3D(i,j,1:Nmode)
        !   Enddo
        !Enddo
        var_ln="MCA for drought index and Mutivar - V-matrix:"// &
               " 4: slp, 1: u_srf, 2: u850, 3: u200"
        Call wr_ncdf_3d (Fnv,xlon1,xvar,xmod,Nlon_slp,Nvar,Nmode,V3D,"v", &
             t_unit, var_ln, "-", "lon", "var", "mod", Miss)

        !var_ln="Percentage of contribution of singular values"
        !Call wr_ncdf_1d_apd (Fnv,Percent,"per", Nmode, var_ln,"-","mod", Miss)
        !var_ln="Accumulated Percent. of contribution of singular values"
        !Call wr_ncdf_1d_apd (Fnv,Accum,"cum", Nmode, var_ln,"-","mod", Miss)

        ! A & B matrix
        var_ln="Timeseries of MCA of drought index (Ak)"
        Call wr_ncdf_2d (Fnab,xmod,xtime,Nmode,Nt,Ak, "a", var_ln, &
                         "-", "mod", "mon", Miss)
        var_ln="Timeseries of MCA of SST (Bk)"
        Call wr_ncdf_2d_apd (Fnab,Bk,"b",Nmode,Nt, var_ln, "-", &
                         "mod", "mon", Miss)
        var_ln="Correlation coefficients between A & B"
        Call wr_ncdf_1d_apd (Fnab,Corre,"cor", Nmode, var_ln,"-","mod", Miss)
        var_ln="p-values of Correlation coefficients between A & B"
        Call wr_ncdf_1d_apd (Fnab,probs,"prob", Nmode, var_ln,"-","mod", Miss)

        var_ln="Percentage of contribution of singular values"
        Call wr_ncdf_1d_apd (Fnab,Percent,"per", Nmode, var_ln,"-","mod", Miss)
        var_ln="Accumulated Percent. of contribution of singular values"
        Call wr_ncdf_1d_apd (Fnab,Accum,"cum", Nmode, var_ln,"-","mod", Miss)

      End Subroutine Write_to_NC
      !--------------------------------------------------------------------
      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

        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 Centered_Data (Normalize)
        Integer :: i,j
        Logical :: Normalize
        Real, Allocatable :: Std(:), X1d(:)
       
        ! Centering:
        print*,"Centering: Nx, Ny, Nt: ",Nx, Ny, Nt
        Allocate (Std(Nx), X1d(Nt))
        Do j=1,Nx
           X1d(1:Nt)=X(j,1:Nt)
           Xbar(j)=Ave    (X1d,Nt,Miss)
           Std (j)=Stdev_1(X1d,Nt,Miss)
           Do i=1,Nt
              Xc(j,i)=X(j,i)-Xbar(j)
              if (Normalize) Xc(j,i)=Xc(j,i)/Std(j) !/Stdev_1(X(j,1:Nt),Nt,Miss)
           Enddo
        Enddo
        Deallocate (Std); Allocate (Std(Ny))
        Do j=1,Ny
           X1d(1:Nt)=Y(j,1:Nt)
           Ybar(j)=Ave    (X1d,Nt,Miss)
           Std (j)=Stdev_1(X1d,Nt,Miss)
           Do i=1,Nt
              Yc(j,i)=Y(j,i)-Ybar(j)
              if (Normalize) Yc(j,i)=Yc(j,i)/Std(j) ! Stdev_1(Y(j,1:Nt),Nt,Miss)
           Enddo
        Enddo
      End Subroutine Centered_Data
      !--------------------------------------------------------------------
      Subroutine Anomaly (X,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 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
      !----------------------------------------------------------------------------------------------
      Subroutine Read_SPEI (Fname)
      Character (Len=*) :: Fname
      Integer, parameter :: Yrs_Start=1959
      Integer :: ilon,ilat,it,it1, Nt, Nt1, ix, it_skp
      Real, Allocatable :: spe3d(:,:,:), X3d(:,:,:), Y3d(:,:,:), xtime(:), x1d(:)
      

      !Lon_spe(:), Lat_spe(:)
      !Call rd_ncdf_3d (Fname,lon_spe,lat_spe,xtime,nlon_spe,nlat_spe,nt,X3d,"spei12")
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlon,nlat,Nt,X3d,"spei12")
      !Deallocate (lon,lat,xtime)
      !Call rd_ncdf_2d (Fname,lon,lat,nlon,nlat,Mask,"land")
      Where (IsNaN(X3d)) X3d=Miss
      it_skp=(Yr1-Yrs_Start)*Nmon     ! In case Yr1>Yrs_Start
      it1=it_skp+1; Nt1=Nt-it_skp

      !it1=Nmon+1; Nt1=Nt-Nmon

      if (Nt1/=Ntime) then; print*,"Check: Nt,Ntime:",Nt1,Ntime; stop; endif

      Allocate (spe3d(nlon,nlat,Nt1),mask(nlon,nlat))
      Allocate (x1d(Nt1))
      spe=Miss; mask=Miss; Landpoint=0
      Do ilon=1,nlon
         Do ilat=1,nlat
            x1d(1:Nt1)=X3d(ilon,ilat,it1:Nt)
            if (.not.Ar1_Empty (X1d,Nt1,Miss)) then
               mask(ilon,ilat)=1
               Landpoint=Landpoint+1
               spe3d(ilon,ilat,1:Nt1)=X3d(ilon,ilat,it1:Nt)
            endif
         Enddo
      Enddo
      ! Reshape to 2D arr: Ntime=Nt1
      Allocate (X(Landpoint,Ntime))
      !Allocate (spe2d(LandSpe,Ntime))
      ix=0
      Do ilat=1,Nlat
         Do ilon=1,Nlon
            if (mask(ilon,ilat)/=Miss) then
               ix=ix+1; X(ix,1:Ntime)=spe3d(ilon,ilat,1:Ntime)
            endif
         Enddo
      Enddo
     End Subroutine Read_SPEI
     !----------------------------------------------------------------------------------------------
      Logical Function Ar1_Empty(X,N, Miss)
        Integer :: N, i
        Real :: X(N), Miss
        Logical :: Empty
        Empty=.false.
        Do i=1,N
           if (X(i)==Miss) then
              Empty=.true.; exit
           endif
        Enddo
        Ar1_Empty=Empty
      End Function Ar1_Empty
      !--------------------------------------------------------------------
      Include "/home/tanpv/Mylibs/Sub_Trend_MK_Linear.f90"
      Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
      Include "/home/tanpv/Mylibs/Ave.f90"
      Include "/home/tanpv/Mylibs/Stdev_1.f90"
      !Include "Sub_MCA.f90"
      !--------------------------------------------------------------------
      END
      Include "/home/tanpv/Mylibs/Student_Prob.f90"