Use NetCDF

   Real, Parameter :: Miss=-99.0
   Integer, Parameter :: IMiss=-99, Nmem=11, NLead=46
   Character (Len=2), Parameter :: xmem(Nmem)=(/"01","02","03","04","05","06","07","08","09","10","11"/)
   Integer, Allocatable :: time(:)

   Integer :: i, LenSt, sys
   Character (Len=255) :: Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2
   Character (Len=255) :: DirInEC, DirInTRMM, DirInStn 
   Character (Len=12) :: date0, timeres, idate0, idate1
   Character (Len=1) :: Ch
   Integer :: yy1,mm1,dd1,hh1,yy2,mm2,dd2,hh2
   Real*8 :: factor, offset, Missing
   Real, Allocatable :: lon(:), lat(:)
   Integer, Allocatable :: lead(:), fcsdate(:)
   Real, Allocatable :: EC(:,:,:,:,:)    ! (Nlon,Nlat,NLead,Nfcs,Nmem)
   Real, Allocatable :: TRMM(:,:,:), xlon(:), xlat(:), xtime(:)
   Integer, Allocatable :: Init_EC(:), Time_OBS(:)
   Real, Allocatable :: X4d(:,:,:,:), X3d(:,:,:) 
   Integer, Allocatable, Dimension (:) :: date_m1   ! (Nfcs)
   Integer :: Nlon, Nlat, Nfcs, date_fcs, Nlon0, Nlat0    ! Nfcs = So lan du bao, date_fcs = yyyymmdd cua ngay DB
   Integer :: Nlo, Nla, Ntime    ! Nday from 2000-01-01 to 2019-12-31

   Integer :: imem, ifcs, iy, im, id, NLd

   DirInEC  ="/home/tanpv/Progs/S2S/Data_EC/Daily/"
   DirInTRMM="/home/tanpv/Progs/S2S/Data_TRMM/"
   DirInStn ="/home/tanpv/Progs/S2S/Data_Stn/"

   ! Read ECMWF S2S data: 11 members, from 2000-01-02 to 2019-12-31
   ! Call Read_EC
   ! Read TRMM data:
   Fname=trim(DirInTRMM)//"TRMM_VN_2000-2019_Daily.nc"
   Call Read_TRMM (FName)
   Do imem=1,Nmem
      Fname=trim(DirInEC)//"Prec_daily_s2s_ecmwf_"//xmem(imem)//".nc"
      Call Read_EC (Fname)
      if (imem==1) Allocate (EC(Nlon,Nlat,NLead,Nfcs,Nmem))
      EC(:,:,:,:,imem) = X4d(:,:,:,:)
      if (imem<Nmem) Deallocate (lon,lat,lead,fcsdate,X4d)
   Enddo
   CONTAINS
   !----------------------------------------------------------------------------------------------
   Subroutine Read_EC (Fname)
      Character (Len=*) :: Fname
      Integer :: imem, ifcs

!      Do imem=1,Nmem !imem=1,Nmem
!         Fname=trim(DirInEC)//"Prec_daily_s2s_ecmwf_"//xmem(imem)//".nc" 
         Call rd_ncdf_4d_new (Fname,lon,lat,lead,fcsdate,nlon,nlat,NLd,Nfcs,X4d,"pr")
         ! lead(:): 1-46 days; fcsdate(:): in format of yyyymmdd
         print*,imem,Nlon,Nlat,NLead,Nfcs,Nmem 
!         if (imem==1) Allocate (EC(Nlon,Nlat,NLead,Nfcs,Nmem))
!         print*,"===="
!         EC(:,:,:,:,imem) = X4d(:,:,:,:)
!         if (imem<Nmem) Deallocate (lon,lat,lead,fcsdate,X4d)
!      Enddo
      !Allocate (Init_EC(Nfcs))
      !Do ifcs=1,Nfcs
      !   Init_EC(ifcs)=Int(fcsdate(ifcs))   ! Convert to Integer
      !Enddo
