Use NetCDF

   Real, Parameter :: Miss=-99.0
   Real, Parameter :: Slon1=110, Slon2=289, Slat1=-30, Slat2=30
   Integer, Parameter :: Nmon=12 

   Character (Len=255) :: Fname
   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
   !----------------------------------------------------------------------------------------------
   ! For PCA:
      Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:)
      Real*8, Allocatable :: EOFs(:,:), PCs(:,:)
      Real*8, Allocatable   :: Eigenvalues(:), ECumulate(:)
      Integer :: Nt, Nx   ! Kich thuoc thuc cua mang X
      Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K
      Real*8 :: Factor

      Character (Len=50) :: pca_out, eof_out

      Namelist /Model_Info/ pca_out, eof_out, &
                            switch, standard, PC_Stard, Nmode
   !----------------------------------------------------------------------------------------------
   Yr1=1979; Yr2=2004

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

      Open (1,file="Namelist_PCA_New.txt",status="Old")
      Read (1,Model_Info)
      Close (1)

      Nt = Ntime; Nx = SeaPoint
      !---------------------------------------------------------------------------------
      if (switch==1) then
         Mcount=Nt; Ncount=Nx
      else if (switch==0) then
         Mcount=Nx; Ncount=Nt
      endif
      Factor = Sqrt(Real(Nt))

      Allocate (xdat(Nt,Nx), ydat(Nt,Nx))
      Allocate (Xtb(Nx), Sx(Nx))   ! TB, DLC, ma tran Covar hoac Corre
      Allocate (Eigenvalues(Mcount))          ! Tri rieng va Vecto rieng (1D)
      Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode))
      !---------------------------------------------------------------------------------
      Xdat(1:Nt,1:Nx) = sst2d(1:Nt,1:Nx)
      Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, &
                        Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode)
      Call Write_EOF_PCA (pca_out, eof_out)
   CONTAINS
   !----------------------------------------------------------------------------------------------
      Subroutine Write_EOF_PCA (F_pca, F_eof)
         Character (Len=*) :: F_eof, F_pca
         Character (Len=255) :: var_ln, t_unit
         Integer :: i,j,k,iidx,it
         Real :: tmp
         Real, Allocatable :: PCA(:,:,:), xmode(:), EOF1(:,:), X1(:), Xtmp(:), xtime(:)

         Allocate (PCA(Nlon,Nlat,Nmode), xmode(Nmode), Xtmp(Nmode)) !, EOF1(Nmode,Nt))
         Allocate (EOF1(Nt,Nmode))
         !-------------------------------------------------------------------------         
         ! Write to PCs file:
         PCA(:,:,:)=Miss; k=0
         Do j=1, Nlat
            Do i=1, Nlon
               if (Mask_sst(i,j)>0) then
                  k=k+1; PCA(i,j,1:Nmode) = PCs(k,1:Nmode)
               endif
            Enddo
         Enddo
         Do i=1, Nmode; xmode(i)=i; Enddo         ! PCA's modes
         var_ln="Spatial distribution of PCAs of SST for period 1961-2015"
         t_unit="months since 1961-01-01 "
         !if (iidx==1) then
            Call wr_ncdf_3d (F_pca,lon,lat,xmode,nlon,nlat,Nmode, PCA, &
              "pc", trim(t_unit),trim(var_ln), &
               "-", "lon", "lat","mode", Miss)
         !else
         !   Call wr_ncdf_3d_apd (F_pca,PCA,trim(IDX_nm(iidx)),nlon,nlat, &
         !        Nmode, trim(var_ln), "-", "lon", "lat", "mode", Miss)
         !endif
         !-------------------------------------------------------------------------
         Xtmp(1:Nmode) = Eigenvalues(1:Nmode)
         Call wr_ncdf_1d_apd (F_pca, Xtmp, &
                "eig", Nmode, &
                "Eigen values", "-", "mode", Miss)

         tmp=0.0
         Do i=1,Mcount
            tmp=tmp+Eigenvalues(i)
            if (I<=Nmode) then
               ECumulate(i)=0
               Do j=1,I
                  ECumulate(i)=ECumulate(i)+Eigenvalues(j)
               Enddo
            endif
         Enddo
         Do i=1,Nmode
            ECumulate(i)=ECumulate(i)/tmp; Xtmp(i)=Eigenvalues(i)/tmp
         Enddo

         Call wr_ncdf_1d_apd (F_pca, Xtmp, &
                "var",Nmode, &
                "Percent Contribution ", &
                "-", "mode", Miss )

         Xtmp(1:Nmode) = ECumulate(1:Nmode)
         Call wr_ncdf_1d_apd (F_pca, Xtmp, &
                "cum",Nmode, &
                "Cumuluted Eigenvalues ", &
                "-", "mode", Miss )
         !-------------------------------------------------------------------------
         ! Write to EOF file:
         if (Allocated(xtime)) Deallocate (xtime); Allocate (xtime(1:Nt))
         Do it=1, Nt; xtime(it)=it; Enddo     ! Months
         !Do i=1, Nmode; EOF1(i,1:Nt)=EOFs(1:Nt,i); Enddo
         Do i=1, Nmode; EOF1(1:Nt,i)=EOFs(1:Nt,i); Enddo
         var_ln="Timeseries of eofs "
         !if (iidx==1) then
            !Call wr_ncdf_2d (F_eof,xmode,xtime, Nmode,Nt, EOF1, &
            !         trim(IDX_nm(iidx)), &
            !         trim(var_ln), "-", "mode", "mon", Miss)
            Call wr_ncdf_2d (F_eof,xtime,xmode, Nt,Nmode, EOF1, &
                     "eof", &
                     trim(var_ln), "-", "mon","mode", Miss)
         !else
            !Call wr_ncdf_2d_apd (F_eof,EOF1, trim(IDX_nm(iidx)), Nmode,Nt, &
            !         trim(var_ln), "-", "mode", "mon", Miss)
         !   Call wr_ncdf_2d_apd (F_eof,EOF1, trim(IDX_nm(iidx)), Nt,Nmode, &
         !            trim(var_ln), "-", "mon", "mode", Miss)
         !endif
         print*,"Write_EOF_PCA: Done !!!"
      End Subroutine Write_EOF_PCA

   !----------------------------------------------------------------------------------------------
   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 > days.txt')
         Open (1,file="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(Ntime,SeaPoint))
      ix=0
      Do ilat=1,Nlat
         Do ilon=1,Nlon
            if (mask_sst(ilon,ilat)/=Miss) then
               ix=ix+1; sst2d(1:Ntime,ix)=sst3d(ilon,ilat,1:Ntime)
            endif
         Enddo
      Enddo
      ! Anomaly:
      Do ix=1,SeaPoint
         Call Anomaly (sst2d(1:Ntime,ix),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

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