PROGRAM Combine_all_BC_1982_2020   ! Combine_all_BC_1982-2020.f90 
   Use Netcdf

   Integer, Parameter :: Yr1=1982, Yr2=2011, Nmon=12, Nreg=7, NLead_BL=6, Nvar=5, Nfile=9
   Integer, Parameter :: Nyrs=(Yr2-Yr1+1), Ntime_BL=Nyrs*Nmon, NLead=NLead_BL, Nmem=7
   Integer, Parameter :: Yr1_test=2011, Yr2_test=2021, Nyrs_test=(Yr2_test-Yr1_test+1), Ntime_test=Nyrs_test*Nmon
   Integer, Parameter :: Ntime=(Yr2_test-Yr1+1)*Nmon
   Real, Parameter    :: Miss=-99.
   Character (Len=3)  :: vname(Nvar)=(/"obs","ann","mod","cli","mlr"/)

   Character (Len=255) :: Fname
   Real, Allocatable   :: lon(:),lat(:),Regs(:,:), X3d(:,:,:)

   Real, Allocatable   :: xlon(:),xlat(:),xlev(:),xtime(:)
   Integer             :: Nlo,Nla,Nz,Nt 
   Integer             :: Nlon, Nlat, Ngrd, ilon, ilat, it
   Integer             :: Mvars

   Real, Allocatable   :: Dat_BL  (:,:,:,:,:)    ! (Nlon,Nlat,Ntime_BL,Nvar,NLead_BL)
   Real, Allocatable   :: Dat_test(:,:,:,:,:,:)  ! (Nlon,Nlat,Ntime_test,Nvar,Nmem,NLead)
   Real, Allocatable   :: Dat_All (:,:,:,:,:,:)  ! (Nlon,Nlat,Ntime,Nvar,Nmem,NLead)
   Real, Allocatable   :: Dat_Reg (:,:,:,:,:)    ! (Nreg,Ntime,Nvar,Nmem,NLead)
   Real, Allocatable   :: Dat_Reg_box (:,:,:,:,:,:)

   Character (Len=5)   :: yvarn, Model_nm

   !--------------------------------------------------------------------------------
   Integer :: Nyear, imon      ! Tong so nam co SL va chi so thang
   !--------------------------------------------------------------------------------
   character (len=250) :: Dir, Dir_Cli, Dir_4boxpl, Dir_Corrected, Dir_Combined

   Integer :: Nsample, Nepoc, Nsample1
   integer :: i,j,k
   Real :: Q25(Nreg,Nmon,Nvar,Nmem,NLead), Q75(Nreg,Nmon,Nvar,Nmem,NLead)

   Real, Allocatable :: X4d(:,:,:,:), X5d(:,:,:,:,:), F1d(:),O1d(:), Score_Cor(:,:,:,:), Score_AV(:,:,:,:) 
   Namelist /Model_info/ Model_nm, yvarn

   !--------------------------------------------------------------------------------

   Open (1, file="Namelist_Comb_for_Dr.txt",Status="old")
   Read (1, Model_info)
   Close (1)

   Fname="/home/tanpv/Data/AWC/Mask_VN_0.25_Corrected.nc"
   Call rd_ncdf_2d (Fname,lon,lat,Nlon,Nlat,Regs,"reg") ! 7 sub-regs mask

   Dir_Corrected ="/home/tanpv/Data/562_VNU/0.25_VN/VAnh/Corrected/"
   Dir_Combined  ="/home/tanpv/Data/562_VNU/0.25_VN/VAnh/Combined/"

