! 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
   Integer             :: Yr1, Yr2, Nyrs, Nspl, Lw, Nws
   Integer, Allocatable :: time(:)

   Character (Len=255) :: Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2
   Character (Len=255) :: DirInEC, DirInObs, DirInStn 
   Real, Allocatable :: lon(:), lat(:), Regs(:,:)
   Integer :: Nlon, Nlat, Nlead
   !----------------------------------------------------------------------------------------------
   Real, Allocatable :: Tx(:,:,:), T2(:,:,:)
   !----------------------------------------------------------------------------------------------
   Real, Allocatable :: xlon(:), xlat(:), xtime(:)
   Real, Allocatable :: X4d(:,:,:,:), X3d(:,:,:) 
   Integer :: 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 :: Nfcs               ! Number of forecast over all years
   !----------------------------------------------------------------------------------------------
   Integer :: imem, ifcs, iy, im, id, NLd, iobs
   !--------------------------------------------------------------------------------
   ! ANN:
   character (len=250) :: input_file, output_file, save_file
   character (len=250) :: Dir_home, Dir_Mask, 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 :: yfc,dfc,mfc
   !-------------------------------------------------------------------------------
   yfc=2025; mfc=04; dfc=11
   Call Read_data
   Fname="Tx.nc"
   Call Write_BC_Outputs (Fname,yfc,mfc,dfc)
   CONTAINS
   !----------------------------------------------------------------------------------------------
   Subroutine Write_BC_Outputs (Fname,yfc,mfc,dfc)
      Character (Len=*) :: Fname
      Character (Len=255) :: t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, var_nm
      Real, Allocatable :: lead(:),sampl(:)
      Integer :: ilon,ilat,imon,ilead,imem,ifcs,im,iobs, i,j,iyr, sys1,sys2
      Integer :: yfc,dfc,mfc   ! forecast date and month - hdate(Nfcs) - Yr1=2001, Yr2=2020
      Integer :: days_diff, jd_ref, jd_target

      !t_unit="days since 1980-01-01 "; 
      write (t_unit,"(A,I4.4,A,I2.2,A,I2.2)") "days since ",yfc,"-",mfc,"-",dfc
      var_unit="C"; lon_nm="lon"; lat_nm="lat"; z_nm="fcs_time"; t_nm="day"

      if (Allocated (lead)) Deallocate (lead); if (Allocated (sampl)) Deallocate (sampl)
      Allocate (lead(Nlead),sampl(Nspl))
      Do ilead=1,Nlead; lead (ilead)=ilead-1; Enddo

      !Do i=1,Nlead
      !   Do j=1,Nfcs
      !      jd_ref = julian_day(1980, 1, 1) 
      !      jd_target = julian_day(yfc, mfc, dfc)
      !      sampl(i)=jd_target - jd_ref !iyr*10000+mfc(j)*100+dfc(j)
      !   Enddo
      !Enddo

      var_ln="Areal mean of Tx"
      Call wr_ncdf_3d (Fname,Lon,Lat,lead,Nlon,Nlat,Nlead,Tx,"tx",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

!      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="Moving weekly sum with a 1-day step of S2S corrected 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)
!      if (.not.training_mode) then
!         var_ln="Moving weekly sum with a 1-day step of S2S un-corrected precipitation of ecmwf forecast"
!         Call wr_ncdf_4d_apd (Fname,EC_Org_w,"org", Nlon,Nlat,Nspl,Nlead, &
!                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
!      endif
   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 :: ncid, varid, xtype, len, attnum
      Real, Allocatable :: xlon(:), xlat(:), xtime(:), X3d(:,:,:)
      Integer :: Nlo,Nla,NLt,i,ifc,ispl,iyr
      Logical :: status
      Character (Len=100) :: t_unit
      Character (Len=1) :: Ch

      ! Real, Allocatable, Dimension(:,:,:,:) :: LH, SH, OLR, TX, TM
      ! Read precipitation forecast data:
      if (Allocated (TX)) Deallocate (TX); if (Allocated (T2)) Deallocate (T2)

      Fname="/home/tanpv/Data/From_Nam/SFC_out/clwrf/2025041100/t2.nc"
      Call rd_ncdf_3d (Fname,Lon,Lat,xtime,Nlon,Nlat,Nlead,X3d,"t2")
      Deallocate (xtime)
      Fname="tx1.nc"
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlon,Nlat,Nlead,Tx,"T2MAX")
      Where (Tx==ec_miss) Tx = Miss
   End Subroutine Read_data
   !----------------------------------------------------------------------------------------------
   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_julian_day.f90"
   END
   !---------------------------------------------------------------------------