! 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(:)                    ! (Ntime)
      !---------------------------------------------------------------------------------
      ! Composite maps of SLP:
      Real, allocatable :: lon_slp(:), lat_slp(:), Mask_slp(:,:)
      Integer :: Nlon_slp, Nlat_slp
      Real, allocatable :: slp_dry(:,:), slp_wet(:,:), slp_nor(:,:)
      Real, allocatable :: slp_cli(:,:)
      Real, allocatable :: slp_3d(:,:,:)
      !---------------------------------------------------------------------------------
      ! UV_srf:
      Real, allocatable :: us_dry(:,:),vs_dry(:,:)
      Real, allocatable :: us_wet(:,:),vs_wet(:,:)
      Real, allocatable :: us_nor(:,:),vs_nor(:,:)
      Real, allocatable :: us_cli(:,:),vs_cli(:,:)
      Real, allocatable :: us_3d(:,:,:), vs_3d(:,:,:) 
      Real, allocatable :: lon_uvs(:), lat_uvs(:), Mask_uvs(:,:)
      Integer :: Nlon_uvs, Nlat_uvs
      !---------------------------------------------------------------------------------
      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 
      !---------------------------------------------------------------------------------
      ! Composite SST maps:
      !---------------------------------------------------------------------------------
      ! SST info:
      Yr_sta     = 1891
      Miss_sst   = 1.e+20         ! Missing values
      Slon1      = 110
      Slon2      = 290
      Slat1      = -30
      !Slat2      =  30
      Slat2      =  50

      ! SST:
      Fname="/work/projects/G09/SEACLID/SST/sst.mon.mean_COBE.nc"
      Call Read_sst_data (Fname)
      Fname="/home/tanpv/Progs/2021/NC/SST_composite_maps_Large.nc"
      !!Fname="/home/tanpv/Progs/2021/NC/SST_composite_maps_Small.nc"
      Call SST_Composite_maps (Fname)

      ! SLP:
      Fname="/work/projects/G09/SEACLID/NCEP/slp.mon.mean.nc"
      Call Read_SLP (Fname)
      Fname="/home/tanpv/Progs/2021/NC/SLP_composite_maps_Large.nc"
      Call SLP_Composite_maps (Fname)

      ! UV-srf:
      Fname ="/work/projects/G09/SEACLID/NCEP/uwnd_srf.mon.mean.nc"
      Fname1="/work/projects/G09/SEACLID/NCEP/vwnd_srf.mon.mean.nc"
      Call Read_UV_Srf (Fname,Fname1)
      Fname="/home/tanpv/Progs/2021/NC/UV_srf_composite_maps_Large.nc"
      Call UV_srf_Composite_maps (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 SLP_Composite_maps (Fout)
        Character (Len=*) :: Fout
        Character (Len=255) :: Fname, Dir
        Real, Allocatable :: X2d(:,:), xreg(:),xscal(:), time(:), X3d(:,:,:)
        Real, Allocatable :: Pval_slp(:,:)
        Real, parameter :: Dr_thres=-0.5, W_thres=0.5
        Integer :: Nregs, Nt, nscal
        Integer :: ist, i, j, ilon, ilat, ipoint, it, nlon, nlat
        Integer :: ic,it1, it2, ilon1, ilon2,ilat1,ilat2, imon, M, Md2
        Real :: T_Dist, tmp1, tmp2, ntmp1, ntmp2, Sx, Sy
        Real :: tmp3, tmp4, ntmp3, ntmp4

        Allocate (slp_dry(Nlon_slp,Nlon_slp),slp_wet(Nlon_slp,Nlon_slp))
        Allocate (slp_nor(Nlon_slp,Nlon_slp),slp_cli(Nlon_slp,Nlon_slp))
        Allocate (Pval_slp(Nlon_slp,Nlon_slp))
        slp_dry=Miss; slp_wet=Miss; slp_nor=Miss; slp_cli=Miss; Pval_slp=Miss

        ic=0
        Do ilon=1,Nlon_slp
           Do ilat=1,Nlat_slp
               if (.not.Empty_Ar1(slp_3d(ilon,ilat,1:Ntime),Ntime)) then
                  tmp1=0; tmp2=0; ntmp1=0; ntmp2=0; Sx=0; Sy=0
                  tmp3=0; tmp4=0; ntmp3=0; ntmp4=0; ic=ic+1
                  Do it=1,Ntime
                     ntmp4=ntmp4+1; tmp4=tmp4+slp_3d(ilon,ilat,it)    ! Climatology
                     if (Areal_mean_Dr(it)<=Dr_thres) then            ! Dry
                        ntmp1=ntmp1+1
                        tmp1=tmp1+slp_3d(ilon,ilat,it)
                        Sx=Sx+slp_3d(ilon,ilat,it)*slp_3d(ilon,ilat,it)
                     else if (Areal_mean_Dr(it)>=W_thres) then        ! Wet
                        ntmp2=ntmp2+1
                        tmp2=tmp2+slp_3d(ilon,ilat,it)
                        Sy=Sy+slp_3d(ilon,ilat,it)*slp_3d(ilon,ilat,it)
                     else if (Areal_mean_Dr(it)>Dr_thres.AND.Areal_mean_Dr(it)<W_thres) then
                        ntmp3=ntmp3+1
                        tmp3=tmp3+slp_3d(ilon,ilat,it)
                     endif
                  Enddo
                  tmp1=tmp1/ntmp1; tmp2=tmp2/ntmp2
                  Sx=Sx/(ntmp1-1); Sy=Sy/(ntmp2-1)
                  Sx=SQRT(Sx-tmp1*tmp1); Sy=SQRT(Sy-tmp2*tmp2)
                  slp_dry(ilon,ilat)=tmp1
                  slp_wet(ilon,ilat)=tmp2
                  Pval_slp(ilon,ilat)=P_value(tmp1, tmp2, Sx, Sy, ntmp1, ntmp2)
                  slp_nor(ilon,ilat)=tmp3/ntmp3
                  slp_cli(ilon,ilat)=tmp4/ntmp4
                  !print*,ntmp1,ntmp2,ntmp3,ntmp4
               endif
           Enddo
        Enddo
        Call wr_ncdf_2d (Fout,Lon_slp,Lat_slp, Nlon_slp,Nlat_slp,slp_dry, &
                     "slp_dry","Mean SLP in case Dr_idx<-0.5", &
                     "mb", "lon","lat", Miss)
        Call wr_ncdf_2d_apd (Fout,slp_wet, "slp_wet", Nlon_slp,Nlat_slp, &
                 "Mean SLP in case Dr_idx>0.5", "mb", "lon", "lat", Miss)
        Call wr_ncdf_2d_apd (Fout,Pval_slp,"pval",Nlon_slp,Nlat_slp, &
                   "P-value of testing if slp_dry /= slp_wet","-","lon", "lat", Miss)
        Call wr_ncdf_2d_apd (Fout,slp_nor, "slp_nor", Nlon_slp,Nlat_slp, &
                 "Mean SLP in -0.5<case Dr_idx<0.5", "mb", "lon", "lat", Miss)
        Call wr_ncdf_2d_apd (Fout,slp_cli, "slp_cli", Nlon_slp,Nlat_slp, &
                 "Climatological Mean SLP", "mb", "lon", "lat", Miss)
      End Subroutine SLP_Composite_maps
      !--------------------------------------------------------------------------------------
      Subroutine UV_srf_Composite_maps (Fout)
        Character (Len=*) :: Fout
        Character (Len=255) :: Fname, Dir
        Real, Allocatable :: X2d(:,:), xreg(:),xscal(:), time(:), X3d(:,:,:)
        Real, Allocatable :: Pval_slp(:,:)
        Real, parameter :: Dr_thres=-0.5, W_thres=0.5
        Integer :: Nregs, Nt, nscal
        Integer :: ist, i, j, ilon, ilat, ipoint, it, nlon, nlat
        Integer :: ic,it1, it2, ilon1, ilon2,ilat1,ilat2, imon, M, Md2
        Real :: T_Dist, tmp11, tmp12,tmp21, tmp22, ntmp1, ntmp2, Sx, Sy
        Real :: tmp31, tmp41,tmp32, tmp42, ntmp3, ntmp4

        Allocate (us_dry(Nlon_uvs,Nlon_uvs),us_wet(Nlon_uvs,Nlon_uvs))
        Allocate (us_nor(Nlon_uvs,Nlon_uvs),us_cli(Nlon_uvs,Nlon_uvs))
        Allocate (vs_dry(Nlon_uvs,Nlon_uvs),vs_wet(Nlon_uvs,Nlon_uvs))
        Allocate (vs_nor(Nlon_uvs,Nlon_uvs),vs_cli(Nlon_uvs,Nlon_uvs))
!        Allocate (Pval_slp(Nlon_slp,Nlon_slp))
        us_dry=Miss; us_wet=Miss; us_nor=Miss; us_cli=Miss;
        vs_dry=Miss; vs_wet=Miss; vs_nor=Miss; vs_cli=Miss;

        ic=0
        Do ilon=1,Nlon_uvs
           Do ilat=1,Nlat_uvs
               if (.not.Empty_Ar1(slp_3d(ilon,ilat,1:Ntime),Ntime)) then
                  tmp11=0; tmp12=0;tmp21=0; tmp22=0; ntmp1=0; ntmp2=0; Sx=0; Sy=0
                  tmp31=0; tmp32=0; tmp41=0; tmp42=0; ntmp3=0; ntmp4=0; ic=ic+1
                  Do it=1,Ntime
                     ntmp4=ntmp4+1; 
                     tmp41=tmp41+us_3d(ilon,ilat,it)    ! Climatology
                     tmp42=tmp42+vs_3d(ilon,ilat,it)
                     if (Areal_mean_Dr(it)<=Dr_thres) then            ! Dry
                        ntmp1=ntmp1+1
                        tmp11=tmp11+us_3d(ilon,ilat,it)
                        tmp12=tmp12+vs_3d(ilon,ilat,it)
                     else if (Areal_mean_Dr(it)>=W_thres) then        ! Wet
                        ntmp2=ntmp2+1
                        tmp21=tmp21+us_3d(ilon,ilat,it)
                        tmp22=tmp22+vs_3d(ilon,ilat,it)
                     else if (Areal_mean_Dr(it)>Dr_thres.AND.Areal_mean_Dr(it)<W_thres) then
                        ntmp3=ntmp3+1
                        tmp31=tmp31+us_3d(ilon,ilat,it)
                        tmp32=tmp32+vs_3d(ilon,ilat,it)
                     endif
                  Enddo
                  us_dry(ilon,ilat)=tmp11/ntmp1
                  vs_dry(ilon,ilat)=tmp12/ntmp1
                  us_wet(ilon,ilat)=tmp21/ntmp2
                  vs_wet(ilon,ilat)=tmp22/ntmp2
                  us_nor(ilon,ilat)=tmp31/ntmp3
                  vs_nor(ilon,ilat)=tmp32/ntmp3
                  us_cli(ilon,ilat)=tmp41/ntmp4
                  vs_cli(ilon,ilat)=tmp42/ntmp4
                  !print*,ntmp1,ntmp2,ntmp3,ntmp4
               endif
           Enddo
        Enddo
        Call wr_ncdf_2d (Fout,Lon_uvs,Lat_uvs, Nlon_uvs,Nlat_uvs,us_dry, &
                     "u_dry","U-srf in case Dr_idx<-0.5", &
                     "m/s", "lon","lat", Miss)
        Call wr_ncdf_2d_apd (Fout,us_wet, "u_wet", Nlon_uvs,Nlat_uvs, &
                 "U-srf in case Dr_idx>0.5", "m/s", "lon", "lat", Miss)
        Call wr_ncdf_2d_apd (Fout,vs_dry, "v_dry", Nlon_uvs,Nlat_uvs, &
                 "V-srf in case Dr_idx>0.5", "m/s", "lon", "lat", Miss)
        Call wr_ncdf_2d_apd (Fout,vs_wet, "v_wet", Nlon_uvs,Nlat_uvs, &
                 "V-srf in case Dr_idx>0.5", "m/s", "lon", "lat", Miss)

        Call wr_ncdf_2d_apd (Fout,us_nor, "u_nor", Nlon_uvs,Nlat_uvs, &
                 "U-srf in -0.5<case Dr_idx<0.5", "m/s", "lon", "lat", Miss)
        Call wr_ncdf_2d_apd (Fout,vs_nor, "v_nor", Nlon_uvs,Nlat_uvs, &
                 "V-srf in -0.5<case Dr_idx<0.5", "m/s", "lon", "lat", Miss)

        Call wr_ncdf_2d_apd (Fout,us_cli, "u_cli", Nlon_uvs,Nlat_uvs, &
                 "Climatological U-srf", "m/s", "lon", "lat", Miss)
        Call wr_ncdf_2d_apd (Fout,vs_cli, "v_cli", Nlon_uvs,Nlat_uvs, &
                 "Climatological V-srf", "m/s", "lon", "lat", Miss)
      End Subroutine UV_srf_Composite_maps
      !--------------------------------------------------------------------------------------
      Subroutine SST_Composite_maps (Fout)
        Character (Len=*) :: Fout
        Character (Len=255) :: Fname, Dir
        Real, Allocatable :: X2d(:,:), xreg(:),xscal(:), time(:), X3d(:,:,:)
        Real, Allocatable :: Pval(:,:)
        Real, parameter :: Dr_thres=-0.5, W_thres=0.5
        Integer :: Nregs, Nt, nscal
        Integer :: ist, i, j, ilon, ilat, ipoint, it, nlon, nlat
        Integer :: ic,it1, it2, ilon1, ilon2,ilat1,ilat2, imon, M, Md2
        Real :: T_Dist, tmp1, tmp2, ntmp1, ntmp2, Sx, Sy
        Real :: tmp3, tmp4, ntmp3, ntmp4

!      Real, allocatable :: sst_dry(:,:), sst_wet(:,:)   ! (Nlon_sst, Nlat_sst): SST with Areal mean of SPEI<-0.5 and >-0.5
!      Real, allocatable :: Areal_mean_Dr(:)                    ! (Ntime)

        Dir="/work/projects/G09/SEACLID/GPCC_CRU_SEA/"
        Fname=trim(Dir)//"SPEI_SPI_GPCC_CRU_Reg_59_19_New.nc"
        Call rd_ncdf_3d (Fname, xreg, xscal, time, Nregs, nscal, Nt,X3d, "spe12")
        !Call rd_ncdf_2d (Fname,xreg,time,Nregs,Nt,X2d,"dridx")
        Allocate (Areal_mean_Dr(Ntime))
        Areal_mean_Dr(1:Ntime) = X3d(Nregs,1,Nmon+1:Nt)      ! Chi lay trung binh toan SEA
        Deallocate (X3d,xreg,xscal,time)

        Allocate (sst_dry(Nlon_sst,Nlon_sst),sst_wet(Nlon_sst,Nlon_sst))
        Allocate (sst_nor(Nlon_sst,Nlon_sst),sst_cli(Nlon_sst,Nlon_sst))
        Allocate (Pval(Nlon_sst,Nlon_sst))
        sst_dry=Miss; sst_wet=Miss; sst_nor=Miss; sst_cli=Miss; Pval=Miss
        ic=0
        Do ilon=1,Nlon_sst
           Do ilat=1,Nlat_sst
               if (Mask_sst(ilon,ilat)==1.AND..not.Empty_Ar1(sst_3d(ilon,ilat,1:Nt),Nt)) then
                  tmp1=0; tmp2=0; ntmp1=0; ntmp2=0; Sx=0; Sy=0
                  tmp3=0; tmp4=0; ntmp3=0; ntmp4=0; ic=ic+1
                  Do it=1,Ntime
                     ntmp4=ntmp4+1; tmp4=tmp4+sst_3d(ilon,ilat,it)    ! Climatology
                     if (Areal_mean_Dr(it)<=Dr_thres) then            ! Dry
                        ntmp1=ntmp1+1 
                        tmp1=tmp1+sst_3d(ilon,ilat,it)
                        Sx=Sx+sst_3d(ilon,ilat,it)*sst_3d(ilon,ilat,it)
                     else if (Areal_mean_Dr(it)>=W_thres) then        ! Wet
                        ntmp2=ntmp2+1
                        tmp2=tmp2+sst_3d(ilon,ilat,it)
                        Sy=Sy+sst_3d(ilon,ilat,it)*sst_3d(ilon,ilat,it)
                     else if (Areal_mean_Dr(it)>Dr_thres.AND.Areal_mean_Dr(it)<W_thres) then
                        ntmp3=ntmp3+1
                        tmp3=tmp3+sst_3d(ilon,ilat,it)
                     endif
                  Enddo
                  tmp1=tmp1/ntmp1; tmp2=tmp2/ntmp2
                  Sx=Sx/(ntmp1-1); Sy=Sy/(ntmp2-1)
                  Sx=SQRT(Sx-tmp1*tmp1); Sy=SQRT(Sy-tmp2*tmp2)
                  sst_dry(ilon,ilat)=tmp1
                  sst_wet(ilon,ilat)=tmp2
                  Pval(ilon,ilat)=P_value(tmp1, tmp2, Sx, Sy, ntmp1, ntmp2)
                  sst_nor(ilon,ilat)=tmp3/ntmp3
                  sst_cli(ilon,ilat)=tmp4/ntmp4
                  !print*,ntmp1,ntmp2,ntmp3,ntmp4
               endif
           Enddo
        Enddo

        ! Lon_sst(Nlon_sst), Lat_sst(Nlat_sst)
        
        Call wr_ncdf_2d (Fout,Lon_sst,Lat_sst, Nlon_sst,Nlat_sst,sst_dry, &
                     "sst_dry","Mean SST in case Dr_idx<-0.5", &
                     "C", "lon","lat", Miss)
        Call wr_ncdf_2d_apd (Fout,sst_wet, "sst_wet", Nlon_sst,Nlat_sst, &
                 "Mean SST in case Dr_idx>0.5", "C", "lon", "lat", Miss)
        Call wr_ncdf_2d_apd (Fout,Pval,"pval",Nlon_sst,Nlat_sst, &
                   "P-value of testing if sst_dry /= sst_wet","-","lon", "lat", Miss)
        Call wr_ncdf_2d_apd (Fout,sst_nor, "sst_nor", Nlon_sst,Nlat_sst, &
                 "Mean SST in -0.5<case Dr_idx<0.5", "C", "lon", "lat", Miss)
        Call wr_ncdf_2d_apd (Fout,sst_cli, "sst_cli", Nlon_sst,Nlat_sst, &
                 "Climatological Mean SST", "C", "lon", "lat", Miss)
      End Subroutine SST_Composite_maps
      !---------------------------------------------------------------------------------
      Real Function P_value(Xtb, Ytb, Sx, Sy, Nx, Ny)   ! Testing whether Xmean and Ymean are different
        Real :: t, Xtb, Ytb, Sx, Sy, Nx, Ny, df, Sx2, Sy2
        Real :: T_Dist, A, tmp1, tmp2, tmp3, tmp4, pval

        df=nx+ny-2
        if (Xtb/=Miss .AND. Ytb/=Miss .AND. Sx>0 .AND. Sy>0 .AND. Nx>1 .AND. Ny>1) then
           Sx2=Sx*Sx; Sy2=Sy*Sy     ! S^2
           tmp1=1.0/real(nx)+1.0/real(ny); tmp1=1.0/SQRT(tmp1) ! 1/Sqrt(1/n1 + 1/n2)
           tmp2=SQRT(df)
           tmp3=SQRT((nx-1)*Sx2 + (ny-1)*Sy2)
           tmp4=tmp3/tmp2
           A=tmp1/tmp4
           t=ABS(A*(Xtb-Ytb))
           pval=T_Dist (t, df)
        else
           pval=Miss
        endif
        P_value = pval
      End Function P_value
      !--------------------------------------------------------------------
      Subroutine Read_UV_Srf (Fu,Fv)
        Character (Len=*) :: Fu,Fv
        Real, Parameter :: Miss_uv=-9.96921e+36
        Integer, Parameter :: Yr_sta=1948
        Real, Allocatable :: XX(:,:,:),YY(:,:,:)
        Real, Allocatable :: X1d(:), Y1d(:), tlon(:), tlat(:), xtime(:)
        Integer :: ist, i, j, ilon, ilat, ipoint, nt, it, nlon, nlat
        Integer :: it1, it2, ilon1, ilon2,ilat1,ilat2, imon, M, Md2

        print*, "Reading SLP ..."
        Call rd_ncdf_3d (Fu,tlon,tlat,xtime,nlon,nlat,nt,XX,"uwnd") !trim(varn))
        Deallocate (tlon,tlat,xtime)
        Call rd_ncdf_3d (Fv,tlon,tlat,xtime,nlon,nlat,nt,YY,"vwnd")

        Where (XX==Miss_uv) XX=Miss; Where (YY==Miss_uv) YY=Miss;
        Where (IsNaN(XX)) XX=Miss;   Where (IsNaN(YY)) YY=Miss;

        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
        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

        Nlon_uvs=ilon2-ilon1+1; Nlat_uvs=ilat1-ilat2+1 
        Allocate (us_3d(Nlon_uvs,Nlat_uvs,Ntime),vs_3d(Nlon_uvs,Nlat_uvs,Ntime))
        Allocate (lon_uvs(Nlon_uvs),lat_uvs(Nlat_uvs))

        i=0
        Do ilon=1,Nlon
           if (ilon>=ilon1.AND.ilon<=ilon2) then
              i=i+1; Lon_uvs(i)=tlon(ilon)
              j=Nlat_uvs+1
              Do ilat=1,Nlat
                 if (ilat>=ilat2.AND.ilat<=ilat1) then
                    j=j-1; Lat_uvs(j)=tlat(ilat)
                    us_3d(i,j,1:Ntime) = XX(ilon,ilat,it1:it2)
                    vs_3d(i,j,1:Ntime) = YY(ilon,ilat,it1:it2)
                 endif
              Enddo
           endif
        Enddo

