Use NetCDF

   Real, Parameter :: Miss=-99.0
   Integer, Parameter :: Nreg=7, Nmon=12, IMiss=-99, Nmem=11, 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"/)
   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(:), Regs(:,:)
   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)

   !----------------------------------------------------------------------------------------------
   Real, Allocatable :: RMEreg(:,:,:,:), RMAEreg(:,:,:,:), CORreg(:,:,:,:)  ! (Nreg,Mon,LT,Mem)
   Real, Allocatable :: RMEwreg(:,:,:,:), RMAEwreg(:,:,:,:), CORwreg(:,:,:,:)  ! (Nreg,Mon,Nws,Mem)

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

   Call rd_ncdf_2d ("/home/tanpv/Data/S2S/Mask/Mask_VN_0.125.nc",lon,lat,nlo,nla,Regs,"mask") ! 7 sub-regs mask
   Deallocate (lon,lat)

   DirInEC  ="/home/tanpv/Data/S2S/EC/"
   DirInTRMM="/home/tanpv/Data/S2S/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")

   Call Very_Score_Reg
   Call Write_Very_Scores_Reg ("Score_daily_reg.nc","Score_weekly_reg.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))
      RME=Miss; RMAE=Miss; COR=Miss; RMEw=Miss; RMAEw=Miss; CORw=Miss
      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
               Do imon=1,Nmon
                  tmp1=0; tmp2=0; tmp3=0; ntmp1=0; ntmp2=0; ntmp3=0
                  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
                        !print*,ifcs,fcsdate(ifcs,imem),iobs,Time_OBS(iobs)
                        !Do ilead=1,Nlead
                        !   tmp1(ilead)=tmp1(ilead)+(EC(ilon,ilat,ilead,ifcs,imem)-TRMM(ilon,ilat,iobs))
                        !Enddo
                        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 Very_Score_Reg
      Integer :: ireg,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)
      Real :: X(Nlon,Nlat,NLead), Y(Nlead)