!      if (Allocated(X4d)) Deallocate (X4d)
   End Subroutine Read_EC
   !----------------------------------------------------------------------------------------------
   Subroutine Read_TRMM (FName)
      Character (Len=*) :: Fname
      Integer :: id

      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,Ntime,TRMM,"precipitation")
      sys = system ('rm -rf days.txt')
      Do id=1,Ntime
         write (Str,"(A,I5,A)") "2000-1-1 + ", Int(xtime(id))," days"
         sys = system ('date -u -d "'//trim(Str)//' " +%Y%m%d | cat >> days.txt')
      Enddo
      Allocate (Time_OBS(Ntime))
      Open (1,file="days.txt",Status="old")
      Do id=1,Ntime; Read (1,*) Time_OBS(id); Enddo
   End Subroutine Read_TRMM
   !----------------------------------------------------------------------------------------------
   Subroutine GetTime (fname, nt, date_fcs, scale_factor, add_offset,missing_value)
      Character (Len=*) :: fname
      Integer, Intent (Out) :: nt, date_fcs
      Integer :: ncid, id_t_dim, varid, xtype, len, attnum
      Integer, Allocatable :: time(:)
      Logical :: status
      Character (Len=13) :: date0, timeres, idate0, idate1
      Character (Len=1) :: Ch
      Character (Len=100) :: t_unit
      Integer :: sys,yy1,mm1,dd1,hh1,yy2,mm2,dd2,hh2, LenSt
      Real*8 :: scale_factor, add_offset, missing_value

      status = nf90_open(fname, nf90_NoWrite, ncid)
      ! Time dimension:
      status = nf90_inq_dimid(ncid,"time", id_t_dim)
      status = nf90_inquire_dimension(ncid, id_t_dim,len=Nt)
      Allocate (Time(Nt))
      status = nf90_inq_varid(ncid,"time", varid)
      status = nf90_get_var(ncid, varid, Time)
      write(idate0,"(A,I10)") " ",time(1); write(idate1,"(A,I10)") " ",time(Nt)
      Deallocate (Time)

      status = nf90_inquire_attribute(ncid, varid, "time", xtype, len, attnum) 
      status = nf90_get_att(ncid, varid, "units", t_unit)

      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+12)
      timeres=" "//t_unit(1:4)
      sys = system ('rm -rf days.txt')
      sys = system ('date -u -d "'//trim(date0)//' '//trim(idate0)//trim(timeres)//' " +%Y%m%d | cat > days.txt')
      Open (1,file="days.txt")
      Read (1,*) date_fcs
      Close(1)

      status = nf90_inq_varid(ncid,"tp", varid)
      status = nf90_get_att(ncid, varid, "scale_factor", scale_factor)
      status = nf90_get_att(ncid, varid, "add_offset", add_offset)
      status = nf90_get_att(ncid, varid, "missing_value", missing_value)
      status = nf90_close(ncid)
   End Subroutine GetTime
   !----------------------------------------------------------------------------------------------
Subroutine wr_ncdf_4d_new (filenm,lon,lat,z,t,nlon,nlat,nz,nt,var,var_nm, &
                       t_unit, var_long_name, var_unit, lon_nm, lat_nm, z_nm, t_nm, flag)
   character (len = *), intent (in) :: filenm, var_nm, t_unit, var_long_name, var_unit
   character (len = *) :: lon_nm, lat_nm, z_nm, t_nm
   integer, parameter :: NDIMS = 4
   integer :: ncid, var_id, dims_id(NDIMS)
   integer :: lon_id, lat_id, z_id, t_id
   integer, intent (in) :: nlon, nlat, nz, nt
   real, dimension (nlon,nlat,nz,nt), intent (in) :: var
   real, intent (in) :: lon(nlon), lat(nlat), flag !, z(nz), t(nt), flag
   integer, intent (in) :: z(nz), t(nt)

   print*,"   Writing netCDF file: ",filenm,"..."
   call hulk ( nf90_create (filenm, nf90_write, ncid) )
        call hulk ( nf90_def_dim (ncid, lon_nm, nlon, dims_id(1)) )
        call hulk ( nf90_def_dim (ncid, lat_nm, nlat, dims_id(2)) )
        call hulk ( nf90_def_dim (ncid, z_nm,   nz,   dims_id(3)) )
        call hulk ( nf90_def_dim (ncid, t_nm,   nt,   dims_id(4)) )

        call hulk ( nf90_def_var (ncid,lon_nm,nf90_float,dims_id(1),lon_id) )
        call hulk ( nf90_def_var (ncid,lat_nm,nf90_float,dims_id(2),lat_id) )
        call hulk ( nf90_def_var (ncid,z_nm,NF90_INT,dims_id(3),z_id) )
        call hulk ( nf90_def_var (ncid,t_nm,NF90_INT,dims_id(4),t_id) )
        call hulk ( nf90_def_var (ncid,var_nm,nf90_float,dims_id(1:4),var_id) )
             call hulk ( nf90_put_att (ncid,lon_id,"long_name","Longitude") )
             call hulk ( nf90_put_att (ncid,lon_id,"units","degree_east") )
             call hulk ( nf90_put_att (ncid,lon_id,"axis","X") )
             call hulk ( nf90_put_att (ncid,lat_id,"long_name","Latitude") )
             call hulk ( nf90_put_att (ncid,lat_id,"units","degree_north") )
             call hulk ( nf90_put_att (ncid,lat_id,"axis","Y") )
             call hulk ( nf90_put_att (ncid,z_id,  "long_name","z") )
             call hulk ( nf90_put_att (ncid,z_id,  "units","Lev") )
             call hulk ( nf90_put_att (ncid,z_id,  "axis","Z") )

             call hulk ( nf90_put_att (ncid,t_id,"long_name","Time") )
             call hulk ( nf90_put_att (ncid,t_id,"units",trim(t_unit)//" 00:00:00") )
             call hulk ( nf90_put_att (ncid,t_id,"calendar","proleptic_gregorian") )
             call hulk ( nf90_put_att (ncid,t_id,"axis","T") )

             call hulk ( nf90_put_att (ncid,var_id,"long_name",trim(var_long_name)))
             call hulk ( nf90_put_att (ncid,var_id,"units", trim(var_unit)))
             call hulk ( nf90_put_att (ncid,var_id,"missing_value",flag) )


        call hulk ( nf90_enddef (ncid) )
        call hulk ( nf90_put_var (ncid,lon_id,lon) )
        call hulk ( nf90_put_var (ncid,lat_id,lat) )
        call hulk ( nf90_put_var (ncid,z_id,  z) )
        call hulk ( nf90_put_var (ncid,t_id,  t) )
        call hulk ( nf90_put_var (ncid,var_id,var) )
   call hulk ( nf90_close (ncid) )
   print*,"   Completed writing file ",filenm

End Subroutine wr_ncdf_4d_new
   !----------------------------------------------------------------------------------------------
Subroutine rd_ncdf_4d_new (filenm,lon,lat,lev,t,nlon,nlat,nlev,nt,var,var_nm)
!   use netcdf
   character (len = *), intent (in) :: filenm, var_nm
   character (len = 15) :: lon_nm, lat_nm, lev_nm, t_nm
   integer, parameter :: NDIMS = 4
   integer :: ncid, var_id, dims_id(NDIMS)
   integer :: lon_id, lat_id, lev_id, t_id
   integer, intent (in out) :: nlon, nlat,  nlev,nt
   real, allocatable, dimension (:,:,:,:), intent (in out) :: var
   real, allocatable, dimension (:), intent (in out) :: lon, lat !, lev, t
   integer, allocatable, dimension (:), intent (in out) :: lev, t

   print*,"   Reading netCDF file: ",filenm,"..."
   status =  nf90_open (filenm,nf90_Nowrite,ncid)
   status =  nf90_inq_varid (ncid,var_nm,var_id)
   status =  nf90_inquire_variable (ncid,var_id,dimids = dims_id(1:4))
   status =  nf90_inquire_dimension (ncid, dims_id(1),name = lon_nm, len = nlon)
   status =  nf90_inquire_dimension (ncid, dims_id(2),name = lat_nm, len = nlat)
   status =  nf90_inquire_dimension (ncid, dims_id(3),name = lev_nm, len = nlev)
   status =  nf90_inquire_dimension (ncid, dims_id(4),name = t_nm, len = nt)

        print*,"      - Longitude's name = ",lon_nm,"; length = ", nlon
        print*,"      - Latitude's name  = ",lat_nm,"; length = ", nlat
        print*,"      - Level's name     = ",lev_nm,"; length = ", nlev
        print*,"      - Time's name      = ",t_nm,"; length = ", nt

   ALLOCATE (lon(nlon),lat(nlat),lev(nlev),t(nt),var(nlon,nlat,nlev,nt))
   status =  nf90_inq_varid (ncid,lon_nm,lon_id)
   status =  nf90_get_var (ncid,lon_id,lon)
   status =  nf90_inq_varid (ncid,lat_nm,lat_id)
   status =  nf90_get_var (ncid,lat_id,lat)

   status =  nf90_inq_varid (ncid,lev_nm,lev_id)
   status =  nf90_get_var (ncid,lev_id,lev)

   status =  nf90_inq_varid (ncid,t_nm,t_id)
   status =  nf90_get_var (ncid,t_id,t)

   status =  nf90_get_var(ncid,var_id,var)
   status =  nf90_close(ncid)
   print*,"   Completed reading file ",filenm

   return

End Subroutine rd_ncdf_4d_new
!-------------------------------------------------------
! <-~~~=| Nf90 checking errors ...
SUBROUTINE hulk (stat)
   use netcdf
   integer :: stat
   if ( stat .ne. nf90_noerr) then
      print*, trim(nf90_strerror(stat))
      stop
   endif
END SUBROUTINE hulk
!----------------------------------------------------------------------------------

   Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
   !Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
   END