Use NetCDF

   Real, Parameter :: Miss=-99.0
   Integer, Parameter :: Nmon=12, IMiss=-99, Nmem=1, NLead=46, Nws=9, Lw=5  ! Nws=So tuan; Lw = Do dai tuan
   !Character (Len=2), Parameter :: xmem(Nmem)=(/"01","02","03","04","05","06","07","08","09","10","11"/)
   Character (Len=2), Parameter :: xmem(Nmem)=(/"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 :: Nspl(Nmon,Nmem)  !  Nsamples
   Integer :: Nlon, Nlat, Nfcs, date_fcs, Nlon0, Nlat0    ! Nfcs = So lan du bao, date_fcs = yyyymmdd cua ngay DB
   Integer :: Num_Of_Fcs(Nmem)                            ! So lan du bao (Nfc) cua tung Mem
   Integer :: Nlo, Nla, Ntime    ! Nday from 2000-01-01 to 2019-12-31
   !----------------------------------------------------------------------------------------------
   ! Verification Scores:
   Real, Allocatable :: RME(:,:,:,:,:), RMAE(:,:,:,:,:), COR(:,:,:,:,:)  ! (Lon,Lat,Mon,LT,Mem)
   Real, Allocatable :: RMEw(:,:,:,:,:), RMAEw(:,:,:,:,:), CORw(:,:,:,:,:)  ! (Lon,Lat,Mon,Nws,Mem)

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

   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)
   !
   Nspl=0
   !Do imem=1,Nmem
   !   Do ifcs=1,Nfcs
   !      if (fcsdate(ifcs,imem)/=IMiss) then
   !         iobs=0
   !         Do
   !            iobs=iobs+1
   !            if (Time_OBS(iobs)==fcsdate(ifcs,imem)) exit
   !         Enddo
   !         !print*, ifcs,fcsdate(ifcs,imem),iobs,Time_OBS(iobs); ! pause
   !         im=fcsdate(ifcs,imem)/100-fcsdate(ifcs,imem)/10000*100
   !         Nspl(im,imem)=Nspl(im,imem)+1
   !      endif
   !   Enddo
   !   print*,imem,Nspl(:,imem)
   !Enddo
   Call Very_Score
   Call Write_Very_Scores ("Score_daily.nc","Score_weekly.nc")
   CONTAINS
   !----------------------------------------------------------------------------------------------
   Subroutine Very_Score
      Integer :: ilon,ilat,imon,ilead,imem,ifcs,im,iobs,Nspl
      Real :: tmp1(NLead),tmp2(NLead),tmp3(NLead),ntmp1(NLead),ntmp2(NLead),ntmp3(NLead)
      Real,Allocatable :: Ob(:,:),Fc(:,:) !(Nspl,NLead)
      Real :: ME(NLead),MAE(NLead),CO(NLead), MEw(Nws),MAEw(Nws),COw(Nws)
      

      Allocate (RME(Nlon,Nlat,Nmon,NLead,Nmem), RMAE(Nlon,Nlat,Nmon,NLead,Nmem), COR(Nlon,Nlat,Nmon,NLead,Nmem))
      Allocate (RMEw(Nlon,Nlat,Nmon,Nws,Nmem),RMAEw(Nlon,Nlat,Nmon,Nws,Nmem),CORw(Nlon,Nlat,Nmon,Nws,Nmem))
      Allocate (Ob(300,NLead),Fc(300,NLead))   ! So mau trong mot thang toi da la 300
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            Do imem=1,Nmem
               !print*,"ilon,ilat,imem:", ilon,ilat,imem
               Do imon=1,Nmon
                  Nspl=0
                  Do ifcs=1,Num_Of_Fcs(imem) !Nfcs
                     if (fcsdate(ifcs,imem)/=IMiss) then
                        iobs=0
                        Do
                           iobs=iobs+1
                           if (Time_OBS(iobs)==fcsdate(ifcs,imem)) exit
                        Enddo
                        im=fcsdate(ifcs,imem)/100-fcsdate(ifcs,imem)/10000*100   ! Mon
                     endif
                     if (im==imon) then
                        if (iobs<=Ntime-Nlead) then
                           Nspl=Nspl+1;
                           Fc(Nspl,1:NLead)=EC(ilon,ilat,1:NLead,ifcs,imem)
                           Ob(Nspl,1:NLead)=TRMM(ilon,ilat,iobs:iobs+NLead-1)
                        endif
                     endif
                  Enddo
                  Call Scores (Ob,Fc,Nspl,Nlead,ME,MAE,CO,MEw,MAEw,COw)
                  RME (ilon,ilat,imon,1:NLead,imem)=ME (1:NLead)
                  RMAE(ilon,ilat,imon,1:NLead,imem)=MAE(1:NLead)
                  COR (ilon,ilat,imon,1:NLead,imem)=CO(1:NLead)
                  RMEw (ilon,ilat,imon,1:Nws,imem) =MEw (1:Nws)
                  RMAEw(ilon,ilat,imon,1:Nws,imem) =MAEw(1:Nws)
                  CORw (ilon,ilat,imon,1:Nws,imem) =COw (1:Nws)
                  print*,"ilon,ilat,imem,imon: ",ilon,ilat,imem,imon
                  !pause
               Enddo
            Enddo
         Enddo
      Enddo
   End Subroutine Very_Score   
   !----------------------------------------------------------------------------------------------
   Subroutine Write_Very_Scores (Fname1,Fname2)
      Character (Len=*) :: Fname1,Fname2
      Character (Len=255) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, var_nm
      Real, Allocatable :: X4d(:,:,:,:),xmon(:),lead_day(:),lead_week(:)
      Integer :: ilon,ilat,imon,ilead,imem,ifcs,im,iobs,Nspl,Nlead_day
      Logical :: Add

      t_unit="days since 2000-01-01 "; var_unit="-";lon_nm="lon"; lat_nm="lat"; z_nm="mon"; t_nm="day"
     
      Nlead_day=Nlead-Lw+1   ! So ngay leadtime sau khi trung binh truot Lw tinh tu ngay du bao (LT=0)
      Allocate (xmon(Nmon),lead_day(Nlead_day),lead_week(Nws))
      Allocate (X4d(Nlon,Nlat,Nmon,Nlead_day))
      Do imon=1,Nmon; xmon(imon)=imon; Enddo
      Do ilead=1,Nlead_day; lead_day(ilead)=ilead; Enddo
      Do ilead=1,Nws; lead_week(ilead)=ilead; Enddo
      ! Relative Mean Errors:
      Add=.False.
      Do imem=1,Nmem
         X4d(:,:,:,1:Nlead_day)=RME(:,:,:,1:Nlead_day,imem); 
         write(var_nm,"(A,I2.2)") "rme_",imem; write(var_ln,"(A,I2.2)") "Relative mean error of member ",imem
         if (.not.Add) then
            Call wr_ncdf_4d (Fname1,lon,lat,xmon,lead_day,nlon,nlat,Nmon,Nlead_day,X4d,trim(var_nm), &
                       trim(t_unit), trim(var_ln), trim(var_unit), trim(lon_nm), trim(lat_nm), trim(z_nm), trim(t_nm), Miss)
            Add=.True.
         else
            Call wr_ncdf_4d_apd (Fname1,X4d,trim(var_nm), nlon,nlat,Nmon,Nlead_day, &
                           trim(var_ln), trim(var_unit), trim(lon_nm), trim(lat_nm), trim(z_nm), trim(t_nm), Miss)
         endif
      Enddo
      ! Relative Absolute Mean Errors:
      Do imem=1,Nmem
         X4d(:,:,:,1:Nlead_day)=RMAE(:,:,:,1:Nlead_day,imem); 
         write(var_nm,"(A,I2.2)") "rmae_",imem; write(var_ln,"(A,I2.2)") "Relative absolute mean error of member ",imem
         Call wr_ncdf_4d_apd (Fname1,X4d,trim(var_nm), nlon,nlat,Nmon,Nlead_day, &
                           trim(var_ln), trim(var_unit), trim(lon_nm), trim(lat_nm), trim(z_nm), trim(t_nm), Miss)
      Enddo
      ! Correlation Coefficients:
      Do imem=1,Nmem
         X4d(:,:,:,1:Nlead_day)=COR(:,:,:,1:Nlead_day,imem);
         write(var_nm,"(A,I2.2)") "cor_",imem; write(var_ln,"(A,I2.2)") "Correlation coefficients of member ",imem
         Call wr_ncdf_4d_apd (Fname1,X4d,trim(var_nm), nlon,nlat,Nmon,Nlead_day, &
                           trim(var_ln), trim(var_unit), trim(lon_nm), trim(lat_nm), trim(z_nm), trim(t_nm), Miss)
      Enddo
      !--------------------------------------
      ! Relative Mean Errors for weely data:
      Deallocate(X4d); Allocate (X4d(Nlon,Nlat,Nmon,Nws))
      Add=.False.; t_nm="week"
      Do imem=1,Nmem
         X4d(:,:,:,1:Nws)=RMEw(:,:,:,1:Nws,imem);
         write(var_nm,"(A,I2.2)") "rme_",imem; write(var_ln,"(A,I2.2)") "Relative mean error of member ",imem
         if (.not.Add) then
            Call wr_ncdf_4d (Fname2,lon,lat,xmon,lead_week,nlon,nlat,Nmon,Nws,X4d,trim(var_nm), &
                       trim(t_unit), trim(var_ln), trim(var_unit), trim(lon_nm), trim(lat_nm), trim(z_nm), trim(t_nm), Miss)
            Add=.True.
         else
            Call wr_ncdf_4d_apd (Fname2,X4d,trim(var_nm), nlon,nlat,Nmon,Nws, &
                           trim(var_ln), trim(var_unit), trim(lon_nm), trim(lat_nm), trim(z_nm), trim(t_nm), Miss)
         endif
      Enddo
      ! Relative Absolute Mean Errors for weely data:
      Do imem=1,Nmem
         X4d(:,:,:,1:Nws)=RMAEw(:,:,:,1:Nws,imem);
         write(var_nm,"(A,I2.2)") "rmae_",imem; write(var_ln,"(A,I2.2)") "Relative absolute mean error of member ",imem
         Call wr_ncdf_4d_apd (Fname2,X4d,trim(var_nm), nlon,nlat,Nmon,Nws, &
                           trim(var_ln), trim(var_unit), trim(lon_nm), trim(lat_nm), trim(z_nm), trim(t_nm), Miss)
      Enddo
      ! Correlation Coefficients for weely data:
      Do imem=1,Nmem
         X4d(:,:,:,1:Nws)=CORw(:,:,:,1:Nws,imem);
         write(var_nm,"(A,I2.2)") "cor_",imem; write(var_ln,"(A,I2.2)") "Correlation coefficients of member ",imem
         Call wr_ncdf_4d_apd (Fname2,X4d,trim(var_nm), nlon,nlat,Nmon,Nws, &
                           trim(var_ln), trim(var_unit), trim(lon_nm), trim(lat_nm), trim(z_nm), trim(t_nm), Miss)
      Enddo
   End Subroutine Write_Very_Scores
   !----------------------------------------------------------------------------------------------
   Subroutine Scores (O,F,N,M,ME,MAE,COR,MEw,MAEw,CORw)
      Integer :: N,M, i,j,k,i1,j1, ic ! N=Nspl; M=Nlead
      Real :: O(N,M),F(N,M),Ow(N,Nws),Fw(N,Nws), X(M), Y(M)
      Real :: Ob(N,M),Fc(N,M),X1d(N),Y1d(N)
      Real :: ME(NLead),MAE(NLead),COR(NLead), MEw(Nws),MAEw(Nws),CORw(Nws)
      Real :: Corr,prob,tmp1,tmp2
      Do i=1,N
         X(1:M)=O(i,1:M)
         Call Moving_Ave (X,N,Lw,Y); Ob(i,1:M)=Y(1:M)
         X(1:M)=F(i,1:M)
         Call Moving_Ave (X,N,Lw,Y); Fc(i,1:M)=Y(1:M)
      Enddo
      Do j=1,M-Lw+1
         ME(j)=0; MAE(j)=0; ic=0
         Do i=1,N
            if (Ob(i,j)/=0.AND.Ob(i,j)/=Miss) then
               ic=ic+1
               ME(j)=ME(j)+(Fc(i,j)-Ob(i,j))/Ob(i,j)
               MAE(j)=MAE(j)+ABS((Fc(i,j)-Ob(i,j))/Ob(i,j))
            endif
         Enddo
         ME(j)=ME(j)/real(ic); MAE(j)=MAE(j)/real(ic)
         X1d(1:N)=Fc(1:N,j); Y1d(1:N)=Ob(1:N,j)
         Call Correlation(X1d,Y1d,N,Corr,prob,Miss)
         COR(j)=Corr
      Enddo
      Do i=1,N
         X(1:M)=O(i,1:M); Y(1:M)=F(i,1:M)
         k=0
         Do j=2,M,Lw  ! 2:Nlead,5
            tmp1=0; tmp2=0
            Do j1=j,j+Lw-1
               tmp1=tmp1+X(j1); tmp2=tmp2+Y(j1)
            Enddo
            k=k+1; Ow(i,k)=tmp1/Lw; Fw(i,k)=tmp2/Lw
         Enddo
      Enddo
      Do j=1,Nws
         MEw(j)=0; MAEw(j)=0; ic=0
         Do i=1,N
            if (Ow(i,j)/=0.AND.Ow(i,j)/=Miss) then
               ic=ic+1
               MEw(j)=MEw(j)+(Fw(i,j)-Ow(i,j))/Ow(i,j)
               MAEw(j)=MAEw(j)+ABS((Fw(i,j)-Ow(i,j))/Ow(i,j))
            endif
         Enddo
         MEw(j)=MEw(j)/real(ic); MAEw(j)=MAEw(j)/real(ic)
         X1d(1:N)=Fw(1:N,j); Y1d(1:N)=Ow(1:N,j)
         Call Correlation(X1d,Y1d,N,Corr,prob,Miss)
         CORw(j)=Corr
      Enddo
   End Subroutine Scores
   !----------------------------------------------------------------------------------------------
   Subroutine Read_EC
      Integer :: imem, ifcs
      Real, Allocatable :: xlon(:), xlat(:)
      Integer, Allocatable :: xlead(:), xfcsdate(:)
      Integer :: Nfc

      Do imem=1,Nmem !imem=1,Nmem
         Fname=trim(DirInEC)//"Prec_daily_s2s_ecmwf_"//xmem(imem)//".nc" 
         Call rd_ncdf_4d_new (Fname,xlon,xlat,xlead,xfcsdate,Nlon,Nlat,NLd,Nfc,X4d,"pr")   ! lead(:): 1-46 days; fcsdate(:): in format of yyyymmdd
         print*,imem,Nlon,Nlat,NLd,Nfc,Nmem 
         if (imem==1) then
            Nfcs=Nfc
            Allocate (EC(Nlon,Nlat,NLd,Nfcs,Nmem))
            Allocate (Lon(Nlon),Lat(Nlat),Lead(NLd),fcsdate(Nfcs,Nmem))
            EC=Miss; fcsdate=IMiss
            Lon(:)=xlon(:); Lat(:)=xlat(:); Lead(:)=xlead(:)
         endif
         Num_Of_Fcs(imem)=Nfc
         print*,"====> Nfcs,Num_Of_Fcs(imem): ",Nfcs,Num_Of_Fcs(imem)
         fcsdate(1:Nfc,imem)=xfcsdate(1:Nfc)
         EC(:,:,:,1:Nfc,imem) = X4d(:,:,:,1:Nfc)
         Deallocate (xlon,xlat,xlead,xfcsdate,X4d)
      Enddo
   End Subroutine Read_EC
   !----------------------------------------------------------------------------------------------
   Subroutine Read_TRMM (FName)
      Character (Len=*) :: Fname
      Real, Allocatable :: X3d(:,:,:),xlon(:),xlat(:),xtime(:)
      Integer :: id,ilon,ilat,it

      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,Ntime,X3d,"precipitation")
      !print*,"TRMM - nlo,nla,Ntime: ", nlo,nla,Ntime; pause
      !print*,shape(xlon)
      !print*,shape(xlat);
      !print*,shape(TRMM); pause
      Allocate (TRMM(Nla,Nlo,Ntime))
      Do ilon=1,Nlo
         Do ilat=1,Nla
            TRMM(ilat,ilon,:)=X3d(ilon,ilat,:)
         Enddo
      Enddo
      print*,"Reading "//trim(Fname)//" ... done!"
      sys = system ('rm -rf /home/tanpv/Progs/S2S/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 >> /home/tanpv/Progs/S2S/days.txt')
         print*,id,Ntime
      Enddo
      Allocate (Time_OBS(Ntime))
      Open (1,file="/home/tanpv/Progs/S2S/days.txt",Status="old")
      Do id=1,Ntime; Read (1,*) Time_OBS(id); Enddo
      print*,"Calling Read_TRMM ... done!"
   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
