Use NetCDF
      Real, Parameter :: Miss=-99.0
      Integer, Parameter :: Nmon=12, Nyrs=20, Ntime=Nyrs*Nmon, Nvar=7
      Real :: Year(Nyrs), To(Nyrs,Nmon), Po(Nyrs,Nmon), Tm(Nyrs,Nmon), Pm(Nyrs,Nmon)
      Real :: xtime(Ntime), Tobs(Ntime), Pobs(Ntime), Tmod(Ntime), Pmod(Ntime)
      Real :: PDoo(Nyrs,Nmon), PDom(Nyrs,Nmon), PDmo(Nyrs,Nmon)
      Real :: PDSI_oo(Ntime), PDSI_om(Ntime), PDSI_mo(Ntime)
      Real :: xvar(Nvar)
      Integer :: i,j,k
      ! 1: Tobs, 2: Pobs, 3: Tmod, 4: Pmod, 5: PDSI_oo, 6: PDSI_om, 7: PDSI_mo
      Real, Allocatable :: X(:,:), XX(:,:,:), xlon(:)

      Call Read_T_P
      Open (1, file="To_Po.txt",status="old")
      Call Read_Data (PDoo,PDSI_oo,Nyrs,Nmon)
      Close(1)
      Open (1, file="To_Pm.txt",status="old")
      Call Read_Data (PDom,PDSI_om,Nyrs,Nmon)
      Close(1)
      Open (1, file="Tm_Po.txt",status="old")
      Call Read_Data (PDmo,PDSI_mo,Nyrs,Nmon)
      Close(1)
      ! Convert inch to mm and F to C:
      Do i=1, Ntime
         Pobs(i)=Pobs(i)*25.4; Pmod(i)=Pmod(i)*25.4
         Tobs(i)=(Tobs(i)-32.)/1.8; Tmod(i)=(Tmod(i)-32.)/1.8
      Enddo
      !
      !Allocate (X(Ntime,Nvar))
      !X(:,1)=Tobs(:); X(:,2)=Pobs(:); X(:,3)=Tmod(:); X(:,4)=Pmod(:)
      !X(:,5)=PDSI_oo(:); X(:,6)=PDSI_om(:); X(:,7)=PDSI_mo(:)

      Allocate (XX(Ntime,1,Nvar),xlon(1))
      XX(:,1,1)=Tobs(:); XX(:,1,2)=Pobs(:); XX(:,1,3)=Tmod(:); XX(:,1,4)=Pmod(:)
      XX(:,1,5)=PDSI_oo(:); XX(:,1,6)=PDSI_om(:); XX(:,1,7)=PDSI_mo(:)

      Do i=1, Ntime; xtime(i)=i; Enddo
      Do i=1, Nvar; xvar(i)=i; Enddo
      xlon(1)=1
    
      !Call wr_ncdf_2d ("Sen_test.nc",xtime,xvar,Ntime,Nvar,X,"drou", &
      !   "1: Tobs, 2: Pobs, 3: Tmod, 4: Pmod, 5: PDSI_oo, 6: PDSI_om, 7: PDSI_mo", &
      !                  "-", "mon", "var", Miss)

      ! wr_ncdf_3d (filenm,lon,lat,t,nlon,nlat,nt,var,var_nm,t_unit, var_long_name, var_unit, lon_nm, lat_nm, t_nm, flag)
      Call wr_ncdf_3d ("Sen_test_3D.nc",xtime,xlon,xvar,Ntime,1,Nvar,XX, &
           "dr","months since 1980-01-01 ", &
           "1: Tobs, 2: Pobs, 3: Tmod, 4: Pmod, 5: PDSI_oo,"// &
           " 6: PDSI_om, 7: PDSI_mo", &
           "-", "mon", "dump", "var", Miss)

      CONTAINS
      Subroutine Read_T_P
        Integer :: i,j
        
        Open (1, file="T_P.txt", status="old")
        Read (1,*)
        Call Read_Data (To,Tobs,Nyrs,Nmon)
        Read (1,*)
        Call Read_Data (Po,Pobs,Nyrs,Nmon)
        Read (1,*)
        Call Read_Data (Tm,Tmod,Nyrs,Nmon)
        Read (1,*)
        Call Read_Data (Pm,Pmod,Nyrs,Nmon)
        Close(1)
      End Subroutine Read_T_P

      Subroutine Read_Data (X,Y,N,M)
        Integer :: i,j, it
        Real :: X(Nyrs,Nmon), Y(N*M)

        Do i=1,N
           Read (1,*) Year(i), (X(i,j),j=1,M)
        Enddo
        it=0
        Do i=1,N
           Do j=1,M
              it=it+1; Y(it)=X(i,j)
           Enddo
        Enddo
      End Subroutine Read_Data
      !
      Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
      END