!      Real, allocatable :: lon_uvs(:), lat_uvs(:), Mask_uvs(:,:)
!      Integer :: Nlon_uvs, Nlat_uvs

      End Subroutine Read_UV_Srf
      !--------------------------------------------------------------------
      Subroutine Read_SLP (Fname)
        Character (Len=*) :: Fname
        Real, Parameter :: Miss_slp=-9.96921e+36
        Real, Parameter :: Miss_uv=-9.96921e+36
        Integer, Parameter :: Yr_sta=1948
        Real, Allocatable :: XX(:,:,:),YY(:,:,:)
        Real, Allocatable :: X1d(:), Y1d(:), tlon(:), tlat(:), xtime(:)
        Integer :: ist, i, j, ilon, ilat, ipoint, nt, it, nlon, nlat
        Integer :: it1, it2, ilon1, ilon2,ilat1,ilat2, imon, M, Md2

        print*, "Reading SLP ..."
        Call rd_ncdf_3d (Fname,tlon,tlat,xtime,nlon,nlat,nt,XX,"slp") !trim(varn))
        Where (XX==Miss_slp) XX=Miss
        Where (IsNaN(XX)); XX=Miss; Endwhere
        Where (XX==Miss_slp); 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
        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 ! + Nmon         ! Bo nam dau tien
        it2=it1+Ntime-1
        print*,"SST: Yr_sta,Yr1,it1,it2,Ntime: ", Yr_sta,Yr1,it1,it2,Ntime
        !print*,tlon(ilon1),tlon(ilon2),tlat(ilat1),tlat(ilat2), it1/Nmon, it2/Nmon

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

        Allocate (YY(Nlon_slp,Nlat_slp,Ntime),Lon_slp(Nlon_slp), Lat_slp(Nlat_slp))
        Allocate (slp_3d(Nlon_slp,Nlat_slp,Ntime))         ! Luu lai de su dung trong composite maps, Nt = do dai chuoi chua them Md2

        i=0
        Do ilon=1,Nlon
           if (ilon>=ilon1.AND.ilon<=ilon2) then
              i=i+1; Lon_slp(i)=tlon(ilon)
              j=Nlat_slp+1
              Do ilat=1,Nlat
                 if (ilat>=ilat2.AND.ilat<=ilat1) then
                    j=j-1; Lat_slp(j)=tlat(ilat)
                    YY(i,j,1:Ntime) = XX(ilon,ilat,it1:it2)
                 endif
              Enddo
           endif
        Enddo
        slp_3d(:,:,1:Ntime) = YY(:,:,1:Ntime) ! YY(:,:,Md2+1:Ntime-Md2)    ! Luu lai

      End Subroutine Read_SLP
      !--------------------------------------------------------------------
      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(:)
        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)
                 endif
              Enddo
           endif
        Enddo
        sst_3d(:,:,1:Nt) = YY(:,:,1:Nt) ! YY(:,:,Md2+1:Ntime-Md2)    ! Luu lai
        !print*,"check: Nt,Ntime=",Nt,Ntime,Md2+1,"-",Ntime-Md2; pause
        ! Processing for MCA:
        SeaPoint=0
        Do i=1,Nlon_sst
           Do j=1,Nlat_sst
              X1d(1:Ntime) = YY(i,j,1:Ntime)
              if (.not.Empty_Ar1(X1d,Ntime)) then
                 SeaPoint=SeaPoint+1; Mask_sst(i,j)=1
              endif
           Enddo
        Enddo
        ! Anomaly:
        print*,"SST: ",Nt,SeaPoint
        Allocate (sst(Nt,SeaPoint), Y1d(Nt))
        ipoint=0
        Do ilat=1, Nlat_sst
           Do ilon=1, Nlon_sst
              if (Mask_sst(ilon,ilat)/=Miss) then
                 X1d(1:Ntime)=YY(ilon,ilat,1:Ntime)
                 Call Anomaly (X1d,Ntime)
                 ipoint=ipoint+1
                 sst(1:Nt,ipoint)=X1d(1:Nt)
              endif
           Enddo
        Enddo
        print*,"SST: ",ipoint
        ! Moving Average:
        !M=Mov_Ave_Step               ! Buoc truot = 13 thang
        !Nt=Ntime-M+1 !--> Da bo sung
        !Allocate (sst(Nt,SeaPoint), Y1d(Nt))
        !ipoint=0
        !Do ilat=1, Nlat_sst
        !   Do ilon=1, Nlon_sst
        !      if (Mask_sst(ilon,ilat)/=Miss) then
        !         ipoint=ipoint+1
        !         X1d(1:Ntime)=YY(ilon,ilat,1:Ntime)
        !         Call Moving_Ave (X1d,Ntime,M,Y1d)
        !         sst(1:Nt,ipoint)=Y1d(1:Nt) ! YY(ilon,ilat,1:Ntime)
        !      endif
        !   Enddo
        !Enddo
        !Ntime = Nt
        Deallocate (YY, XX, X1d,Y1d,tlat,tlon,xtime)
      End Subroutine Read_sst_data
      !---------------------------------------------------------------------------------
      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))>1300) 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"