Use netCDF
      Integer, Parameter :: Nmon=12, Yr1=1981, Yr2=2014, Nyrs=(Yr2-Yr1)+1
      Real, Parameter :: Miss=-99.0, missing_value = 9.999e+20
      Real, Parameter :: Dlon=0.5, Dlat=0.5, Lon1=60.0, Lon2=290.0, Lat1=-20.0, lat2=60.0

      Integer :: Ntime,Nst,Nlon,Nlat, ilon,ilat,ist, Nlag      

      Real, Allocatable :: sst(:,:,:,:)   ! (Nlon,Nlat,Nyrs,Nmon)
      Real, Allocatable :: TC(:,:)        ! (Nyrs,Nmon)
      Real, Allocatable :: xlon(:),xlat(:),xyrs(:),xmon(:)  ! Nlon,Nlat,Nyrs,Nmon
      Character (Len=250) :: Fname, pca_out, eof_out
      Character (Len=10)  :: vname, vartc, varf
      !-------------------------------------------------------------------------------------------
      Real, Allocatable :: Corre(:,:,:,:), Pval(:,:,:,:) ! (Nlon,Nlat,Nmon,0:Nlag)
      !-------------------------------------------------------------------------------------------
      Fname="/home/tanpv/Data/TCs/VES_Monthly_TC_NTC_Count_jp_new.nc"
      vartc="tc"
      Call Read_TC (Fname,trim(vartc))
      vname="sst" !"vor" !"prmsl" ! "sst"
      !vname="vor"
      !vname="prmsl"
      !vname="vws"
      !vname="h850"
      Fname="/home/tanpv/Data/CFSR0.5/Vars/"//trim(vname)//".nc"

      !vname="vws"
      !Fname="/home/tanpv/Data/CFSR0.5/Vars/uv850_200.nc"

      varf=vname
      if (trim(varf)=="h850") then
         varf="hgt"
         Call Read_var_4d (Fname,varf)
      else
         Call Read_var (Fname,varf)
      endif
      !Call Read_uv (Fname,vname)

      !---------------------------------------------------------------------------------
      Nlag=13  ! 0-12
      Fname="Corr_"//trim(vartc)//"_"//trim(vname)//"_Reanalysis.nc"
      Call Calc_TQ (Fname,vname)

      CONTAINS
      !-------------------------------------------------------------------------------------------
      Subroutine Read_TC (Fname,vname)
        Character (Len=*) :: Fname, vname
        Real, Allocatable :: xyear(:),xmon(:), xtc(:,:)
        Integer :: Ny, Nm, Nt, it, ix, iy, imon, iyr

        Call rd_ncdf_2d (Fname,xyear,xmon,Ny,Nm,xtc,trim(vname))
        Allocate (TC(Nyrs,Nmon))
        iy=0
        Do iyr=Yr1,Yr2
           Do it=1,Ny
              if (int(xyear(it))==iyr) then
                 iy=iy+1; TC(iy,1:Nmon)=xtc(it,1:Nmon)
              endif
           Enddo
           !print"(I5,12F5.1)",iy,(TC(iy,1:Nmon))
        Enddo
      End Subroutine Read_TC
      !-------------------------------------------------------------------------------------------
      Subroutine Calc_TQ (Fname, vname)
        Character (Len=*) :: Fname, vname
        Integer :: Lag1, lag,N,imtc,imvar,ytc1,ytc2,yva1,yva2
        Integer :: iyr, imon,ilon,ilat
        Character (Len=50) :: var_nm, t_unit, var_ln, var_unit, &
                              lon_nm, lat_nm, z_nm, t_nm
        Real :: Corr, Prob
        Real, Allocatable :: X(:), Y(:), xmon(:), xlag(:)

        Allocate (X(Nyrs),Y(Nyrs))
        Allocate (Corre(1:Nlon,1:Nlat,1:Nmon,1:Nlag), Pval(1:Nlon,1:Nlat,1:Nmon,1:Nlag)) 
        Do imtc=1,Nmon  ! TC mon
           Do lag1=1,Nlag
              lag=lag1-1
              Call Calc_init_mon_year (lag,Nyrs,imtc,imvar,ytc1,ytc2,yva1,yva2)
              !print "(20I5)",imtc,imvar,ytc1,ytc2,yva1,yva2
              X(:)=Miss
              X(ytc1:ytc2)=TC(ytc1:ytc2,imtc)
              !print*,X
              Do ilat=1,Nlat
                 Do ilon=1,Nlon
                    Y(:)=Miss; Y(yva1:yva2)=sst(ilon,ilat,yva1:yva2,imvar)
                    !print "(35F7.1)", X; print "(35F7.1)", Y
                    Call Correlation (X, Y, Nyrs, Corr, Prob, Miss)
                    !print*,Corr, Prob; pause
                    Corre(ilon,ilat,imtc,lag1)=Corr
                    Pval (ilon,ilat,imtc,lag1)=Prob
                 Enddo
              Enddo
           Enddo
        Enddo
        Allocate (xmon(Nmon),xlag(Nlag))
        Do lag=1,Nlag; xlag(lag)=lag-1; Enddo
        Do imon=1,Nmon; xmon(imon)=imon; Enddo

        lon_nm="lon"; lat_nm="lat"; z_nm="mon"; t_nm="lag"; var_unit="-";
        t_unit="months since 1981-01-01 "
        var_nm="cor"; var_ln="Lag Correlation Coefficients - TC-"//trim(vname) 
        Call wr_ncdf_4d (Fname,xlon,xlat,xmon,xlag,nlon,nlat,nmon,nlag,Corre,var_nm, &
                       t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
        var_nm="pval"; var_ln="P-values of t-test"
        Call wr_ncdf_4d_apd (Fname,Pval, var_nm, nlon,nlat,nmon,nlag, &
                           var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
      End Subroutine Calc_TQ
      !-------------------------------------------------------------------------------------------
      Subroutine Calc_init_mon_year (Lag,N,imtc,imvar,ytc1,ytc2,yva1,yva2)
        ! Input: Lag, N, imtc
        Integer :: Lag,N,imtc,imvar,ytc1,ytc2,yva1,yva2

        if (Lag==0) then
           imvar=imtc; ytc1=1; ytc2=N; yva1=1; yva2=N
        else
           if (imtc>Lag) then
              imvar=imtc-Lag; ytc1=1; ytc2=N; yva1=1; yva2=N
           else
              imvar=Nmon+imtc-Lag; ytc1=2; ytc2=N; yva1=1; yva2=N-1
           endif
        endif
      End Subroutine Calc_init_mon_year
      !-------------------------------------------------------------------------------------------
      Subroutine Read_uv (Fname,vname)
        Character (Len=*) :: Fname,vname
        Integer :: Ntime, Nx, Ny, Nlev, it, ix, iy, imon, iyr
        Real, Allocatable :: xtmp(:,:,:), u(:,:,:,:), v(:,:,:,:)
        Real, Allocatable :: lon(:),lat(:), lev(:), xt(:)

        Call rd_ncdf_4d (Fname,lon,lat,lev,xt,nx,ny,nlev,ntime,u,"u")

        !Call rd_ncdf_3d (Fname,lon,lat,xt,nx,ny,ntime, X, trim(vname)) !"sst")

        Nlon=Int((Lon2-Lon1)/Dlon)+1; Nlat=Int((lat2-Lat1)/Dlat)+1
        ix=Int(Lon1/Dlon)+1; iy=Int((Lat1+90)/Dlat)+1

        Allocate (xtmp(Nlon,Nlat,Ntime))
        Allocate (xlon(Nlon),xlat(Nlat))
        xtmp(1:Nlon,1:Nlat,1:Ntime) = u(ix:ix+Nlon,iy:iy+Nlat,1,1:Ntime) &
                                     -u(ix:ix+Nlon,iy:iy+Nlat,2,1:Ntime)
        xlon(1:Nlon) = lon(ix:ix+Nlon); xlat(1:Nlat) = lat(iy:iy+Nlat)

        Allocate (sst(Nlon,Nlat,Nyrs,Nmon),xyrs(Nyrs),xmon(Nmon))
        Do iyr=Yr1,Yr2; xyrs(iyr-Yr1+1)=iyr; Enddo
        Do imon=1,Nmon; xmon(imon)=imon; Enddo

        it=0
        Do iyr=1,Nyrs
           Do imon=1,Nmon
              it=it+1
              sst(1:Nlon,1:Nlat,iyr,imon)=xtmp(1:Nlon,1:Nlat,it)
           Enddo
        Enddo
        Where (sst==missing_value) sst=Miss

        print*,nx,ny,ntime,Nlon,Nlat,ix,iy
        print*,xlon(1),xlon(Nlon),xlat(1),xlat(Nlat)
        print*,ix,ix+Nlon,iy,iy+Nlat
      End Subroutine Read_uv
      !-------------------------------------------------------------------------------------------
      Subroutine Read_var (Fname,vname)
        Character (Len=*) :: Fname,vname
        Integer :: Ntime, Nx, Ny, it, ix, iy, imon, iyr
        Real, Allocatable :: xtmp(:,:,:), X(:,:,:),lon(:),lat(:),xt(:)

        !if (trim(vname)=="h850") vname="hgt"
        Call rd_ncdf_3d (Fname,lon,lat,xt,nx,ny,ntime, X, trim(vname)) !"sst")
        Nlon=Int((Lon2-Lon1)/Dlon)+1; Nlat=Int((lat2-Lat1)/Dlat)+1
        ix=Int(Lon1/Dlon)+1; iy=Int((Lat1+90)/Dlat)+1
        Allocate (xtmp(Nlon,Nlat,Ntime))
        Allocate (xlon(Nlon),xlat(Nlat))
        xtmp(1:Nlon,1:Nlat,1:Ntime) = X(ix:ix+Nlon,iy:iy+Nlat,1:Ntime)
        xlon(1:Nlon) = lon(ix:ix+Nlon); xlat(1:Nlat) = lat(iy:iy+Nlat)
        
        Allocate (sst(Nlon,Nlat,Nyrs,Nmon),xyrs(Nyrs),xmon(Nmon))
        Do iyr=Yr1,Yr2; xyrs(iyr-Yr1+1)=iyr; Enddo
        Do imon=1,Nmon; xmon(imon)=imon; Enddo

        it=0
        Do iyr=1,Nyrs
           Do imon=1,Nmon
              it=it+1
              sst(1:Nlon,1:Nlat,iyr,imon)=xtmp(1:Nlon,1:Nlat,it)
           Enddo
        Enddo
        if (trim(vname)=="sst") then
           Where (sst/=missing_value) sst=sst-273.15
        endif
        Where (sst==missing_value) sst=Miss
        
        print*,nx,ny,ntime,Nlon,Nlat,ix,iy
        print*,xlon(1),xlon(Nlon),xlat(1),xlat(Nlat)
        print*,ix,ix+Nlon,iy,iy+Nlat
      End Subroutine Read_var
      !-------------------------------------------------------------------------------------------
      Subroutine Read_var_4d (Fname,vname)
        Character (Len=*) :: Fname,vname
        Integer :: Ntime, Nx, Ny, Nz, it, ix, iy, imon, iyr
        Real, Allocatable :: xtmp(:,:,:), X(:,:,:,:),lon(:),lat(:),lev(:),xt(:)

        !Call rd_ncdf_3d (Fname,lon,lat,xt,nx,ny,ntime, X, trim(vname)) !"sst")
        Call rd_ncdf_4d (Fname,lon,lat,lev,xt,nx,ny,nz,ntime,X,trim(vname))

        Nlon=Int((Lon2-Lon1)/Dlon)+1; Nlat=Int((lat2-Lat1)/Dlat)+1
        ix=Int(Lon1/Dlon)+1; iy=Int((Lat1+90)/Dlat)+1
        Allocate (xtmp(Nlon,Nlat,Ntime))
        Allocate (xlon(Nlon),xlat(Nlat))
        xtmp(1:Nlon,1:Nlat,1:Ntime) = X(ix:ix+Nlon,iy:iy+Nlat,1,1:Ntime)
        xlon(1:Nlon) = lon(ix:ix+Nlon); xlat(1:Nlat) = lat(iy:iy+Nlat)

        Allocate (sst(Nlon,Nlat,Nyrs,Nmon),xyrs(Nyrs),xmon(Nmon))
        Do iyr=Yr1,Yr2; xyrs(iyr-Yr1+1)=iyr; Enddo
        Do imon=1,Nmon; xmon(imon)=imon; Enddo

        it=0
        Do iyr=1,Nyrs
           Do imon=1,Nmon
              it=it+1
              sst(1:Nlon,1:Nlat,iyr,imon)=xtmp(1:Nlon,1:Nlat,it)
           Enddo
        Enddo
        if (trim(vname)=="sst") then
           Where (sst/=missing_value) sst=sst-273.15
        endif
        Where (sst==missing_value) sst=Miss

        print*,nx,ny,ntime,Nlon,Nlat,ix,iy
        print*,xlon(1),xlon(Nlon),xlat(1),xlat(Nlat)
        print*,ix,ix+Nlon,iy,iy+Nlat
      End Subroutine Read_var_4d

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