Use NetCDF

   Real, Parameter :: Miss=-99.0
   Real, Parameter :: Slon1=110, Slon2=289, Slat1=-30, Slat2=30
   !Real, Parameter :: Lonvn1=100, Lonvn2=111, Latvn1=8, Latvn2=24
   Real, Parameter :: Lonvn1=89.50, Lonvn2=146.51, Latvn1=-14.81, Latvn2=26.95
   Integer, Parameter :: Nmon=12 

   Character (Len=255) :: Fname, Fu, Fv, Fab
   Character (Len=1) :: Ch
   Real, Allocatable :: sst3d(:,:,:), sst2d(:,:), mask_sst(:,:), Lon(:), Lat(:)
   !Integer :: time(:)
   Integer :: i, sys, ncid, varid, timeid, Nlon, Nlat, Ntime, Yr1, Yr2, SeaPoint
   Real :: missing_value
   !----------------------------------------------------------------------------------------------
   Real, Allocatable :: Pre3d(:,:,:), pre2d(:,:), mask_pre(:,:), Lon_pre(:), Lat_pre(:)
   Integer :: Nlon_pre, Nlat_pre, LandPoint
   !----------------------------------------------------------------------------------------------
   Real, Allocatable :: spe3d(:,:,:), spe2d(:,:), mask_spe(:,:), Lon_spe(:), Lat_spe(:)
   Integer :: Nlon_spe, Nlat_spe, LandSpe
   !----------------------------------------------------------------------------------------------
   ! MCA:
         Real, Allocatable :: X(:,:), Y(:,:)
         Real, Allocatable :: C(:,:),  U(:,:), V(:,:), WW(:,:), W(:)
         Real, Allocatable :: Uk(:,:), Vk(:,:), Wk(:), Percent(:) !, WWk(:,:)
         Real, Allocatable :: Ak(:,:), Bk(:,:)
         Character (Len=3) :: VarX, VarY
         Integer :: Nmode !, i,j,k
         Logical :: Switch   ! =.true. if switch Prec and SST
         !Real :: temp1
         !Real, Allocatable :: temp2(:),temp3(:)
   !----------------------------------------------------------------------------------------------
      Yr1=1961; Yr2=2015

      Fname = "/work/projects/G09/SEACLID/SST/HadISST_sst.nc"
      Call Read_sst (Fname)

      !Fname = "/work/projects/G09/SEACLID/APH_2021/APH_025_SEA_Pre_V1101_T_V1808_1961-2015.nc"
      !Call Read_Prec (Fname)

      Fname = "/work/projects/G09/SEACLID/GPCC_CRU_SEA/SPEI_SPI_GPCC_CRU.nc"
      Call Read_SPEI (Fname)

      Nt=Ntime; Nmode=5

      !if (LandPoint<=SeaPoint) then
      !   Nx=LandPoint; Ny=SeaPoint; Switch=.False. !VarX="pre"; VarY="sst"
      !else
      !   Nx=SeaPoint; Ny=LandPoint; Switch=.True.
      !endif

      if (LandPoint<=LandSpe) then
         Nx=LandSpe; Ny=SeaPoint; Switch=.False.
      else
         Nx=SeaPoint; Ny=LandSpe; Switch=.True.
      endif
      print*,"Land point=",LandSpe !LandPoint

      !print*,"Land point=",LandPoint
      print*,"Sea  point=",SeaPoint

      Allocate (X(Nx,Nt),Y(Ny,Nt))
      Allocate (C(Nx,Ny),W(Ny)) !,U(Nx,Ny),V(Ny,Ny),W(Ny),WW(Ny,Ny))
      Allocate (Uk(Nx,Nmode),Vk(Ny,Nmode),Wk(Nmode),Percent(Nmode))
      Allocate (Ak(Nt,Nmode),Bk(Nt,Nmode))

      if (Switch) then
         !X(:,:)=sst2d(:,:); Y(:,:)=pre2d(:,:)
         X(:,:)=sst2d(:,:); Y(:,:)=spe2d(:,:)
      else
         !X(:,:)=pre2d(:,:); Y(:,:)=sst2d(:,:)
         X(:,:)=spe2d(:,:); Y(:,:)=sst2d(:,:)
      endif
      !Call MCA_SVD (X,Y,C,Nx,Ny,Nt,Uk,Vk,Wk,Ak,Bk,Percent,Nmode,Miss)

      Call Covar_Mtrx(X,Y,C,Nt,Nx,Ny,Miss)
      Call Write_C_Mtrx ("C_Matrix.nc")

      sys = system("sleep 1")
      sys = system("./run_mlab.sh")

      Call Calc_UVW_AB ("U.nc","V.nc","S.nc")
      print*,"Calc_UVW_AB ... Done"
      Deallocate (U,V,W,WW)

      Fu="/home/tanpv/Progs/2021/NC/U_Dr_SST.nc"; 
      Fv="/home/tanpv/Progs/2021/NC/V_Dr_SST.nc"; 
      Fab="/home/tanpv/Progs/2021/NC/AB_Dr_SST.nc"
      Call Write_MCA (Fu,Fv,Fab)
   CONTAINS
   !----------------------------------------------------------------------------------------------
   Subroutine Write_C_Mtrx (Fc)
      Character (Len=*) :: Fc
      Integer :: ix, iy
      Real :: x1d(Nx), y1d(Ny)

      Do ix=1,Nx; x1d(ix)=ix; Enddo
      Do iy=1,Ny; y1d(iy)=iy; Enddo
      Call wr_ncdf_2d (Fc,x1d,y1d,Nx,Ny,C,"cor","Covariance Matrix", "-","x","y", Miss)
   End Subroutine Write_C_Mtrx
   !----------------------------------------------------------------------------------------------
   Subroutine Calc_UVW_AB (Fu,Fv,Fw)
      Character (Len=*) :: Fu,Fv,Fw
      Real, Allocatable :: xdim(:), ydim(:)
      Integer :: N1,N2, i
      Real :: TotalVar, tmp

      print*,"Reading U,V,W files:"
      !Fname="MCA_U.nc"
      Call rd_ncdf_2d (Fu,xdim,ydim,N1,N2,U,"U")
      Do i=1,Nmode; Uk(1:Nx,i)=U(1:Nx,i); Enddo
      Deallocate (xdim, ydim, U)
      print*,"U: xdim, ydim=",N1,N2; pause
      !Fname="MCA_V.nc"
      Call rd_ncdf_2d (Fv,xdim,ydim,N1,N2,V,"V")
      Do i=1,Nmode; Vk(1:Nx,i)=V(1:Ny,i); Enddo
      Deallocate (xdim, ydim, V)
      print*,"V: xdim, ydim=",N1,N2
      !Fname="MCA_S.nc"
      Call rd_ncdf_2d (Fw,xdim,ydim,N1,N2,WW,"S")
      Deallocate (xdim, ydim)
      print*,"WW: xdim, ydim=",N1,N2
      if (Allocated(W)) Deallocate (W)
      Allocate (W(Ny))
      Do i=1,Ny; W(i)=WW(i,i); Enddo

      TotalVar = Sum(W(1:Nx))
      Do i=1,Nmode
         !Uk(1:Nx,i)=U(1:Nx,i)
         !Vk(1:Nx,i)=V(1:Ny,i)
         Wk(i)     =W(i)
         tmp=Sum(W(1:i))
         Percent(i)=tmp/TotalVar
      Enddo
      Ak = MatMul(Transpose(X),Uk) 
      Bk = MatMul(Transpose(Y),Vk)
   End Subroutine Calc_UVW_AB
   !----------------------------------------------------------------------------------------------
   Subroutine Write_MCA (Fu,Fv,Fab)
      Character (Len=*) :: Fu,Fv,Fab
      Real, Allocatable :: U(:,:,:),V(:,:,:),xmode(:), lon_land(:), lat_land(:)
      Real, Allocatable :: A(:,:), B(:,:), xtime(:), mask_land(:,:)
      Integer :: i, ilon,ilat,ix, Nlon_land, Nlat_land
      Character (Len=5) :: var_land
      
      ! For Prec:
      !Nlon_land=Nlon_pre; Nlat_land=Nlat_pre
      !Allocate (lon_land(Nlon_land),lat_land(Nlat_land),mask_land(Nlon_land,Nlat_land))
      !lon_land=lon_pre; lat_land=lat_pre; mask_land=mask_pre
      !var_land="prec"

      ! For SPEI:
      Nlon_land=Nlon_spe; Nlat_land=Nlat_spe
      Allocate (lon_land(Nlon_land),lat_land(Nlat_land),mask_land(Nlon_land,Nlat_land))
      lon_land=lon_spe; lat_land=lat_spe; mask_land=mask_spe
      var_land="spei"
      !--------------------------------------------------------------

      Allocate (xmode(Nmode),A(Nmode,Nt),B(Nmode,Nt),xtime(Nt))
      Do ix=1,Nmode; xmode(ix)=ix; Enddo
      Do ix=1,Nt;    xtime(ix)=ix; Enddo

      if (Switch) then  ! Left SST, Right Prec
         Allocate (U(Nlon,Nlat,Nmode),V(Nlon_pre,Nlat_pre,Nmode))
         ix=0; U=Miss; V=Miss
         Do ilat=1,Nlat
            Do ilon=1,Nlon
               if (mask_sst(ilon,ilat)/=Miss) then
                  ix=ix+1; 
                  Do i=1,Nmode
                     U(ilon,ilat,i)=Uk(ix,i)
                  Enddo
               endif
            Enddo
         Enddo
         ix=0
         Do ilat=1,Nlat_land
            Do ilon=1,Nlon_land
               if (mask_land(ilon,ilat)/=Miss) then
                  ix=ix+1; V(ilon,ilat,1:Nmode)=Vk(ix,1:Nmode)
               Endif
            Enddo
         Enddo
         Do i=1,Nmode
            A(i,1:Nt)=Ak(1:Nt,i); B(i,1:Nt)=Bk(1:Nt,i)
         Enddo
         Call wr_ncdf_3d (Fu,lon,lat,xmode,nlon,nlat,Nmode, U,"sst", "months since 1961-1-1 ","SST MCAs","-", "lon", "lat","mode", Miss)
         Call wr_ncdf_3d (Fv,lon_land,lat_land,xmode,nlon_land,nlat_land,Nmode, V,trim(var_land), "months since 1961-1-1 ",trim(var_land)//" MCAs","-", "lon", "lat","mode", Miss)
         Call wr_ncdf_2d (Fab,xmode,xtime,Nmode,Nt,A,"sst","MCA timeseries of SST", "-","mode","mon", Miss)
         Call wr_ncdf_2d_apd (Fab,B,trim(var_land),Nmode,Nt, "MCA timeseries of "//trim(var_land), "-", "mode", "mon", Miss)
         Call wr_ncdf_1d_apd (Fab, Percent,"cum", Nmode,"Accumulated percentage of variance", "-", "mode", Miss)
      else              ! Left Prec, Right SST
         Allocate (V(Nlon,Nlat,Nmode),U(Nlon_pre,Nlat_pre,Nmode))
         ix=0; U=Miss; V=Miss
         Do ilat=1,Nlat
            Do ilon=1,Nlon
               if (mask_sst(ilon,ilat)/=Miss) then
                  ix=ix+1; V(ilon,ilat,1:Nmode)=Vk(ix,1:Nmode)
               endif
            Enddo
         Enddo
         ix=0
         Do ilat=1,Nlat_land
            Do ilon=1,Nlon_land
               if (mask_land(ilon,ilat)/=Miss) then
                  ix=ix+1; U(ilon,ilat,1:Nmode)=Uk(ix,1:Nmode)
               Endif
            Enddo
         Enddo
         Do i=1,Nmode
            A(i,1:Nt)=Ak(1:Nt,i); B(i,1:Nt)=Bk(1:Nt,i)
         Enddo
         Call wr_ncdf_3d (Fv,lon,lat,xmode,nlon,nlat,Nmode, V,"sst", "months since 1961-1-1 ","SST MCAs","-", "lon", "lat","mode", Miss)
         Call wr_ncdf_3d (Fu,lon_land,lat_land,xmode,nlon_land,nlat_land,Nmode, U,trim(var_land), "months since 1961-1-1 ",trim(var_land)//" MCAs","-", "lon", "lat","mode", Miss)
         Call wr_ncdf_2d (Fab,xmode,xtime,Nmode,Nt,A,trim(var_land),"MCA timeseries of "//trim(var_land), "-","mode","mon", Miss)
         Call wr_ncdf_2d_apd (Fab,B,"sst",Nmode,Nt, "MCA timeseries of SST", "-", "mode", "mon", Miss)
         Call wr_ncdf_1d_apd (Fab, Percent,"cum", Nmode,"Accumulated percentage of variance", "-", "mode", Miss)
      endif
   End Subroutine Write_MCA
   !----------------------------------------------------------------------------------------------
   Subroutine Read_Prec (Fname)
      Character (Len=*) :: Fname
      Integer :: Ntime, ix,iy,nx, ny, nt, ilon, ilat, iday, year1, year2
      Integer :: it, it1, it2, imon, ilon1, ilon2, ilat1, ilat2
      Real, Allocatable :: y3d(:,:,:), x3d(:,:,:), x2d(:,:), x1d(:), y1d(:), xtime(:)
      Real, Allocatable :: xlon(:), xlat(:), tmp(:)


!   Real, Allocatable :: Pre3d(:,:,:), pre2d(:,:), mask_pre(:,:), Lon_pre(:), Lat_pre(:)
!   Integer :: Nlon_pre, Nlat_pre

      Call rd_ncdf_3d (Fname,x1d,y1d,xtime,nx,ny,Nt, x3d, "pr")
      year1=1961; year2=2015
      it=0
      Do iy=year1,year2
         Do imon=1,Nmon
            it=it+1
            if (iy==Yr1.AND.imon==1)    it1=it
            if (iy==Yr2.AND.imon==Nmon) it2=it
         Enddo
      Enddo
      Ntime=it2-it1+1
      print*,"Nt,it1,it2,Nx,Ny=",Nt,it1,it2,Nx,Ny
     ! Determine the VN domain:
      Do ilon=1,Nx-1
         if (x1d(ilon)<=Lonvn1.AND.x1d(ilon+1)>=Lonvn1) ilon1=ilon
         if (x1d(ilon)<=Lonvn2.AND.x1d(ilon+1)>=Lonvn2) ilon2=ilon+1
         !print*,ilon,x1d(ilon),ilon1,ilon2
      Enddo
      Do ilat=1,Ny-1
         if (y1d(ilat)<=Latvn1.AND.y1d(ilat+1)>=Latvn1) ilat1=ilat
         if (y1d(ilat)<=Latvn2.AND.y1d(ilat+1)>=Latvn2) ilat2=ilat+1
      Enddo

      ! Extract to domain: Mask sea points
      Nlon_pre=ilon2-ilon1+1; Nlat_pre=ilat2-ilat1+1
      print*,"ilon1, ilon2, ilat1, ilat2=",ilon1, ilon2, ilat1, ilat2
      print*,"Nlon_pre,Nlat_pre,Ntime=",Nlon_pre,Nlat_pre,Ntime

      Allocate (Pre3d(Nlon_pre,Nlat_pre,Ntime),mask_pre(Nlon_pre,Nlat_pre))
      Allocate (Lon_pre(Nlon_pre),Lat_pre(Nlat_pre))
      Pre3d=Miss; mask_pre=Miss; 
      Lon_pre(1:Nlon_pre)=x1d(ilon1:ilon2); Lat_pre(1:Nlat_pre)=y1d(ilat1:ilat2)
      !print*,Lon_pre(1:Nlon_pre),Lat_pre(1:Nlat_pre)

      Deallocate (x1d,y1d)
      Allocate (x1d(Ntime))

      ix=0; LandPoint=0
      Do ilon=ilon1,ilon2 !1,Nlon_pre
         ix=ix+1; iy=0 
         Do ilat=ilat1,ilat2 !1,Nlat_pre
            iy=iy+1
            x1d(1:Ntime)=x3d(ilon,ilat,it1:it2)
            if (.not.Empty_Ar1(x1d, Ntime, Miss)) then
               LandPoint=LandPoint+1
               mask_pre(ix,iy)=1 ! (ilon,ilat)=1
               Pre3d(ix,iy,1:Ntime)=x3d(ilon,ilat,it1:it2)
            endif
         Enddo
      Enddo
      !print*,"ix,iy,LandPoint=",ix,iy,LandPoint
      Allocate (pre2d(LandPoint,Ntime))
      ix=0
      Do ilat=1,Nlat_pre
         Do ilon=1,Nlon_pre 
            if (mask_pre(ilon,ilat)==1) then
               ix=ix+1
               x1d(1:Ntime)=Pre3d(ilon,ilat,1:Ntime)
               Call Anomaly (x1d,Ntime)
               pre2d(ix,1:Ntime)=x1d(1:Ntime)
            endif
         Enddo
      Enddo
      !print*,"LandPoint=",LandPoint
      !pause
   End Subroutine Read_Prec
   !----------------------------------------------------------------------------------------------
   Subroutine Read_SPEI (Fname)
      Character (Len=*) :: Fname
      Integer :: ilon,ilat,it,it1, Nt, Nt1, ix
      Real, Allocatable :: 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")
      !Deallocate (lon,lat,xtime)
      !Call rd_ncdf_2d (Fname,lon,lat,nlon,nlat,Mask,"land")
      Where (IsNaN(X3d)) X3d=Miss
      it1=Nmon+1; Nt1=Nt-Nmon

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

      Allocate (spe3d(nlon_spe,nlat_spe,Nt1),mask_spe(nlon_spe,nlat_spe))
      Allocate (x1d(Nt1))
      spe=Miss; mask_spe=Miss; LandSpe=0
      Do ilon=1,nlon_spe
         Do ilat=1,nlat_spe
            x1d(1:Nt1)=X3d(ilon,ilat,it1:Nt)
            if (.not.Empty_Ar1 (X1d,Nt1,Miss)) then
               mask_spe(ilon,ilat)=1
               LandSpe=LandSpe+1
               spe3d(ilon,ilat,1:Nt1)=X3d(ilon,ilat,it1:Nt)
            endif
         Enddo
      Enddo
      ! Reshape to 2D arr: Ntime=Nt1
      Allocate (spe2d(LandSpe,Ntime))
      ix=0
      Do ilat=1,Nlat_spe
         Do ilon=1,Nlon_spe
            if (mask_spe(ilon,ilat)/=Miss) then
               ix=ix+1; spe2d(ix,1:Ntime)=spe3d(ilon,ilat,1:Ntime)
            endif
         Enddo
      Enddo
   End Subroutine Read_SPEI
   !----------------------------------------------------------------------------------------------
   Subroutine Read_sst (Fname)
      Character (Len=*) :: Fname
      Integer :: ix,iy,nx, ny, nt, ilon, ilat, iday
      Integer :: it, it1, it2, ilon1, ilon2, ilat1, ilat2
      Real, Allocatable :: y3d(:,:,:), x3d(:,:,:), x2d(:,:), x1d(:), y1d(:), xtime(:)
      Real, Allocatable :: xlon(:), xlat(:), tmp(:)
      Logical :: status
      Integer :: ym1, ym2, xtype, lenst, attnum, date_fcs
      Character (Len=100) :: t_unit
      Character (Len=13) :: date0, timeres, idate0, idate1
      !-------------------------------------------------------------------------------------------
      status = nf90_open(fname, nf90_NoWrite, ncid)
      status = nf90_inq_varid(ncid,"sst", varid)
      status = nf90_get_att(ncid, varid, "missing_value", missing_value)

      status = nf90_inq_varid(ncid,"time", varid)
      status = nf90_inquire_attribute(ncid, varid, "time", xtype, lenst, attnum)
      status = nf90_get_att(ncid, varid, "units", t_unit)
      status =  nf90_close(ncid)

      LenSt=Len_trim(t_unit)
      Do i=1,LenSt
         Ch=t_unit(i:i)
         if (Ch=="1".OR.Ch=="2") exit
      Enddo
      date0=t_unit(i:i+8)
      timeres=" "//t_unit(1:4)

      Call rd_ncdf_3d (Fname,x1d,y1d,xtime,nx,ny,Nt, x3d, "sst")
      Where (x3d==missing_value) x3d=Miss

      ! Search it1, it2:
      ym1=Yr1*100+1; ym2=Yr2*100+12
      Do it=1,Nt
         iday=Int(xtime(it))
         sys = system ('rm -rf days.txt')
         write (idate0,*) iday
         sys = system ('date -u -d "'//trim(date0)//' '//trim(idate0)//trim(timeres)//' " +%Y%m | cat > /tmp/days.txt')
         Open (1,file="/tmp/days.txt"); Read (1,*) date_fcs; Close(1)
         if (date_fcs==ym1) it1=it
         if (date_fcs==ym2) it2=it  
         !print*,it,trim(idate0)," ",trim(timeres), date_fcs
      Enddo
      Ntime=it2-it1+1
      print*, ym1, ym2, it1, it2

      ! Search for lat, lon:
      Allocate (y3d(nx,ny,nt), xlon(nx), xlat(ny))
      Do ilat=1,Ny
         xlat(ilat)=y1d(Ny-ilat+1)
         y3d(:,ilat,:)=x3d(:,Ny-ilat+1,:)
      Enddo
      Do ilon=1,Nx
         if (ilon<=180) then
            xlon(ilon)=x1d(ilon+180)
            x3d(ilon,:,:)=y3d(ilon+180,:,:)
         else
            xlon(ilon)=x1d(ilon-180)+360
            x3d(ilon,:,:)=y3d(ilon-180,:,:)
         endif
      Enddo
      !Call wr_ncdf_3d ("sst.nc",xlon,xlat,xtime,nx,ny,nt,x3d,"sst","days since 1870-01-15","sst","C","lon","lat","mon", Miss)
      Do ilon=1,Nx-1
         if (xlon(ilon)>=Slon1.AND.xlon(ilon)<=Slon1+1) ilon1=ilon
         if (xlon(ilon)>=Slon2.AND.xlon(ilon)<=Slon2+1) ilon2=ilon 
      Enddo
      Do ilat=1,Ny-1
         if (xlat(ilat)>=Slat1.AND.xlat(ilat)<=Slat1+1) ilat1=ilat
         if (xlat(ilat)>=Slat2.AND.xlat(ilat)<=Slat2+1) ilat2=ilat
      Enddo

      ! Extract to domain: Mask sea points
      Nlon=ilon2-ilon1+1; Nlat=ilat2-ilat1+1
      Allocate (sst3d(Nlon,Nlat,Ntime),mask_sst(Nlon,Nlat), Lon(Nlon),Lat(Nlat), tmp(Ntime))

      ix=0; SeaPoint=0; mask_sst=Miss
      Do ilon=ilon1,ilon2
         ix=ix+1; Lon(ix)=xlon(ilon)
         iy=0
         Do ilat=ilat1,ilat2
            iy=iy+1; Lat(iy)=xlat(ilat)
            sst3d(ix,iy,1:Ntime)=x3d(ilon,ilat,it1:it2)
            tmp(1:Ntime)=sst3d(ix,iy,1:Ntime)
            if (.not.Empty_Ar1(tmp, Ntime, Miss)) then
               mask_sst(ix,iy)=1
               SeaPoint=SeaPoint+1
            endif
         Enddo
      Enddo
      ! Reshape to 2D arr:
      Allocate (sst2d(SeaPoint,Ntime))
      ix=0
      Do ilat=1,Nlat
         Do ilon=1,Nlon
            if (mask_sst(ilon,ilat)/=Miss) then
               ix=ix+1; sst2d(ix,1:Ntime)=sst3d(ilon,ilat,1:Ntime)
            endif
         Enddo
      Enddo
      ! Anomaly:
      Do ix=1,SeaPoint
         Call Anomaly (sst2d(ix,1:Ntime),Ntime)
      Enddo
      !print*,ilon1, xlon(ilon1),ilon2, xlon(ilon2)
      !print*,ilat1, xlat(ilat1),ilat2, xlat(ilat2)
   End Subroutine Read_sst
   !------------------------------------------------------------------------------
   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
   !------------------------------------------------------------------------------
   Logical Function Empty_Ar1 (X, N, Miss)
   ! Tra ve .True. neu X(N) khuyet toan bo SL
     Integer :: N, i, k
     Real :: X(N), Miss
     Logical :: Empty

     Empty=.False. 
     Do i=1, N
        if (X(i) == Miss) then; Empty=.True.; exit; endif
     Enddo
     Empty_Ar1=Empty 
   End Function Empty_Ar1
      !---------------------------------------------------------------------------------
      SUBROUTINE Covar_Mtrx(X,Y,C,Nt,Nx,Ny,Miss)
      Integer :: Nt, Nx, Ny, N1, L, J, K, IO, I
      Real :: X(Nx,Nt), Y(Ny,Nt), C(Nx,Ny), Xp(Nx,Nt), Yp(Ny,Nt)
      Real :: XBAR(Nx), Stdx(Nx), YBar(Ny), Stdy(Ny)
      Real :: Sx, Xtb, Sy, Ytb, Miss
      !---------------------------------------------------
      ! Anomaly:
      Do j=1, Nx
         Sx=0; Xtb=0; N1=0
         Do i=1, Nt
            if (X(j,i)/=Miss) then
               n1=n1+1; Xtb=Xtb+X(j,i); Sx=Sx+X(j,i)*X(j,i)
            Endif
         Enddo
         if (n1>nt/2) then
            XBAR(j)=Xtb/real(n1); STDx(j)=Sqrt(Sx/real(n1)-XBAR(j)*XBAR(j))
         endif
         Do i=1, Nt
            if (X(j,i)==Miss) X(j,i)=XBAR(j) ! Fill-up Miss by mean values
            Xp(j,i) = (X(j,i)-XBAR(j))
            !if (standard==1) Y(i,j) = Y(i,j)/STD(j)
         Enddo
      Enddo
      !
      Do j=1, Ny
         Sx=0; Xtb=0; N1=0
         Do i=1, Nt
            if (Y(j,i)/=Miss) then
               n1=n1+1; Xtb=Xtb+Y(j,i); Sx=Sx+Y(j,i)*Y(j,i)
            Endif
         Enddo
         if (n1>nt/2) then
            Ybar(j)=Xtb/real(n1); STDy(j)=Sqrt(Sx/real(n1)-Ybar(j)*Ybar(j))
         endif
         Do i=1, Nt
            if (Y(j,i)==Miss) Y(j,i)=Ybar(j) ! Fill-up Miss by mean values
            Yp(j,i) = (Y(j,i)-Ybar(j))
            !if (standard==1) Y(i,j) = Y(i,j)/STD(j)
         Enddo
      Enddo
      C = MatMul(Xp,Transpose(Yp))/Real(Nt)
      End SUBROUTINE Covar_Mtrx
   !----------------------------------------------------------------------------------------------
   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"
   !Include "/home/tanpv/Progs/SST/Sub_MCA.f90"
   END