!   Dir         ="/home/tanpv/Data/562_VNU/0.25_VN/Dat_All/"
!   Dir_4Dr     ="/home/tanpv/Data/562_VNU/0.25_VN/4Drou/"

   !--------------------------------------------------------------------------------

   Call Read_all_data_baseline (Dir_Corrected,Model_nm,yvarn) 
   Call Read_all_data_test     (Dir_Corrected,Model_nm,yvarn)
   Call Combine
   Call Write_for_drought (Dir_Combined,yvarn,Model_nm)
   Call Areal_Mean        ! Sub-regional mean
   Call Write_for_drought_Reg (Dir_Combined,yvarn,Model_nm,Ntime)
   Call Subregions        ! Sub-regional gridpoint values

   Dir_4boxpl="/home/tanpv/Data/562_VNU/0.25_VN/VAnh/4Box_Plot/"
   !Call Write_for_drought_Reg (Dir,yvarn,Model_nm,Nsample)

   Call Write_for_boxplot (Dir_4boxpl,yvarn,Model_nm)

   !--------------------------------------------------------------------------------
      
   CONTAINS
   !--------------------------------------------------------------------------------
   Subroutine Areal_Mean
      Integer :: ilon,ilat,imon,imon1,ilead,ireg, imem, ic, ic0, ic1, ic2, it
      Real :: tmp
      Logical :: Valid
      ! Dat_Reg (:,:,:,:,:)    ! (Nreg,Ntime,Nvar,Nmem,NLead)

      Allocate (Dat_Reg(Nreg,Ntime,Nvar,Nmem,NLead))
      Do ilead=1,NLead
         Do imem=1,Nmem
            Do ivar=1,Nvar
               Do it=1,Ntime
                  Do ireg=1,Nreg
                     ic=0; tmp=0.0; ic0=0
                     Do ilon=1,Nlon
                        Do ilat=1,Nlat
                           if (Int(Regs(ilon,ilat))==ireg) then
                              ic0=ic0+1
                              if (Dat_All(ilon,ilat,it,ivar,imem,ilead)/=Miss) then
                                 ic=ic+1; tmp=tmp+Dat_All(ilon,ilat,it,ivar,imem,ilead)
                              endif
                           endif
                        Enddo
                     Enddo
                     if (ic>0) then 
                        Dat_Reg(ireg,it,ivar,imem,ilead)=tmp/real(ic)
                     else
                        Dat_Reg(ireg,it,ivar,imem,ilead)=Miss
                     endif
                     if (it>ilead.AND.it<469+ilead) then
                        if (ic<10) print"(A,7I5,A)","ilead,imem,ivar,it,ireg,ic0,ic: ",ilead,imem,ivar,it,ireg,ic0,ic," ???"
                     endif
                  Enddo
               Enddo
            Enddo
         Enddo
      Enddo
      ! Make sure all variables have valid values:
      Do ilead=1,NLead
         Do imem=1,Nmem
            Do ireg=1,Nreg
               Do it=1,Ntime
                  Valid=.True.
                  Do ivar=1,Nvar
                     if (Dat_Reg(ireg,it,ivar,imem,ilead)==Miss) Valid=.False.
                  Enddo
                  if (.not.Valid) Dat_Reg(ireg,it,1:Nvar,imem,ilead)=Miss
               Enddo
            Enddo
         Enddo
      Enddo
   End Subroutine Areal_Mean
   !--------------------------------------------------------------------------------
   Integer Function Get_Reg_Points (ireg)
      Integer :: ireg,ilon,ilat,Npoint

      Npoint=0
      Do ilon=1,Nlon
         Do ilat=1,Nlat
            if (Regs(ilon,ilat)==ireg) Npoint=Npoint+1
         Enddo
      Enddo
      Get_Reg_Points = Npoint
   End Function Get_Reg_Points
   !--------------------------------------------------------------------------------
   Subroutine Subregions
      Integer :: Npoint,ilon,ilat,imon,imon1,ilead,ireg, imem, ic, ic0, ic1, ic2, it
      Real :: tmp
      !Real :: Q25(Nreg,Nmon,Nvar,Nmem,NLead), Q75(Nreg,Nmon,Nvar,Nmem,NLead)
      Real, Allocatable :: X1d(:)

      ! Dat_Reg (:,:,:,:,:)    ! (Nreg,Ntime,Nvar,Nmem,NLead)

      Nsample=0
      Do ireg=1,Nreg
         Npoint=Get_Reg_Points (ireg); if (Npoint>Nsample) Nsample=Npoint
      Enddo
      Nsample=Nsample*(Ntime/Nmon)
      Allocate (X1d(Nsample))
      Allocate (Dat_Reg_box(Nreg,Nmon,Nsample,Nvar,Nmem,NLead)); Dat_Reg_box=Miss
      Do ilead=1,NLead
         Do imem=1,Nmem
            Do ivar=1,Nvar
               Do ireg=1,Nreg 
                  Do imon=1,Nmon
                     ic=0
                     Do it=1,Ntime
                        imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon
                        if (imon1==imon) then
                           Do ilon=1,Nlon ! it=1,Ntime 
                              Do ilat=1,Nlat !ilon=1,Nlon
                                 if (Int(Regs(ilon,ilat))==ireg) then
                                    if (Dat_All(ilon,ilat,it,ivar,imem,ilead)/=Miss) then
                                       ic=ic+1
                                       Dat_Reg_box(ireg,imon,ic,ivar,imem,ilead)=Dat_All(ilon,ilat,it,ivar,imem,ilead)
                                       X1d(ic)=Dat_Reg_box(ireg,imon,ic,ivar,imem,ilead)
                                    endif
                                 endif
                              Enddo
                           Enddo
                        endif 
                     Enddo
                     Q25(ireg,imon,ivar,imem,ilead)=percentile(X1d,ic,0.25,Miss)
                     Q75(ireg,imon,ivar,imem,ilead)=percentile(X1d,ic,0.75,Miss) 
                  Enddo
               Enddo
            Enddo
         Enddo
      Enddo
      Deallocate (X1d)
   End Subroutine Subregions
   !--------------------------------------------------------------------------------
   Subroutine Write_for_boxplot (Dir,yvname,Model_nm)
      Character (Len=*) :: Dir,yvname,Model_nm
      Character (Len=250) :: var_nm, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm
      Real, Allocatable, Dimension (:) :: xdim,ydim,zlev, xtime
      INteger :: Ntime, i, Nvar1   ! Dat_Reg(Nreg,Nmon,Nsample,Nvar,Nmem,NLead)

      Ntime=Nsample
      print*,"Nsample,Nreg,Nmem,Nmon=", Ntime,Nreg,Nmem,Nmon
      Allocate (xdim(Nreg),ydim(Nmem),zlev(Nmon),xtime(Ntime))

      Do i=1,Nreg; xdim(i)=i; Enddo; Do i=1,Nmem; ydim(i)=i; Enddo
      Do i=1,Nmon; zlev(i)=i; Enddo; Do i=1,Ntime; xtime(i)=i; Enddo

      t_unit="months since 1982-01-01 "; var_unit="-"
      lon_nm="reg"; lat_nm="mem"; t_nm="sample"; z_nm="mon"

      Allocate (X4d(Nreg,Nmem,Nmon,Ntime)); X4d=Miss
      Do ilead=1,NLead
         write (Fname,"(A,I2.2,A)") trim(Dir)//trim(yvname)//"_"//trim(Model_nm)//"_lead_",ilead,".nc"
         Do ivar=1,Nvar 
            Do ireg=1,Nreg 
               Do imem=1,Nmem 
                  Do imon=1,Nmon
                     Do it=1,Ntime
                        X4d(ireg,imem,imon,it)=Dat_Reg_box(ireg,imon,it,ivar,imem,ilead)
                     Enddo
                  Enddo
               Enddo
            Enddo 
            var_nm=trim(vname(ivar))
            var_ln="Monthly mean values of "//trim(var_nm)//" of "//trim(yvname)
            if (ivar==1) then
               Call wr_ncdf_4d (Fname,xdim,ydim,zlev,xtime,Nreg,Nmem,Nmon,Ntime,X4d,var_nm, &
                          t_unit, var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
            else
               Call wr_ncdf_4d_apd (Fname,X4d,var_nm,Nreg,Nmem,Nmon,Ntime, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
            endif
         Enddo
      Enddo
      Deallocate (xdim,ydim,zlev,xtime,X4d)
   End Subroutine Write_for_boxplot
   !--------------------------------------------------------------------------------
   Subroutine Combine 
      Integer :: ilon,ilat,imon,imon1,ilead,ireg, ic, ic0, ic1, ic2, it
      Real :: X1(Ntime_BL), Y1(Ntime_test), Z1(Ntime) 
     
      Allocate (Dat_All(Nlon,Nlat,Ntime,Nvar,Nmem,NLead))
      Dat_All=Miss
      Do ilead=1,NLead
         Do ivar=1,Nvar
            Do imem=1,Nmem
               Do ilon=1,Nlon
                  Do ilat=1,Nlat
                     if (Int(Regs(ilon,ilat))>0) then
                        X1(1:Ntime_BL)=Dat_BL(ilon,ilat,1:Ntime_BL,ivar,ilead)
                        ! Chu y: X1(1:Ntime_BL) nhu nhau voi tat ca cac imem
                        Y1(1:Ntime_test)=Dat_test(ilon,ilat,1:Ntime_test,ivar,imem,ilead)
                        Call Combine_Arr1 (X1,Y1,Z1)
                        Dat_All(ilon,ilat,1:Ntime,ivar,imem,ilead)=Z1(1:Ntime)
                     endif
                  Enddo
               Enddo
            Enddo
         Enddo
      Enddo
      Where (Dat_All/=Miss.AND.Dat_All<0) Dat_All=0.0
   End Subroutine Combine
   !--------------------------------------------------------------------------------
   Subroutine Combine_Arr1 (X,Y,Z)
      Integer :: i,j,k,it1,it2
      Real :: X(Ntime_BL),Y(Ntime_test),Z(Ntime)
      Real :: A(Nmon),B(Nmon),C(Nmon)
!   1/82       12/2010
!     1 2 3 ...it1-1 |1/2011          12/2021|  ...       Ntime
!     1 2 3 ...it1-1 | it1 it1+1 ... Ntime_BL|  
!                    | 1     2 ...       Nmon| Nmon+1 ... Ntime_test
 
      it1=Ntime_BL-Nmon+1; it2=Nmon
      A(1:Nmon)=X(it1:Ntime_BL); B(1:Nmon)=Y(1:it2)
      Do i=1,Nmon
         if (A(i)/=Miss.AND.B(i)/=Miss) then
            !C(i)=(A(i)+B(i))/2.0
            C(i)=B(i)
         elseif (A(i)/=Miss.AND.B(i)==Miss) then
            C(i)=A(i)
         elseif (A(i)==Miss.AND.B(i)/=Miss) then
            C(i)=B(i)
         else
            C(i)=Miss
         endif
      Enddo
!print"(20F6.1)",A(1:Nmon)
!print"(20F6.1)",B(1:Nmon)
!print"(20F6.1)",C(1:Nmon)
      Z(1:it1-1)=X(1:it1-1); Z(it1:Ntime_BL)=C(1:Nmon); Z(Ntime_BL+1:Ntime)=Y(it2+1:Ntime_test)
   End Subroutine Combine_Arr1
   !--------------------------------------------------------------------------------
   Subroutine Write_for_drought_Reg (Dir,yvname,Model_nm,Nspl)
      Character (Len=*) :: Dir,yvname,Model_nm
      Integer :: Nspl, Ntime, imem, ilead, ivar,it,i
      Character (Len=250) :: var_nm, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm
      Real, Allocatable, Dimension (:) :: xreg, ymem, zlead, xtime
      ! Dat_Reg (:,:,:,:,:)    ! (Nreg,Ntime,Nvar,Nmem,NLead)

      Ntime=Nspl
      Allocate (xreg(Nreg),ymem(Nmem),zlead(NLead),xtime(Ntime))
      Allocate (X4d(Nreg,Nmem,NLead,Ntime))
      Do i=1,Nreg; xreg(i)=i; Enddo; Do i=1,Ntime; xtime(i)=i; Enddo
      Do i=1,Nmem; ymem(i)=i; Enddo; Do i=1,NLead; zlead(i)=i; Enddo
      t_unit="months since 1982-01-01 "; var_unit="-"
      lon_nm="reg"; lat_nm="mem"; t_nm="month"; z_nm="lead"
      write (Fname,"(A)") trim(Dir)//trim(yvname)//"_all_reg_"//trim(Model_nm)//".nc"
      Do ivar=1,Nvar
         Do ireg=1,Nreg
            Do imem=1,Nmem
               Do ilead=1,NLead
                  Do it=1,Ntime
                     X4d(ireg,imem,ilead,it)=Dat_Reg(ireg,it,ivar,imem,ilead)
                  Enddo
               Enddo
            Enddo
         Enddo
         var_nm=trim(vname(ivar))
         var_ln="Regional average of Monthly mean values of "//trim(yvname)//". x:1-Nreg, y:1-Nmem, z:1-NLead, t:1-Ntime"
         if (ivar==1) then
            Call wr_ncdf_4d (Fname,xreg,ymem,zlead,xtime,Nreg,Nmem,NLead,Ntime,X4d,var_nm, &
                            t_unit, var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss)
         else
            Call wr_ncdf_4d_apd (Fname,X4d,var_nm,Nreg,Nmem,NLead,Ntime, &
                                var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         endif
      Enddo
      Deallocate (xreg,ymem,zlead,xtime,X4d)
   End Subroutine Write_for_drought_Reg
   !--------------------------------------------------------------------------------
   Subroutine Write_for_drought (Dir,yvname,Model_nm)
      Character (Len=*) :: Dir,yvname,Model_nm
      Character (Len=250) :: var_nm, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm
      Real, Allocatable, Dimension (:) :: zlev, xtime
      INteger :: i, Nvar1

      Allocate (zlev(Nmem),xtime(Ntime))
      Do i=1,Nmem; zlev(i)=i; Enddo; Do i=1,Ntime; xtime(i)=i; Enddo
      t_unit="months since 1982-01-01 "; var_unit="-"
      lon_nm="lon"; lat_nm="lat"; t_nm="month"; z_nm="mem"

      Allocate (X4d(Nlon,Nlat,Nmem,Ntime))
      Do ilead=1,NLead
         write (Fname,"(A,I2.2,A)") trim(Dir)//trim(yvname)//"_"//trim(Model_nm)//"_lead_",ilead,".nc"
         Do ivar=1,Nvar
            Do imem=1,Nmem
               Do it=1,Ntime
                  X4d(:,:,imem,it)=Dat_All(:,:,it,ivar,imem,ilead)
               Enddo
            Enddo
            var_ln="Monthly mean values of "//trim(yvname)//". z=1..6 member, z=7 Ens. t=1..360: 1-1982..12-2011, t=361..480: 1-2012..12-2021" 
            var_nm=trim(vname(ivar))
            if (ivar==1) then
               Call wr_ncdf_4d (Fname,Lon,Lat,zlev,xtime,Nlon,Nlat,Nmem,Ntime,X4d,var_nm, &
                            t_unit, var_ln, "-", lon_nm, lat_nm, z_nm, t_nm, Miss)
            else
               Call wr_ncdf_4d_apd (Fname,X4d,var_nm,Nlon,Nlat,Nmem,Ntime, &
                                var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) 
            endif
         Enddo
      Enddo
      Deallocate (zlev,xtime,X4d)
   End Subroutine Write_for_drought
   !--------------------------------------------------------------------------------
   Subroutine Read_all_data_baseline (Dir,Model_nm,yvarn)
      Character (Len=*) :: Dir,Model_nm,yvarn
      Character (Len=255) :: Fname
      Real, Allocatable :: xlon(:),xlat(:),xlev(:),xtime(:),X4d(:,:,:,:)
      Integer :: it,ilead,ivar,Nlo,Nla,Nz,Nt

      !print*,"Read_all_data_baseline: NLead_BL,Ntime=",NLead_BL,Ntime; pause
      ! Note: Baseline data timeseries include 1982-2021 (480 months)
      Allocate (Dat_BL(Nlon,Nlat,Ntime,Nvar,NLead_BL)) ! (Nlon,Nlat,Ntime,Nvar,NLead))
      Dat_BL=Miss
      Do ilead=1,NLead_BL
         write (Fname,"(A,I2.2,A)") trim(Dir)//"BC_Training_"//trim(yvarn)//"_"//trim(Model_nm)//"_lead_",ilead,".nc" 
         Call rd_ncdf_4d (Fname,xlon,xlat,xlev,xtime,Nlo,Nla,Nz,Nt,X4d,trim(yvarn))
         if (Nt/=Ntime_BL) then; print*,"Something's wrong in Read_all_data_baseline"; stop; endif
         Do ivar=1,Nvar  ! =5: "obs","ann","mod","cli","mlr"
            Do it=1,Ntime_BL  ! Ntime
               Dat_BL(:,:,it,ivar,ilead)=X4d(:,:,ivar,it)
            Enddo
         Enddo
         Deallocate (xlon,xlat,xlev,xtime,X4d)
      Enddo
      Where (Dat_BL/=Miss.AND.Dat_BL<0) Dat_BL=0.0
   End Subroutine Read_all_data_baseline 
   !--------------------------------------------------------------------------------
   Subroutine Read_all_data_test (Dir,Model_nm,yvarn)
      Character (Len=*) :: Dir,Model_nm,yvarn
      Character (Len=255) :: Fname
      Real, Allocatable :: xlon(:),xlat(:),xlev(:),xtime(:),X4d(:,:,:,:)
      Integer :: iz,it,ilead,ivar,Nlo,Nla,Nz,Nt

      !print*,"Read_all_data_test: Ntime_test",Ntime_test; pause
      Allocate (Dat_test(Nlon,Nlat,Ntime_test,Nvar,Nmem,NLead))
      Dat_test=Miss
      Do ilead=1,NLead
         write (Fname,"(A,I2.2,A)") trim(Dir)//"BC_testing_"//trim(yvarn)//"_"//trim(Model_nm)//"_lead_",ilead,".nc"
         Do ivar=1,Nvar
            Call rd_ncdf_4d (Fname,xlon,xlat,xlev,xtime,Nlo,Nla,Nz,Nt,X4d,trim(vname(ivar)))
            if (Nt/=Ntime_test) then; print*,"Something's wrong in Read_all_data"; Stop; endif
            Do iz=1,Nz
               Do it=1,Ntime_test
                  Dat_test(:,:,it,ivar,iz,ilead)=X4d(:,:,iz,it)
               Enddo
            Enddo
            Deallocate (xlon,xlat,xlev,xtime,X4d)
         Enddo
      Enddo
      Where (Dat_test/=Miss.AND.Dat_test<0) Dat_test=0.0
   End Subroutine Read_all_data_test
   !--------------------------------------------------------------------------------
     
   !Include "/home/tanpv/Mylibs/Sub_HQTT.f90"
   !Include "/home/tanpv/Mylibs/Sub_LAD.f90"
   Include "/home/tanpv/Mylibs/Percentile.f90"
   !Include "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90"  ! Correlation (X, Y, N, Corr, Prob, Miss)
   !Include "/home/tanpv/Mylibs/Sub_Verification.f90"              ! ME, MAE, RME, RMAE
   Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90" 
   Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
   !Include "/home/tanpv/Mylibs/Sub_neural_2023.f90"
   Include "/home/tanpv/Mylibs/Stdev.f90"
   Include "/home/tanpv/Mylibs/Ave.f90"
   !--------------------------------------------------------------------------------
   END
   !Include "/home/tanpv/Mylibs/Student_Prob.f90"
   !-------------------------------------------------------------------------------------------