! Chuong trinh nay dung de hieu chinh mua du bao noi mua (s2s) cua ECMWF
! cho ca hai truong hop training va testing
! Qua trinh training la xay dung ham TP=f(X1,...,Xm) su dung mang ANN mot lop an, mot lop Input, mot lop Output
! Version hien tai chi dung cho 1-2 thang co dinh ma duoc xem la dong nhat ve quy luat hieu chinh
! Neu training thi chuong trinh se tra ve ket qua hieu chinh dong thoi luu lai file trong so va cac tham so cua ANN
! Neu teating thi chuong trinh se doc thong tin trong so va cac tham so cua ANN tu file da luu va tra ve ket qua hieu chinh.
! Cac tham so can chu y:
! 1) Yr1, Yr2: Cac nam duoc tinh
! 2) Lw: Kich thuoc cua so thoi gian se truot (tong hoac trung binh truot)
! 3) So phien du bao trong tap so lieu (vi du, moi thang 9 phien x 2 thang = 18 phien)
!    (cac ngay ung voi cac phien du bao duoc luu trong file hdate.txt)

   Use NetCDF

   Real, Parameter :: Miss=-99.0
   Integer, Parameter :: Nmon=12, IMiss=-99, NLead=46, Lw=5, Nws=NLead-Lw+1  ! Nws=So tuan; Lw = Do dai tuan
   Integer, Parameter :: Yr1=2001, Yr2=2020, Nfcs=18, Nyrs=(Yr2-Yr1+1), Nspl=Nyrs*Nfcs    ! 2 thang x 9 lan x 20 nam
   Integer, Allocatable :: time(:)

   Integer :: LenSt, sys
   Character (Len=255) :: Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2
   Character (Len=255) :: DirInEC, DirInObs, DirInStn 
   Character (Len=12)  :: date0, timeres, idate0, idate1
   Character (Len=5)   :: hdate(Nfcs)     ! forecast dates
   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(:,:,:,:),OBS(:,:,:,:)        ! (Nlon,Nlat,Nspl,NLead)
   Real, Allocatable :: EC_bc(:,:,:,:)                  ! (Nlon,Nlat,Nspl,NLead)
   Real, Allocatable :: EC_w(:,:,:,:),Obs_w(:,:,:,:)    ! (Nlon,Nlat,Nspl,Nws) - Tong/Trung binh tuan

   Real, Allocatable, Dimension(:,:,:)   :: TP_cli, Obs_cli  ! (Nlon,Nlat,NLead) - Clim over all forecasts
   Real, Allocatable, Dimension(:,:,:,:) :: LH, SH, OLR, TX, TM 
   !----------------------------------------------------------------------------------------------
    Integer :: ireg, ilead
    Integer :: Nlayer                             ! Num of layers including Input, output and hidden laysers
    Integer :: Nepoch                             ! Num of epochs
    double precision :: learn_rate ! Learning_Rate
    Integer, Allocatable :: Layer_size(:)                      ! (Nlayer)
   !----------------------------------------------------------------------------------------------   
   Real, Allocatable :: xlon(:), xlat(:), xtime(:)
   Real, Allocatable :: X4d(:,:,:,:), X3d(:,:,:) 
   Integer :: Nlon, Nlat, 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
   !----------------------------------------------------------------------------------------------
   ! Verification Scores:
   Real, Allocatable :: R_ME(:,:,:), R_MAE(:,:,:), COR(:,:,:), Prb(:,:,:)  ! (Lon,Lat,LT)
   Real, Allocatable :: MEe(:,:,:), MAEe(:,:,:), CORe(:,:,:), Prbe(:,:,:)  ! (Lon,Lat,LT)

   !----------------------------------------------------------------------------------------------
   Real, Allocatable :: RMEreg(:), RMAEreg(:), CORreg(:)  ! (LT)

   Integer :: imem, ifcs, iy, im, id, NLd, iobs
   !--------------------------------------------------------------------------------
   ! ANN:
   character (len=250) :: input_file, Out_file, save_file
   character (len=250) :: Dir, Dir_ANN_Out, Dir_Cli, Dir_ANN_Save, Dir_Veri, Dir_Corrected
   logical :: training_mode

   ! Use for calculating scaling factors
   Integer :: Nsample, Nepoc, Nsample1
   Character (Len=5)   :: yvarn, Model_nm

   Real,allocatable, Dimension(:,:) :: Inputs, Outputs, Inputs_X, Outputs_Y
   Real,allocatable, Dimension(:)   :: min_in, max_in, min_out, max_out
   Real :: mean_error, mean_error_check
   !-------------------------------------------------------------------------------
    integer ::  Nin, Nhidden, Nout, iin, ihid, iout
    integer :: i,j,k

    Real,allocatable, Dimension(:,:) :: weight_ih, weight_ho, moment_ih, moment_ho
    Real,allocatable, Dimension(:)   :: fin, fhid, fout, ftarget, delta_o, delta_h
 
    Real, Allocatable   :: X_inp(:,:), Y_inp(:), Yest(:,:), Yhq(:,:)

   ! Training
    Real :: learning_rate0, learning_rate,momentum, error, tmp, mintmp
    logical :: debug
   !-------------------------------------------------------------------------------
   ! controller_save_info:
   ! Nin,Nhidden,Nout   - Nin=MXvar; Nout=MYvar
   Integer :: MXvar, MYvar, Mx, Mvars 
   Integer :: Ncase,icase
   Real,allocatable, Dimension(:,:,:) :: weight_ip    ! (Ncase,Nin,Nhidden) - Ncase=So lan tinh, 
   Real,allocatable, Dimension(:,:,:) :: weight_op    ! (Ncase,Nhidden,Nout)
   Real,allocatable, Dimension(:,:)   :: min_ip, max_ip, min_op, max_op  ! (Ncase)
   ! For saving in ANN subroutine:
   Real,allocatable, Dimension(:,:)   :: wei_i, wei_o
   Real,allocatable, Dimension(:)     :: min_i, max_i, min_o, max_o
   !------------------- For ANN -------------------------------------------------------------
   !Namelist /ann_info/   Nlayer,learn_rate,Nepoch
   !Namelist /layer_info/ Layer_size
   Namelist /ann_info/ Nhidden,learning_rate,momentum,training_mode,Nepoc,debug
   Namelist /Model_info/ Model_nm,yvarn
   !--------------------------------------------------------------------------------

    Open (1, file="Namelist_S2S_ecmwf.txt",Status="old")
    Read (1, ann_info)
    Read (1, Model_info)
   ! Allocate (Layer_size(Nlayer))
   ! Read (1, layer_info)
    Close (1)

   Fname="/home/tanpv/Data/S2S/2024/Mask_R6_0.125_S2S.nc"
   Call rd_ncdf_2d (trim(Fname),lon,lat,Nlon,Nlat,Regs,"r6") ! 7 sub-regs mask
   Fname="/home/cordex_sea/Chi/hdate.txt"
   Open (1,file=trim(Fname), status="old")
   Do i=1,Nfcs
      read (1,*) hdate(i)
   Enddo
   Close (1)

   DirInEC       ="/home/cordex_sea/Chi/TayNguyen/Store/"
   DirInObs      ="/home/cordex_sea/Chi/TayNguyen/CHIRPS/"
   Dir_ANN_Save  ="/home/cordex_sea/Chi/TayNguyen/ANN_Save/"
   Dir_Cli       ="/home/cordex_sea/Chi/TayNguyen/Cli_Out/" !"/home/tanpv/Data/562_VNU/0.25_VN/Cli_Out/"
   Dir_Corrected ="/home/cordex_sea/Chi/Res/"

   Call Read_data   ! Read all data of TP, LH, SH, OLR, TX, TM from entire domain (no mask for Central Highland), and CHIRPs
   if (training_mode) then
      MXvar = 22; MYvar = 1
      Mvars = MXvar
      Ncase = Npoint()   ! Obtain from Function Npoint
      Allocate (weight_ip(1:Ncase,0:MXvar,1:Nhidden),weight_op(Ncase,0:Nhidden,MYvar))
      Allocate (min_ip(Ncase,MXvar),max_ip(Ncase,MXvar),min_op(Ncase,MYvar),max_op(Ncase,MYvar))
   else
      Fname=trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//".txt"
      Ncase = Npoint()
      Call Read_Ann_info (Fname)
      Mvars = MXvar
   endif
   Allocate (wei_i(0:MXvar,1:Nhidden), wei_o(0:Nhidden,MYvar))     ! To be used in ANN Sub
   Allocate (min_i(MXvar),max_i(MXvar),min_o(MYvar),max_o(MYvar))  ! To be used in ANN Sub

   Allocate (EC_bc(Nlon,Nlat,Nspl,Nlead))
   icase = 0
   Do ilon=22,31 !1,Nlon
      Do ilat=18,27 !1,Nlat
         if (Regs(ilon,ilat)>0) then   ! R6
            Do ilead=1,Nlead
               icase = icase+1
               Call Make_data (ilon,ilat,ilead)
               !write (save_file,"(A,I2.2,A)") trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//"_lead_",ilead,".txt"
               Call ANN_2025(Inputs_X,Outputs_Y,Nsample1,Mvars,1,Miss)    !ANN_2023 (Xinp,Youtp,Nspl,MXvar,MYvar,Miss)
               if (training_mode) then
                  Call Save_Ann_info (icase)
               else
                  Call Load_Ann_info (icase)
               endif
               EC_bc(ilon,ilat,1:Nspl,ilead) = Yest(1:Nspl,1)             ! Ann outputs
            Enddo
         endif
      Enddo
   Enddo

   Call Calc_weekly_data 
   Fname=trim(Dir_Corrected)//"S2S_ecmwf_tp_Ann_old_training.nc"
   Call Write_BC_Outputs (Fname)
   if (training_mode) then
      Fname=trim(Dir_ANN_Save)//"Sav_"//trim(yvarn)//".txt"
      Call Write_Ann_info (Fname)
   endif
   CONTAINS
   !----------------------------------------------------------------------------------------------
   Subroutine Save_Ann_info (icase)
      Integer :: icase
      ! Nin=MXvar; Nout=MYvar wei_i, wei_o
      weight_ip(icase,0:MXvar,1:Nhidden) = wei_i(0:MXvar,1:Nhidden)
      weight_op(icase,0:Nhidden,1:MYvar) = wei_o(0:Nhidden,1:MYvar)
      min_ip(icase,1:MXvar) = min_i(1:MXvar); max_ip(icase,1:MXvar) = max_i(1:MXvar)
      min_op(icase,1:MYvar) = min_o(1:MYvar); max_op(icase,1:MYvar) = max_o(1:MYvar)
   End Subroutine Save_Ann_info
   !----------------------------------------------------------------------------------------------
   Subroutine Load_Ann_info (icase)
      Integer :: icase
      ! Nin=MXvar; Nout=MYvar wei_i, wei_o
      wei_i(0:MXvar,1:Nhidden) = weight_ip(icase,0:MXvar,1:Nhidden)
      wei_o(0:Nhidden,1:MYvar) = weight_op(icase,0:Nhidden,1:MYvar)
      min_i(1:MXvar) = min_ip(icase,1:MXvar); max_i(1:MXvar) = max_ip(icase,1:MXvar)
      min_o(1:MYvar) = min_op(icase,1:MYvar); max_o(1:MYvar) = max_op(icase,1:MYvar)
   End Subroutine Load_Ann_info
   !----------------------------------------------------------------------------------------------
   Subroutine Read_Ann_info (Fname)
      Character (Len=*) :: Fname

      Open (1,file=trim(Fname),form="unformatted",status="old")
      read(1) MXvar,Nhidden,MYvar,Ncase

      Allocate (weight_ip(1:Ncase,0:MXvar,1:Nhidden),weight_op(Ncase,0:Nhidden,MYvar))
      Allocate (min_ip(Ncase,MXvar),max_ip(Ncase,MXvar),min_op(Ncase,MYvar),max_op(Ncase,MYvar))

      read(1) weight_ip
      read(1) weight_op
      read(1) min_ip
      read(1) max_ip
      read(1) min_op
      read(1) max_op
      Close (1)
   End Subroutine Read_Ann_info
   !----------------------------------------------------------------------------------------------
   Subroutine Write_Ann_info (Fname)
      Character (Len=*) :: Fname

      Open (1,file=trim(Fname),form="unformatted",status="unknown")
      write(1) MXvar,Nhidden,MYvar,Ncase
      write(1) weight_ip
      write(1) weight_op
      write(1) min_ip
      write(1) max_ip
      write(1) min_op
      write(1) max_op
      Close (1)
   End Subroutine Write_Ann_info
   !----------------------------------------------------------------------------------------------
   Integer Function Npoint ()
      Integer :: ilon,ilat,ilead, i

      i=0
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Regs(ilon,ilat)>0) then
               Do ilead=1,Nlead
                  i=i+1
               Enddo
            endif
         Enddo
      Enddo
      Npoint=i
   End Function Npoint
   !----------------------------------------------------------------------------------------------
   Subroutine Make_data (ilon,ilat,ilead)
      Integer :: ilon,ilat,ilead,it,i,j

      ! OBS, EC, LH, SH, OLR, TX, TM
      Nsample1=Nspl  ! Dung luong mau
      Mx=Mvars                 ! So bien doc lap
      if (Allocated(X_inp)) Deallocate (X_inp)
      if (Allocated(Y_inp)) Deallocate (Y_inp)
      Allocate (X_inp(Nsample1,Mx),Y_inp(Nsample1))

      j = 3;  ! Vi tri diem luoi cach tam o luoi hien tai
      Do i=1,Nspl
         Y_inp(i)    = OBS(ilon,ilat,i,ilead)
         X_inp(i,1)  = EC (ilon,ilat,i,ilead)
         X_inp(i,2)  = LH (ilon,ilat,i,ilead)
         X_inp(i,3)  = SH (ilon,ilat,i,ilead)
         X_inp(i,4)  = OLR(ilon,ilat,i,ilead)
         X_inp(i,5)  = TX (ilon,ilat,i,ilead)
         X_inp(i,6)  = TM (ilon,ilat,i,ilead)

         X_inp(i,7)  = EC (ilon-j,ilat,i,ilead)
         X_inp(i,8)  = EC (ilon+j,ilat,i,ilead)
         X_inp(i,9)  = EC (ilon,ilat-j,i,ilead)
         X_inp(i,10) = EC (ilon,ilat+j,i,ilead)

         X_inp(i,11) = LH (ilon-j,ilat,i,ilead)
         X_inp(i,12) = LH (ilon+j,ilat,i,ilead)
         X_inp(i,13) = LH (ilon,ilat-j,i,ilead)
         X_inp(i,14) = LH (ilon,ilat+j,i,ilead)

         X_inp(i,15) = SH (ilon-j,ilat,i,ilead)
         X_inp(i,16) = SH (ilon+j,ilat,i,ilead)
         X_inp(i,17) = SH (ilon,ilat-j,i,ilead)
         X_inp(i,18) = SH (ilon,ilat+j,i,ilead)

         X_inp(i,19) = OLR(ilon-j,ilat,i,ilead)
         X_inp(i,20) = OLR(ilon+j,ilat,i,ilead)
         X_inp(i,21) = OLR(ilon,ilat-j,i,ilead)
         X_inp(i,22) = OLR(ilon,ilat+j,i,ilead)
      Enddo
      Nsample1=Nspl
      if (Allocated(Inputs_X)) Deallocate (Inputs_X)
      if (Allocated(Outputs_Y)) Deallocate (Outputs_Y)
      Allocate (Inputs_X(Mx,Nsample1), Outputs_Y(1,Nsample1))
      i=0
      Do it=1,Nspl
         i=i+1
         Outputs_Y(1,i)=Y_inp(it)
         Do j=1,Mx
            Inputs_X(j,i)=X_inp(it,j)
         Enddo
      Enddo
   End Subroutine Make_data
   !----------------------------------------------------------------------------------------------
   Subroutine Write_BC_Outputs (Fname)
      Character (Len=*) :: Fname
      Character (Len=255) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, var_nm
      Real, Allocatable :: X4d(:,:,:,:),xmon(:),lead(:),sampl(:)
      Integer :: ilon,ilat,imon,ilead,imem,ifcs,im,iobs, i,j,iyr, sys1,sys2
      Integer :: dfc(Nfcs), mfc(Nfcs)   ! forecast date and month - hdate(Nfcs) - Yr1=2001, Yr2=2020
      Integer :: days_diff, jd_ref, jd_target

      t_unit="days since 1980-01-01 "; 
      var_unit="mm/day"; lon_nm="lon"; lat_nm="lat"; z_nm="fcs_time"; t_nm="week"

      Allocate (lead(Nlead),sampl(Nspl))
      Do ilead=1,Nlead; lead (ilead)=ilead; Enddo
      Do ispl =1,Nspl;  sampl(ilead)=ispl; Enddo

      Do i=1,Nfcs
         read(hdate(i),"(I2)") mfc(i); read(hdate(i),"(3x,I2)") dfc(i)
      Enddo
      i=0
      Do iyr=Yr1,Yr2
         Do j=1,Nfcs
            i=i+1;
            jd_ref = julian_day(1980, 1, 1) 
            jd_target = julian_day(iyr, mfc(j), dfc(j))
            sampl(i)=jd_target - jd_ref !iyr*10000+mfc(j)*100+dfc(j)
         Enddo
      Enddo
      var_ln="Observation precipitation (CHIRPS)"
      Call wr_ncdf_4d (Fname,Lon,Lat,sampl,lead,Nlon,Nlat,Nspl,Nlead,Obs_w,"obs", &
                       t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
      var_ln="S2S precipitation of ecmwf forecast"
      Call wr_ncdf_4d_apd (Fname,EC_w,"fcs", Nlon,Nlat,Nspl,Nlead, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
   End Subroutine Write_BC_Outputs
   !----------------------------------------------------------------------------------------------
   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
         if (ntmp>0) Y(ilead)=tmp/ntmp
      Enddo
   End Subroutine Areal_Mean
   !----------------------------------------------------------------------------------------------
   Subroutine Read_data
      Real, Parameter :: ec_miss=-32767.0, obs_miss=-9999.0
      Integer :: imem, ifcs
      Real, Allocatable :: xlon(:), xlat(:), xtime(:), X3d(:,:,:)
      Integer :: Nlo,Nla,NLt,i,ifc,ispl,iyr

      ! Real, Allocatable, Dimension(:,:,:,:) :: LH, SH, OLR, TX, TM
      ! Read precipitation forecast data:
      Allocate (EC(Nlon,Nlat,Nspl,Nlead))
      Allocate (LH(Nlon,Nlat,Nspl,Nlead), SH(Nlon,Nlat,Nspl,Nlead), OLR(Nlon,Nlat,Nspl,Nlead))
      Allocate (TX(Nlon,Nlat,Nspl,Nlead), TM(Nlon,Nlat,Nspl,Nlead))
      ispl=0
      Do iyr=Yr1,Yr2
         Do ifc=1,Nfcs
            write(Fname,"(A,I4.4,A)") trim(DirInEC)//"LH_SH_OLR_TP_cf_",iyr,"-"//hdate(ifc)//".nc"
            Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"tp") 
            Where (X3d==ec_miss) X3d=Miss
            ispl=ispl+1
            EC(:,:,ispl,:) = X3d(:,:,:)
            Deallocate (xlon,xlat,xtime,X3d)
            Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"slhf")  ! LH
            Where (X3d==ec_miss) X3d=Miss
            LH(:,:,ispl,:) = X3d(:,:,:)
            Deallocate (xlon,xlat,xtime,X3d)
            Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"sshf")  ! SH
            Where (X3d==ec_miss) X3d=Miss
            SH(:,:,ispl,:) = X3d(:,:,:)
            Deallocate (xlon,xlat,xtime,X3d)
            Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"ttr")  ! OLR
            Where (X3d==ec_miss) X3d=Miss
            OLR(:,:,ispl,:) = X3d(:,:,:)
            Deallocate (xlon,xlat,xtime,X3d)
         Enddo
      Enddo
      ispl=0
      Do iyr=Yr1,Yr2
         Do ifc=1,Nfcs
            write(Fname,"(A,I4.4,A)") trim(DirInEC)//"TX_TM_cf_",iyr,"-"//hdate(ifc)//".nc"
            Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"mx2t6")
            Where (X3d==ec_miss) X3d=Miss
            ispl=ispl+1
            TX(:,:,ispl,:) = X3d(:,:,:)
            Deallocate (xlon,xlat,xtime,X3d)
            Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"mn2t6")  ! LH
            Where (X3d==ec_miss) X3d=Miss
            TM(:,:,ispl,:) = X3d(:,:,:)
            Deallocate (xlon,xlat,xtime,X3d)
         Enddo
      Enddo

      ! Read precipitation obs data:
      Allocate (OBS(Nlon,Nlat,Nspl,Nlead))
      ispl=0
      Do iyr=Yr1,Yr2
         Do ifc=1,Nfcs
            write(Fname,"(A,I4.4,A)") trim(DirInObs)//"chirps_",iyr,"-"//hdate(ifc)//".nc"
            Call rd_ncdf_3d (Fname,xlon,xlat,xtime,nlo,nla,nLt,X3d,"precip")
            Where (X3d==obs_miss) X3d=Miss
            ispl=ispl+1
            OBS(:,:,ispl,:) = X3d(:,:,:)
            Deallocate (xlon,xlat,xtime,X3d)
         Enddo
      Enddo
   End Subroutine Read_data
   !----------------------------------------------------------------------------------------------
   Subroutine Calc_weekly_data
      ! Real, Allocatable :: EC_w(:,:,:,:),Obs_w(:,:,:,:)    ! (Nlon,Nlat,Nspl,Nws) - Tong/Trung binh tuan
      Integer :: ilon,ilat,ilead,ispl
      Real    :: X(Nlead), Y(Nlead)

      Allocate (EC_w(Nlon,Nlat,Nspl,Nlead),Obs_w(Nlon,Nlat,Nspl,Nlead))
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Regs(ilon,ilat)>0) then
               Do ispl=1,Nspl
                  X(1:Nlead)=EC_bc(ilon,ilat,ispl,1:Nlead)
                  Call Moving_Ave (X,Nlead,Lw,Y,Miss)
                  EC_w(ilon,ilat,ispl,1:Nlead)=Y(1:Nlead)
!print"(22F6.1)",EC_w(ilon,ilat,ispl,10:30)
                  X(1:Nlead)=OBS(ilon,ilat,ispl,1:Nlead)
                  Call Moving_Ave (X,Nlead,Lw,Y,Miss)
                  Obs_w(ilon,ilat,ispl,1:Nlead)=Y(1:Nlead)
               Enddo
            endif
         Enddo
      Enddo 
   End Subroutine Calc_weekly_data
   !----------------------------------------------------------------------------------------------
      Subroutine Moving_Ave (X,N,m,Y,Miss)
         Integer :: N,m, i,j,k,m2   ! m must be an odd number
         Real :: X(N), Y(N), A(m), tmp, Miss

         Y(:) = Miss; m2=m/2
         Do j=1,N-m+1
            k = j+m-1
            A(1:m) = X(j:k)
            tmp = Ave(A,m,Miss)
            Y(j+m2)=tmp*m           ! Gan cho gia tri giua
         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"
   Include "/home/tanpv/Mylibs/Sub_neural_2025.f90"
   Include "/home/tanpv/Mylibs/Sub_julian_day.f90"
   END
   !---------------------------------------------------------------------------