!   Real, Allocatable :: RMEreg(:,:,:,:), RMAEreg(:,:,:,:), CORreg(:,:,:,:)  ! (Nreg,Mon,LT,Mem)
!   Real, Allocatable :: RMEwreg(:,:,:,:), RMAEwreg(:,:,:,:), CORwreg(:,:,:,:)  ! (Nreg,Mon,Nws,Mem)

      Allocate (RMEreg(Nreg,Nmon,NLead,Nmem), RMAEreg(Nreg,Nmon,NLead,Nmem), CORreg(Nreg,Nmon,NLead,Nmem))
      Allocate (RMEwreg(Nreg,Nmon,Nws,Nmem),RMAEwreg(Nreg,Nmon,Nws,Nmem),CORwreg(Nreg,Nmon,Nws,Nmem))
      RMEreg=Miss; RMAEreg=Miss; CORreg=Miss; RMEwreg=Miss; RMAEwreg=Miss; CORwreg=Miss
      Allocate (Ob(300,NLead),Fc(300,NLead))   ! So mau trong mot thang toi da la 300
      Do ireg=1,Nreg
         Do imem=1,Nmem
            Do imon=1,Nmon
               tmp1=0; tmp2=0; tmp3=0; ntmp1=0; ntmp2=0; ntmp3=0
               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;
                        X(1:Nlon,1:Nlat,1:Nlead)=EC(1:Nlon,1:Nlat,1:NLead,ifcs,imem)
                        Call Areal_Mean (X,Nlon,Nlat,NLead,Y,ireg)
                        Fc(Nspl,1:NLead)=Y(1:NLead)
                        X(1:Nlon,1:Nlat,1:Nlead)=TRMM(1:Nlon,1:Nlat,iobs:iobs+NLead-1)
                        Call Areal_Mean (X,Nlon,Nlat,NLead,Y,ireg)
                        Ob(Nspl,1:NLead)=Y(1:NLead)
                     endif
                  endif
               Enddo
               Call Scores (Ob,Fc,Nspl,Nlead,ME,MAE,CO,MEw,MAEw,COw)
               RMEreg (ireg,imon,1:NLead,imem)=ME (1:NLead)
               RMAEreg(ireg,imon,1:NLead,imem)=MAE(1:NLead)
               CORreg (ireg,imon,1:NLead,imem)=CO(1:NLead)
               RMEwreg (ireg,imon,1:Nws,imem) =MEw (1:Nws)
               RMAEwreg(ireg,imon,1:Nws,imem) =MAEw(1:Nws)
               CORwreg (ireg,imon,1:Nws,imem) =COw (1:Nws)
               print*,"ireg,imem,imon: ",ireg,imem,imon
               !pause
            Enddo
         Enddo
      Enddo
   End Subroutine Very_Score_Reg
   !----------------------------------------------------------------------------------------------
   Subroutine Write_Very_Scores_Reg (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 :: X3d(:,:,:),xreg(:),xmon(:),lead_day(:),lead_week(:)
      Integer :: ireg,ilon,ilat,imon,ilead,imem,ifcs,im,iobs,Nspl,Nlead_day
      Logical :: Add

      t_unit="days since 2000-01-01 "; var_unit="-";lon_nm="reg"; lat_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),xreg(Nreg))
      Allocate (X3d(Nreg,Nmon,Nlead_day))
      Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo
      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
         X3d(:,:,1:Nlead_day)=RMEreg(:,:,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_3d (Fname1,xreg,xmon,lead_day,Nreg,Nmon,Nlead_day,X3d,trim(var_nm), &
                 trim(t_unit),trim(var_ln),trim(var_unit), trim(lon_nm), trim(lat_nm), trim(t_nm), Miss)

            !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_3d_apd (Fname1,X3d,trim(var_nm),Nreg,Nmon,Nlead_day, trim(var_ln), &
                 trim(var_unit),trim(lon_nm), trim(lat_nm), trim(t_nm), Miss)

            !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
         X3d(:,:,1:Nlead_day)=RMAEreg(:,:,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_3d_apd (Fname1,X3d,trim(var_nm),Nreg,Nmon,Nlead_day, trim(var_ln), &
                 trim(var_unit),trim(lon_nm), trim(lat_nm), trim(t_nm), Miss)

         !Call wr_ncdf_3d_apd (Fname1,X3d,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
         X3d(:,:,1:Nlead_day)=CORreg(:,:,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_3d_apd (Fname1,X3d,trim(var_nm),Nreg,Nmon,Nlead_day, trim(var_ln), &
                 trim(var_unit),trim(lon_nm), trim(lat_nm), trim(t_nm), Miss)

         !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(X3d); Allocate (X3d(Nreg,Nmon,Nws))
      Add=.False.; t_nm="week"
      Do imem=1,Nmem
         X3d(:,:,1:Nws)=RMEwreg(:,:,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_3d (Fname2,xreg,xmon,lead_week,Nreg,Nmon,Nws,X3d,trim(var_nm), &
                 trim(t_unit),trim(var_ln),trim(var_unit), trim(lon_nm), trim(lat_nm), trim(t_nm), Miss)

            !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_3d_apd (Fname2,X3d,trim(var_nm),Nreg,Nmon,Nws, trim(var_ln), &
                 trim(var_unit),trim(lon_nm), trim(lat_nm), trim(t_nm), Miss)

            !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
         X3d(:,:,1:Nws)=RMAEwreg(:,:,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_3d_apd (Fname2,X3d,trim(var_nm),Nreg,Nmon,Nws, trim(var_ln), &
                 trim(var_unit),trim(lon_nm), trim(lat_nm), trim(t_nm), Miss)

         !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
         X3d(:,:,1:Nws)=CORwreg(:,:,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_3d_apd (Fname2,X3d,trim(var_nm),Nreg,Nmon,Nws, trim(var_ln), &
                 trim(var_unit),trim(lon_nm), trim(lat_nm), trim(t_nm), Miss)

         !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_Reg
   !----------------------------------------------------------------------------------------------
   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), Ob_Cli(M), Ow_Cli(Nws)
      Real :: ME(NLead),MAE(NLead),COR(NLead), MEw(Nws),MAEw(Nws),CORw(Nws)
      Real :: Corr,prob,tmp1,tmp2

      ME=Miss; MAE=Miss; COR=Miss; MEw=Miss;MAEw=Miss; CORw=Miss
      Do i=1,N   ! Trung binh truot buoc Lw tren tat ca N mau (so lan du bao) 
         X(1:M)=O(i,1:M)
         Call Moving_Ave (X,M,Lw,Y); Ob(i,1:M)=Y(1:M)
         X(1:M)=F(i,1:M)
         Call Moving_Ave (X,M,Lw,Y); Fc(i,1:M)=Y(1:M)
      Enddo
      ! Trung binh quan trac (Ob) tren tat ca N mau cua Trung binh truot:
      Do j=1,M-Lw+1
         X1d(1:N)=Ob(1:N,j)
         Ob_Cli(j)=Ave(X1d, N, Miss)
      Enddo
      ! Sai so tuong doi va he so TQ:
      Do j=1,M-Lw+1
         if (Ob_Cli(j)/=0.AND.Ob_Cli(j)/=Miss) then
            !ME(j)=0; MAE(j)=0; ic=0
            tmp1=0; tmp2=0; ic=0
            Do i=1,N
               if (Fc(i,j)/=Miss.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))
                  tmp1=tmp1+(Fc(i,j)-Ob(i,j))/Ob_Cli(j)
                  tmp2=tmp2+ABS((Fc(i,j)-Ob(i,j))/Ob_Cli(j))
               endif
            Enddo
            !ME(j)=ME(j)/real(ic); MAE(j)=MAE(j)/real(ic)
            ME(j)=tmp1/real(ic); MAE(j)=tmp2/real(ic)
         endif
         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
      ! Trung binh tuan (Lw days in a week):
      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
      ! Trung binh tuan Cli:
      Do j=1,Nws
         X1d(1:N)=Ow(1:N,j)
         Ow_Cli(j)=Ave(X1d, N, Miss)
      Enddo
      ! Sai so va he so TQ:
      Do j=1,Nws
         if (Ow_Cli(j)/=0.AND.Ow_Cli(j)/=Miss) then
            !MEw(j)=0; MAEw(j)=0; ic=0
            tmp1=0; tmp2=0; ic=0
            Do i=1,N
               if (Fw(i,j)/=Miss.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))
                  tmp1=tmp1+(Fw(i,j)-Ow(i,j))/Ow_Cli(j)
                  tmp2=tmp2+ABS((Fw(i,j)-Ow(i,j))/Ow_Cli(j))
               endif
            Enddo
            !MEw(j)=MEw(j)/real(ic); MAEw(j)=MAEw(j)/real(ic)
            MEw(j)=tmp1/real(ic); MAEw(j)=tmp2/real(ic)
         endif
         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 Areal_Mean (X,Nlon,Nlat,NLead,Y,ireg)
      Integer :: Nlon,Nlat,NLead, ireg
      Real :: X(Nlon,Nlat,NLead), Y(NLead)
      ! EC(ilon,ilat,1:NLead,ifcs,imem)
      Integer :: ilon,ilat, ilead
      Real :: tmp, ntmp

      Y(:)=Miss
      Do ilead=1,NLead
         tmp=0; ntmp=0
         Do ilon=1,Nlon
            Do ilat=1,Nlat
               if (Int(Regs(ilon,ilat))==ireg.AND.X(ilon,ilat,ilead)/=Miss) then
                  ntmp=ntmp+1; tmp=tmp+X(ilon,ilat,ilead)
               endif
            Enddo
         Enddo
         Y(ilead)=tmp/ntmp
      Enddo
   End Subroutine Areal_Mean
   !----------------------------------------------------------------------------------------------
   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
      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')
      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
   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"
   Include "/home/tanpv/Mylibs/Ave.f90"
   END
   !---------------------------------------------------------------------------
   Include "/home/tanpv/Mylibs/Student_Prob.f90"