!----------------------------------------------------------------------------------
      Subroutine Correlation(X,Y,N,Corr,prob,Miss)
        Integer :: N, N1, i, j
        Real :: X(n), Y(n), Miss, Xtb, Ytb, Sx, Sy, Rxy, tmp, t

        N1=0; Xtb=0; Ytb=0; Sx=0; Sy=0; Rxy=0
        Do i=1,N
           if (X(i)/=Miss.AND.Y(i)/=Miss) then
              n1=n1+1
              Xtb=Xtb+X(i); Sx=Sx+X(i)*X(i)
              Ytb=Ytb+Y(i); Sy=Sy+Y(i)*Y(i)
              Rxy=Rxy+X(i)*Y(i)
           endif
        Enddo
        if (N1>2) then
           Xtb=Xtb/Real(N1); Ytb=Ytb/Real(N1)
           Sx=Sx/Real(N1)-Xtb*Xtb; Sy=Sy/Real(N1)-Ytb*Ytb
           tmp=Rxy/Real(N1)
           Rxy=(tmp-Xtb*Ytb)/Sqrt(Sx*Sy)
           t = ABS(Rxy/ ( Sqrt(1-Rxy**2) / Sqrt(Real(N1-2)) ) )
           Prob = T_Dist (t, Real(N1-2))
        else
           Rxy=Miss; prob=Miss
        endif
        Corr=Rxy

      End Subroutine Correlation
      !---------------------------------------------------------------------------------
      Subroutine Moving_Ave (X,N,m,Y)
         Integer :: N,m, i,j,k   ! m must be an odd number
         Real :: X(N), Y(N), tmp

         Do j=1,N-m+1
            tmp=0
            Do i=j,j+m-1
               tmp=tmp+X(i)
            Enddo
            Y(j)=tmp/Real(m)
         Enddo
      End Subroutine Moving_Ave
      !---------------------------------------------------------------------------------
   Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
   Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
   END
   !---------------------------------------------------------------------------
   Include "/home/tanpv/Mylibs/Student_Prob.f90"