Use NetCDF
     Integer, parameter :: Nseason=4, Nreg=6, Nexp=43   ! Num of experiments
     Integer, parameter :: Nmon=12, Nday=31, Yr1=2011, Yr2=2015, Nyear=Yr2-Yr1+1
     Integer, parameter :: Nper=8, Ntime=Nyear*365
     Integer, parameter :: NPCs=15
     Integer, parameter :: MonSeason(3,Nseason)=(/12,1,2,3,4,5,6,7,8,9,10,11/)
     Integer, parameter :: EXP_ID(Nexp)=(/ 1, 2, 3, 4, 5, 6, 7, 8, 9,10, &  ! 10
                                          11,   13,14,15,   17,18,19,20, &  ! 8 (11)
                                          21,22,23,   25,26,27,28,29,30, &  ! 9 (23)
                                          31,32,33,34,35,   37,38,39,40, &  ! 9 (35)
                                          41,42,43,44,45,46,47/)            ! 7 (47)
     Real, parameter :: xexp(1:Nexp)=EXP_ID(1:Nexp)
     Integer, parameter :: Nbin=9 ! 16 ! 20
     Real, parameter :: Miss=-99.0
     !Real, Parameter :: Bin(0:Nbin)=(/0.,0.1,0.2,0.3,0.5,0.7,1.0,1.5,2.0, &
     !                                2.5,3.0,3.5,4.0,4.5,5.0,7.0,10.0,20.0, &
     !                               30.0,50.0,99999.0/)
     !Real, Parameter :: xbin(Nbin)=(/0.05,0.15,0.25,0.40,0.60,0.85,1.25,1.75, &
     !                                2.25,2.75,3.25,3.75,4.25,4.75,6.00,8.50, &
     !                                15.0,25.0,40.0,70.0/)

!     Real, Parameter :: Bin(0:Nbin)=(/0.,1.,2.,3.,5.,7.,10.,15.,17., &
!                                     20.,25.,30.,35.,40.,45.,50.,99999.0/)
!     Real, Parameter :: xbin(Nbin)=(/0.5,1.5,2.5,4.0,6.0,8.5,12.5,16.0,17.5, &
!                                     22.5,27.5,32.5,37.5,42.5,47.5,70.0/)

     Real, Parameter :: Bin(0:Nbin)=(/1.,5.,10., &
                                     15.,20.,25.,30.,40.,50.,99999.0/)
     Real, Parameter :: xbin(Nbin)=(/3.0,7.5,12.5,17.5,22.5,27.5,35.,45.,60./)

     Real, Parameter :: p(Nper)=(/10.,20.,40.,50.,60.,80.,90.,95./)

     Real, Parameter :: xreg(Nreg)=(/1.,2.,3.,4.,5.,6./)
     Real, Parameter :: xmon(Nmon)=(/1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12./)
     Real, Parameter :: xseason(Nseason)=(/1.,2.,3.,4./)  ! DJF,MAM,JJA,SON
     Real, Allocatable :: OBS(:,:,:,:,:), EXPs(:,:,:,:,:,:)     
     Real, Allocatable :: Ob3d(:,:,:), EXP3d(:,:,:,:)       ! (nlon,nlat,Ntime, Nexp)
     Real, Allocatable :: Lon(:), Lat(:), Mask(:,:)
     Real, Allocatable :: Mask_reg(:,:)            ! Mask of subreg
     Integer :: Nlon, Nlat !, Ntime
     !--------------------------------------------------------------------------------------
     Real :: FrqO(Nreg,Nseason,Nbin)        ! (Nreg,Nseason,Nbin)
     Real :: FrqE(Nreg,Nseason,Nexp,Nbin)   ! (Nreg,Nseason,Nbin,Nexp)
     Real :: MeanO(Nreg,Nmon)               ! (Nreg,Nmon)
     Real :: MeanE(Nreg,Nmon,Nexp) !,Nmon)          ! (Nreg,Nexp,Nmon)
     Real :: IntO(Nreg,Nmon)                ! (Nreg,Nmon) wet-day
     Real :: IntE(Nreg,Nmon,Nexp) !,Nmon)           ! (Nreg,Nexp,Nmon) wet-day
     Real :: PerO(Nreg,Nper)                ! Percentiles
     Real :: PerE(Nreg,Nper,Nexp) !,Nper)
     Integer :: Npoint(Nreg)                ! Num of gridpoint in each subreg
     !--------------------------------------------------------------------------------------
     Real :: Cor_Frq (Nreg,Nexp,Nseason),P_Frq(Nreg,Nexp,Nseason)  ! Nbin
     Real :: Cor_Mean(Nreg,Nexp), P_Mean(Nreg,Nexp)    ! Nmon
     Real :: Cor_Int (Nreg,Nexp), P_Int (Nreg,Nexp)    ! Nmon
     Real :: Cor_Per (Nreg,Nexp), P_Per (Nreg,Nexp)    ! Nper
     Real :: Cor_PCs (NPCs,Nexp), P_PCs (NPCs,Nexp)    ! Landpoints
     Real :: Cor_EOF (NPCs,Nexp), P_EOF (NPCs,Nexp)    ! Landpoints
     !--------------------------------------------------------------------------------------
     Real :: ME_Frq (Nreg,Nexp,Nseason), MAE_Frq (Nreg,Nexp,Nseason)
     Real :: ME_Mean(Nreg,Nexp), MAE_Mean(Nreg,Nexp)
     Real :: ME_Int (Nreg,Nexp), MAE_Int (Nreg,Nexp)
     Real :: ME_Per (Nreg,Nexp), MAE_Per (Nreg,Nexp)
     
     !--------------------------------------------------------------------------------------
     ! PCA for OBS:
      Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:)
      Real*8, Allocatable :: EOFs(:,:), PCs(:,:)
      Real*8, Allocatable   :: Eigenvalues(:), ECumulate(:)
      Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K
      Integer :: Nx, Nt, Nxx
      Real, Allocatable :: PC_Obs(:,:),  EOF_Obs(:,:)
      Real, Allocatable :: PC_Exp(:,:,:),EOF_Exp(:,:,:)
     !--------------------------------------------------------------------------------------
     ! Myanmar:  92 101  15 29
     ! Thailand: 97 102   7 21
     ! Malay:    95 107 -10  7
     ! Indo:    107 145 -11  7
     ! Phi:     115 130   5 20
     ! VN:      102 110   8 24
      Real, Parameter :: Lon1(Nreg)=(/ 92., 97., 95.,107.,115.,102./)
      Real, Parameter :: Lon2(Nreg)=(/101.,102.,107.,145.,130.,110./)
      Real, Parameter :: Lat1(Nreg)=(/ 15.,  7.,-10.,-11.,  5.,  8./)
      Real, Parameter :: Lat2(Nreg)=(/ 29., 21.,  7.,  7., 20., 24./)

     !--------------------------------------------------------------------------------------
     Character (Len=255) :: Fname
     Character (Len=200) :: Fnm(Nexp) !=(/"Data/pr_exp7_daily.nc", &
                                      !  "Data/pr_exp28_daily.nc"/)
     Integer :: Ndys,ilon,ilat,it,iyr,imon,iday,i,j
     Logical :: OK
     !--------------------------------------------------------------------------------------
     Do i=1,Nexp
        write(Fnm(i),"(A,I2.2,A)") "Data/Exp_Era_2021/pr_exp",EXP_ID(i),"_daily.nc"
        INQUIRE(FILE=trim(Fnm(i)), EXIST=OK)
        !print*,i,trim(Fnm(i)), OK
     Enddo
     !stop
     Call Read_Mask ("Data/RegCM_ST_mask.nc")
     Call Read_TRMM_GPCC    ! Read OBS 
     Call Calc_Npoint
     Call Calc_FrqO !("NC/Freq_O.nc")
     ! Print test:
     !-------------------------------------------------
     Allocate (EXPs(Nlon,Nlat,Nyear,Nmon,Nday,Nexp))
     Allocate (EXP3d(Nlon,Nlat,Ntime,Nexp))
     EXPs=Miss; EXP3d=Miss
     Do j=1,Nexp
        Call Read_Exp (Fnm(j),"pr",j)
     Enddo
     Call Calc_FrqExp ("NC/Freq.nc")
     Call Calc_Mean ("NC/Mean.nc") !,"NC/Mean_E.nc")
     Call Calc_percentiles ("NC/Percentiles.nc") !,"NC/Per_E.nc")
     !-------------------------------------------------
     Call OBS_for_PCA
     !-------------------------------------------------
     Allocate (PC_Obs(Nxx,NPCs),  EOF_Obs(Nt,NPCs))         ! For correlation
     Allocate (PC_Exp(Nxx,NPCs,Nexp),EOF_Exp(Nt,NPCs,Nexp))
     PC_Obs=Miss; EOF_Obs=Miss; PC_Exp=Miss; EOF_Exp=Miss
     !-------------------------------------------------
     switch=1; Standard=1; PC_Stard=0; Nmode=NPCs
      if (switch==1) then
         Mcount=Nt; Ncount=Nx
      else if (switch==0) then
         Mcount=Nx; Ncount=Nt
      endif
      Allocate (ydat(Nt,Nx),Xtb(Nx), Sx(Nx))   ! TB, DLC, ma tran Covar hoac Corre
      Allocate (Eigenvalues(Mcount))          ! Tri rieng va Vecto rieng (1D)
      Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode))
      Factor = Sqrt(Real(Nt))
      print*,"Check Nx, Nt : ",Nx,Nt
      Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, &
                        Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode)
      print*,"Writing PCs..."
      Call Write_EOF_PCA ("NC/pca.nc","NC/eof.nc",0)
      PC_Obs(1:Nx,1:Nmode)=PCs(1:Nx,1:Nmode)
      EOF_Obs(1:Nt,1:Nmode)=EOFs(1:Nt,1:Nmode)
      Deallocate (Xdat,Ydat,Xtb,Sx,Eigenvalues,EOFs,PCs,ECumulate)
      Do j=1,Nexp
         Call EXP_for_PCA (j)
         Allocate (ydat(Nt,Nx),Xtb(Nx), Sx(Nx))
         Allocate (Eigenvalues(Mcount))
         Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode))
         Factor = Sqrt(Real(Nt))
         print*,"Check Exp, Nx, Nt : ",j,Nx,Nt
         Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, &
                        Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode)
         print*,"Writing PCs..."
         Call Write_EOF_PCA ("NC/pca.nc","NC/eof.nc",j)
         PC_Exp(1:Nx,1:Nmode,j)=PCs(1:Nx,1:Nmode)
         EOF_Exp(1:Nt,1:Nmode,j)=EOFs(1:Nt,1:Nmode)
         Deallocate (Xdat,Ydat,Xtb,Sx,Eigenvalues,EOFs,PCs,ECumulate)
      Enddo
      Call Calc_Errors ("NC/ME_MAE_chars.nc")
      Call Calc_Corr ("NC/Cor_chars.nc","NC/Cor_pcs.nc")
     !--------------------------------------------------------------------------------------
     CONTAINS
     !--------------------------------------------------------------------------------------
     Subroutine Calc_Npoint
        Integer :: ilon,ilat,ireg
        Logical :: LLon,LLat

        Allocate (Mask_reg(Nlon,Nlat)); Mask_reg(:,:)=Miss
        Do ireg=1,Nreg
           Npoint(ireg)=0
           Do ilon=1,Nlon
              LLon=Lon(ilon)>=Lon1(ireg).AND.Lon(ilon)<=Lon2(ireg)
              Do ilat=1,Nlat
                 LLat=Lat(ilat)>=Lat1(ireg).AND.Lat(ilat)<=Lat2(ireg)
                 if (Mask(ilon,ilat)>0.AND.LLon.AND.LLat) then
                    Npoint(ireg)=Npoint(ireg)+1
                    Mask_reg(ilon,ilat)=ireg
                 endif
              Enddo
           Enddo
        Enddo
        Do ireg=1,Nreg
           Npoint(ireg)=0
           Do ilon=1,Nlon
              Do ilat=1,Nlat
                 if (int(Mask_reg(ilon,ilat))==ireg) Npoint(ireg)=Npoint(ireg)+1
              Enddo
           Enddo
        Enddo
        print*,Npoint(:)
        ! Print for testing:
        !Call wr_ncdf_2d ("Mask_reg.nc",lon,lat,nlon,nlat,Mask_reg,"land", &
        !             "mask subreg", "-", "lon","lat", Miss)
     End Subroutine Calc_Npoint
     !--------------------------------------------------------------------------------------
     Subroutine Calc_Errors (Fname)
        Character (Len=*) :: Fname
        Integer :: i,j,k, ireg,isea,iexp,imon,iper
        Real, Allocatable :: X1d(:),Y1d(:),XX1d(:),YY1d(:),xexp(:),xmode(:)
        Real :: me1,me2,me3,me4,mae1,mae2,mae3,mae4
        Real :: n11,n12,n13,n14,n21,n22,n23,n24
!     Real :: ME_Frq (Nreg,Nexp,Nseason), MAE_Frq (Nreg,Nexp,Nseason)
!     Real :: ME_Mean(Nreg,Nexp), MAE_Mean(Nreg,Nexp)
!     Real :: ME_Int (Nreg,Nexp), MAE_Int (Nreg,Nexp)
!     Real :: ME_Per (Nreg,Nexp), MAE_Per (Nreg,Nexp)

        ME_Frq=Miss; MAE_Frq=Miss; ME_Mean=Miss; MAE_Mean=Miss
        ME_Int=Miss; MAE_Int=Miss; ME_Per=Miss; MAE_Per=Miss
        Allocate (X1d(Nbin),Y1d(Nbin))
        Do ireg=1,Nreg
           Do isea=1,Nseason
              X1d(1:Nbin)=FrqO(ireg,isea,1:Nbin)
              Do iexp=1,Nexp
                 Y1d(1:Nbin)=FrqE(ireg,isea,iexp,1:Nbin)
                 Call ME_MAE(X1d,Y1d,Nbin,ME1,MAE1,Miss) 
                 ME_Frq(ireg,iexp,isea)=ME1; MAE_Frq(ireg,iexp,isea)=MAE1 
              Enddo
           Enddo
        Enddo
        Deallocate (X1d,Y1d)
        ! Mean and Intensity:
        Allocate (X1d(Nmon),Y1d(Nmon),XX1d(Nmon),YY1d(Nmon))
        Do ireg=1,Nreg
           X1d(1:Nmon)=MeanO(ireg,1:Nmon); XX1d(1:Nmon)=IntO(ireg,1:Nmon)
           Do iexp=1,Nexp
              Y1d (1:Nmon)=MeanE(ireg,1:Nmon,iexp)
              YY1d(1:Nmon)=IntE (ireg,1:Nmon,iexp)
              Call ME_MAE(X1d,Y1d,Nmon,ME1,MAE1,Miss)
              ME_Mean(ireg,iexp)=ME1; MAE_Mean(ireg,iexp)=MAE1
              Call ME_MAE(XX1d,YY1d,Nmon,ME1,MAE1,Miss)
              ME_Int (ireg,iexp)=ME1; MAE_Int (ireg,iexp)=MAE1
           Enddo
        Enddo
        Deallocate (X1d,Y1d,XX1d,YY1d)
        ! Percentiles:
        Allocate (X1d(Nper),Y1d(Nper))
        Do ireg=1,Nreg
           X1d(1:Nper)=PerO(ireg,1:Nper)
           Do iexp=1,Nexp
              Y1d(1:Nper)=PerE(ireg,1:Nper,iexp)
              Call ME_MAE(X1d,Y1d,Nper,ME1,MAE1,Miss)
              ME_Per(ireg,iexp)=ME1; MAE_Per(ireg,iexp)=MAE1
           Enddo
        Enddo
        Deallocate (X1d,Y1d)
        ! 
        print*,"Writing ME, MAE ..."
        Allocate (xexp(Nexp))
        Do iexp=1,Nexp; xexp(iexp)=EXP_ID(iexp); Enddo
        Call wr_ncdf_3d (Fname,xreg,xexp,xseason,Nreg,Nexp,Nseason, &
             ME_Frq,"me_frq","months since 2011-01-01 ", &
           "Mean Errors between frq_E and frq_O", &
           "%", "reg","exps","season", Miss)
        Call wr_ncdf_3d_apd (Fname,MAE_Frq,"mae_frq",Nreg,Nexp,Nseason, &
            "MAE frq_E-frq_O", "%", "reg", "exps", "season", Miss)
        Call wr_ncdf_2d_apd (Fname,ME_Mean,"me_mean",Nreg,Nexp, &
           "ME between mean_O-mean_E", "mm/day", "reg","exps",Miss)
        Call wr_ncdf_2d_apd (Fname,MAE_Mean,"mae_mean",Nreg,Nexp, &
           "MAE mean_O-mean_E", "mm/day", "reg","exps",Miss)
        Call wr_ncdf_2d_apd (Fname,ME_Int,"me_int",Nreg,Nexp, &
           "ME between Int_O-Int_E", "mm/day", "reg","exps",Miss)
        Call wr_ncdf_2d_apd (Fname,MAE_Int,"mae_int",Nreg,Nexp, &
           "MAE Int_O-Int_E", "mm/day", "reg","exps",Miss)
        Call wr_ncdf_2d_apd (Fname,ME_Per,"me_per",Nreg,Nexp, &
           "ME between Per_O-Per_E", "mm/day", "reg","exps",Miss)
        Call wr_ncdf_2d_apd (Fname,MAE_Per,"mae_per",Nreg,Nexp, &
           "MAE Per_O-Per_E", "mm/day", "reg","exps",Miss)
     End Subroutine Calc_Errors
     !--------------------------------------------------------------------------------------
     Subroutine ME_MAE(X,Y,N,ME,MAE,Miss)
        Integer :: N, i,j
        Real :: X(N),Y(N),ME,MAE,Miss
        Real :: tmp1,tmp2,nc

        tmp1=0; tmp2=0; nc=0
        Do i=1,N
           if (X(i)/=Miss.AND.Y(i)/=Miss) then
              nc=nc+1
              tmp1=tmp1+Y(i)-X(i); tmp2=tmp2+ABS(Y(i)-X(i))
           endif
        Enddo
        ME=tmp1/nc; MAE=tmp2/nc
     End Subroutine ME_MAE
     !--------------------------------------------------------------------------------------
     Subroutine Calc_Corr (Fname,Fname1)
        Character (Len=*) :: Fname,Fname1
        Integer :: i,j,k, ireg,isea,iexp,imon,iper
        Real, Allocatable :: X1d(:),Y1d(:),XX1d(:),YY1d(:),xexp(:),xmode(:)
        Real :: Corr,prob

        Cor_Frq=Miss; P_Frq=Miss; Cor_Mean=Miss; P_Mean=Miss
        Cor_Int=Miss; P_Int=Miss; Cor_Per=Miss;  P_Per=Miss
        ! Frq:
        Allocate (X1d(Nbin),Y1d(Nbin))
        Do ireg=1,Nreg
           Do isea=1,Nseason
              X1d(1:Nbin)=FrqO(ireg,isea,1:Nbin)
              Do iexp=1,Nexp
                 Y1d(1:Nbin)=FrqE(ireg,isea,iexp,1:Nbin)
                 Call Correlation(X1d,Y1d,Nbin,Corr,prob,Miss)
                 Cor_Frq(ireg,iexp,isea)=Corr; P_Frq(ireg,iexp,isea)=prob 
              Enddo
           Enddo
        Enddo
        Deallocate (X1d,Y1d)
        ! Mean and Intensity:
        Allocate (X1d(Nmon),Y1d(Nmon),XX1d(Nmon),YY1d(Nmon))
        Do ireg=1,Nreg
           X1d(1:Nmon)=MeanO(ireg,1:Nmon); XX1d(1:Nmon)=IntO(ireg,1:Nmon)
           Do iexp=1,Nexp
              Y1d (1:Nmon)=MeanE(ireg,1:Nmon,iexp)
              YY1d(1:Nmon)=IntE (ireg,1:Nmon,iexp)
              Call Correlation(X1d,Y1d,Nmon,Corr,prob,Miss)
              Cor_Mean(ireg,iexp)=Corr; P_Mean(ireg,iexp)=prob
              Call Correlation(XX1d,YY1d,Nmon,Corr,prob,Miss)
              Cor_Int (ireg,iexp)=Corr; P_Int (ireg,iexp)=prob
           Enddo
        Enddo
        Deallocate (X1d,Y1d,XX1d,YY1d)
        ! Percentiles:
        Allocate (X1d(Nper),Y1d(Nper))
        Do ireg=1,Nreg
           X1d(1:Nper)=PerO(ireg,1:Nper)
           Do iexp=1,Nexp
              Y1d(1:Nper)=PerE(ireg,1:Nper,iexp)
              Call Correlation(X1d,Y1d,Nper,Corr,prob,Miss)
              Cor_Per(ireg,iexp)=Corr; P_Per(ireg,iexp)=prob
           Enddo
        Enddo
        Deallocate (X1d,Y1d)
        ! 
        print*,"Writing Correlations 1"
        Allocate (xexp(Nexp))
        Do iexp=1,Nexp; xexp(iexp)=EXP_ID(iexp); Enddo
        Call wr_ncdf_3d (Fname,xreg,xexp,xseason,Nreg,Nexp,Nseason, &
             Cor_Frq,"cor_frq","months since 2011-01-01 ", &
           "Correlations between frq_E and frq_O", &
           "-", "reg","exps","season", Miss)
        Call wr_ncdf_3d_apd (Fname,P_Frq,"p_frq",Nreg,Nexp,Nseason, &
            "p-value of cor frq_E-frq_O", "-", "reg", "exps", "season", Miss)
        Call wr_ncdf_2d_apd (Fname,Cor_Mean,"cor_mean",Nreg,Nexp, &
           "Correlations between mean_O-mean_E", "-", "reg","exps",Miss)
        Call wr_ncdf_2d_apd (Fname,P_Mean,"p_mean",Nreg,Nexp, &
           "p-value of Corr mean_O-mean_E", "-", "reg","exps",Miss)
        Call wr_ncdf_2d_apd (Fname,Cor_Int,"cor_int",Nreg,Nexp, &
           "Correlations between Int_O-Int_E", "-", "reg","exps",Miss)
        Call wr_ncdf_2d_apd (Fname,P_Int,"p_int",Nreg,Nexp, &
           "p-value of Corr Int_O-Int_E", "-", "reg","exps",Miss)
        Call wr_ncdf_2d_apd (Fname,Cor_Per,"cor_per",Nreg,Nexp, &
           "Correlations between Per_O-Per_E", "-", "reg","exps",Miss)
        Call wr_ncdf_2d_apd (Fname,P_Per,"p_per",Nreg,Nexp, &
           "p-value of Corr Per_O-Per_E", "-", "reg","exps",Miss)
        !
     !Allocate (PC_Obs(Nxx,NPCs),  EOF_Obs(Nt,NPCs))         ! For correlation
     !Allocate (PC_Exp(Nxx,NPCs,Nexp),EOF_Exp(Nt,NPCs,Nexp))
     !Real :: Cor_PCs (NPCs,Nexp), P_PCs (NPCs,Nexp)    ! Landpoints
     !Real :: Cor_EOF (NPCs,Nexp), P_EOF (NPCs,Nexp)
        Cor_PCs=Miss; P_PCs=Miss; Cor_EOF=Miss; P_EOF=Miss
        Allocate (X1d(Nxx),Y1d(Nxx),XX1d(Nt),YY1d(Nt))
        Do i=1,Nmode
           X1d(1:Nxx)=PC_Obs(1:Nxx,i); XX1d(1:Nt)=EOF_Obs(1:Nt,i)
           Do iexp=1,Nexp
              Y1d(1:Nxx)=PC_Exp(1:Nxx,i,iexp); 
              YY1d(1:Nt)=EOF_Exp(1:Nt,i,iexp)
              Call Correlation(X1d,Y1d,Nxx,Corr,prob,Miss)
              Cor_PCs(i,iexp)=Corr; P_PCs(i,iexp)=prob
              Call Correlation(XX1d,YY1d,Nt,Corr,prob,Miss)
              Cor_EOF(i,iexp)=Corr; P_EOF(i,iexp)=prob
           Enddo   
        Enddo
        !
        print*,"Writing Correlations 2"
        Allocate (xmode(Nmode))
        Do i=1, Nmode; xmode(i)=i; Enddo
        Call wr_ncdf_2d (Fname1,xmode,xexp,Nmode,Nexp,Cor_PCs,"cor_pc", &
             "Corr between Obs and Exp PCs", "-", "mode","exps", Miss)
        Call wr_ncdf_2d_apd (Fname1,P_PCs,"p_pc",Nmode,Nexp, &
             "p-value Corr Obs and Exp PCs", "-", "mode","exps", Miss)

        Call wr_ncdf_2d_apd (Fname1,Cor_EOF,"cor_eof",Nmode,Nexp, &
             "Corr between Obs and Exp EOF", "-", "mode","exps", Miss)
        Call wr_ncdf_2d_apd (Fname1,P_EOF,"p_eof",Nmode,Nexp, &
             "p-value Corr Obs and Exp EOF", "-", "mode","exps", Miss)
        Deallocate (X1d,Y1d,XX1d,YY1d,xmode,xexp)
     End Subroutine Calc_Corr
     !--------------------------------------------------------------------------------------
     Subroutine Calc_percentiles (Fname2) !,Fname2)
        Character (Len=*) :: Fname2 !,Fname2
        Integer :: N,ilon,ilat,ireg,ip, iyr,i, it
        Real, Allocatable :: X(:),xexp(:)

        Do ireg=1,Nreg
           N=Npoint(ireg)*Ntime
           Allocate (X(N)); i=0
           Do ilon=1,Nlon
              Do ilat=1,Nlat
                 if (Int(Mask_reg(ilon,ilat))==ireg) then
                    Do it=1,Ntime
                       i=i+1; X(i)=Ob3d(ilon,ilat,it)
                    Enddo
                 endif
              Enddo
           Enddo
           print*,"Check in Calc_percentiles: N,i=", N,i
           Do ip=1,Nper
              PerO(ireg,ip)=percentile(X,N,p(ip)/100.,Miss) ! percentile(X,N,p(ip)/100.,Miss)
              print*,p(ip)
           Enddo
           ! Experiments:
           Do iexp=1,Nexp
              i=0
              Do ilon=1,Nlon
                 Do ilat=1,Nlat
                    if (Int(Mask_reg(ilon,ilat))==ireg) then
                       Do it=1,Ntime
                          i=i+1; X(i)=EXP3d(ilon,ilat,it,iexp)
                       Enddo
                    endif
                 Enddo
              Enddo
              print*,"Check in Calc_percentiles: N,i=", N,i
              Do ip=1,Nper
                 PerE(ireg,ip,iexp)=percentile(X,N,p(ip)/100.,Miss) ! percentile(X,N,p(ip)/100.,Miss)
                 print*,p(ip)
              Enddo
           Enddo
           Deallocate (X)
        Enddo
        !Call wr_ncdf_2d (Fname1,xreg,p,Nreg,Nper,PerO,"per", &
        !   "Daily prec. precentiles: 10.,20.,40.,50.,60.,80.,90.,95. %", &
        !   "mm/day", "reg","percentile", Miss)
        !Call wr_ncdf_2d_apd (Fname1,IntO,"inten",Nreg,Nmon, &
        !         "Intensity of daily prec.", "mm/day", "reg", "mon", Miss)
        Allocate (xexp(Nexp))
        Do iexp=1,Nexp; xexp(iexp)=EXP_ID(iexp); Enddo
        Call wr_ncdf_3d (Fname2,xreg,p,xexp,Nreg,Nper,Nexp,PerE,"per_e", &
           "months since 2011-01-01 ", &
           "Daily prec. precentiles: 10.,20.,40.,50.,60.,80.,90.,95. %", &
           "mm/day", "reg","percentile","exps", Miss)
        Call wr_ncdf_2d_apd (Fname2,PerO,"per_o",Nreg,Nper, &
           "Daily prec. precentiles: 10.,20.,40.,50.,60.,80.,90.,95. %", &
           "mm/day", "reg","percentile",Miss)  
        !Call wr_ncdf_3d_apd (Fname2,IntE,"inten",Nreg,Nexp,Nmon, &
        !    "Intensity of daily prec.", "mm/day", "reg", "exps", "mon", Miss)
     End Subroutine Calc_percentiles
     !--------------------------------------------------------------------------------------
     Subroutine Calc_Mean (Fname2) !,Fname2)
        Character (Len=*) :: Fname2 !,Fname2
        Integer :: ilon,ilat,ireg,imon,iyr
        Integer :: N,Ndys,Ndy,isea,i,j,iexp,N1
        Real, Allocatable :: X(:), Y(:), xexp(:), X3d(:,:,:)
        ! Real :: MeanO(Nreg,Nmon) Real :: MeanE(Nreg,Nexp,Nmon) 

        Do ireg=1,Nreg
           Do imon=1,Nmon
              Ndy=DayOfMon_NoLeap(imon)
              Allocate (X(Npoint(ireg)*Nyear*Ndy))
              Allocate (Y(Npoint(ireg)*Nyear*Ndy))
              N=0; N1=0
              Do ilon=1,Nlon
                 Do ilat=1,Nlat
                    if (Int(Mask_reg(ilon,ilat))==ireg) then
                       Do iyr=1,Nyear
                          Do iday=1,Ndy
                             if (OBS(ilon,ilat,iyr,imon,iday)/=Miss) then
                                N=N+1; X(N)=OBS(ilon,ilat,iyr,imon,iday)
                             endif
                             if (OBS(ilon,ilat,iyr,imon,iday)>=1.0) then   ! wet-day
                                N1=N1+1; Y(N1)=OBS(ilon,ilat,iyr,imon,iday)
                             endif
                          Enddo
                       Enddo
                    endif
                 Enddo
              Enddo 
              MeanO(ireg,imon)=Ave(X,N,Miss)
              IntO (ireg,imon)=Ave(Y,N1,Miss)
              !Deallocate (X)
              Do iexp=1,Nexp
                 N=0; N1=0
                 Do ilon=1,Nlon
                    Do ilat=1,Nlat
                       if (Int(Mask_reg(ilon,ilat))==ireg) then
                          Do iyr=1,Nyear
                             Do iday=1,Ndy
                                if (EXPs(ilon,ilat,iyr,imon,iday,iexp)/=Miss) then
                                   N=N+1; X(N)=EXPs(ilon,ilat,iyr,imon,iday,iexp)
                                endif
                                if (EXPs(ilon,ilat,iyr,imon,iday,iexp)>=1.0) then   ! Wet-day
                                   N1=N1+1; Y(N1)=EXPs(ilon,ilat,iyr,imon,iday,iexp)
                                endif
                             Enddo
                          Enddo
                       endif
                    Enddo
                 Enddo
                 !MeanE(ireg,iexp,imon)=Ave(X,N,Miss)
                 !IntE (ireg,iexp,imon)=Ave(Y,N1,Miss)
                 MeanE(ireg,imon,iexp)=Ave(X,N,Miss)
                 IntE (ireg,imon,iexp)=Ave(Y,N1,Miss)
              Enddo
              Deallocate (X,Y)
           Enddo
        Enddo
        !Call wr_ncdf_2d (Fname1,xreg,xmon,Nreg,Nmon,MeanO,"mean", &
        !         "Monthly mean of daily prec.", "mm/day", "reg","mon", Miss)
        !Call wr_ncdf_2d_apd (Fname1,IntO,"inten",Nreg,Nmon, &
        !         "Intensity of daily prec.", "mm/day", "reg", "mon", Miss)
        Allocate (xexp(Nexp))
        Do iexp=1,Nexp; xexp(iexp)=EXP_ID(iexp); Enddo
        Call wr_ncdf_3d (Fname2,xreg,xmon,xexp,Nreg,Nmon,Nexp,MeanE,"mean_e", &
                   "months since 2011-01-01 ","Monthly mean of daily prec.", &
                   "mm/day", "reg", "mon","exps", Miss)
        Call wr_ncdf_3d_apd (Fname2,IntE,"inten_e",Nreg,Nmon,Nexp, &
            "Intensity of daily prec.", "mm/day", "reg", "mon", "exps", Miss)
        Call wr_ncdf_2d_apd (Fname2,MeanO,"mean_o",Nreg,Nmon, &
             "Monthly mean of daily prec.", "mm/day", "reg","mon", Miss)
        Call wr_ncdf_2d_apd (Fname2,IntO,"inten_o",Nreg,Nmon, &
             "Intensity of daily prec.", "mm/day", "reg", "mon", Miss)
     End Subroutine Calc_Mean
     !--------------------------------------------------------------------------------------
     Subroutine Calc_FrqExp (Fname)
        Character (Len=*) :: Fname
        Integer :: ilon,ilat,ireg,imon,iyr,ibin
        Integer :: N,Ndys,Ndy,isea,i,j,iexp
        Real, Allocatable :: X(:), xexp(:), X4d(:,:,:,:)
        Real :: a,b

        Do iexp=1,Nexp
           Do ireg=1,Nreg
              Do isea=1,Nseason
                 Ndys=0
                 Do i=1,3
                    imon=MonSeason(i,isea) !,MonSeason(3,isea)
                    Ndys=Ndys+DayOfMon_NoLeap(imon);
                 Enddo
                 Allocate (X(Npoint(ireg)*Nyear*Ndys))
                 N=0
                 Do ilon=1,Nlon
                    Do ilat=1,Nlat
                       if (Int(Mask_reg(ilon,ilat))==ireg) then
                          Do iyr=1,Nyear
                             Do j=1,3
                                imon=MonSeason(j,isea) !,MonSeason(3,isea)
                                Ndy=DayOfMon_NoLeap(imon)
                                Do iday=1,Ndy
                                   if (EXPs(ilon,ilat,iyr,imon,iday,iexp)/=Miss) then
                                      N=N+1; X(N)=EXPs(ilon,ilat,iyr,imon,iday,iexp)
                                   endif
                                Enddo
                             Enddo
                          Enddo
                       endif
                    Enddo
                 Enddo
                 Do ibin=1,Nbin
                    a=Bin(ibin-1); b=Bin(ibin)
                    FrqE(ireg,isea,iexp,ibin)=Freq(X,N,a,b)
                 Enddo
                 Deallocate (X)
              Enddo
           Enddo
        Enddo
        Allocate (xexp(Nexp))
        Do iexp=1,Nexp; xexp(iexp)=EXP_ID(iexp); Enddo
        Call wr_ncdf_4d (Fname,xreg,xseason,xexp,xbin,Nreg, &
              Nseason,Nexp,Nbin,FrqE,"frq_e", "days since 2011-01-01 ", &
              "Daily Prec Freq","-", "reg", "season","exp","bin", Miss) 
        Allocate (X4d(Nreg,Nseason,1,Nbin))
        Do ireg=1,Nreg
           Do isea=1,Nseason
              X4d(ireg,isea,1,:)=FrqO(ireg,isea,:)
           Enddo
        Enddo
        Call wr_ncdf_4d_apd (Fname,X4d,"frq_o",Nreg,Nseason,1,Nbin, &
             "Daily Prec Freq","-","reg", "season","exp","bin", Miss)

     End Subroutine Calc_FrqExp
     !--------------------------------------------------------------------------------------
     Subroutine Calc_FrqO !(Fname) 
        !Character (Len=*) :: Fname
        Integer :: ilon,ilat,ireg,imon,iyr,ibin
        Integer :: N,Ndys,Ndy,isea,i,j
        Real, Allocatable :: X(:)
        Real :: a,b

        Do ireg=1,Nreg
           Do isea=1,Nseason
              Ndys=0
              Do i=1,3
                 imon=MonSeason(i,isea) !,MonSeason(3,isea)
                 Ndys=Ndys+DayOfMon_NoLeap(imon);
              Enddo 
              Allocate (X(Npoint(ireg)*Nyear*Ndys))
              N=0
              Do ilon=1,Nlon
                 Do ilat=1,Nlat
                    if (Int(Mask_reg(ilon,ilat))==ireg) then
                       Do iyr=1,Nyear
                          Do j=1,3
                             imon=MonSeason(j,isea) !,MonSeason(3,isea)
                             Ndy=DayOfMon_NoLeap(imon)
                             Do iday=1,Ndy
                                if (OBS(ilon,ilat,iyr,imon,iday)/=Miss) then
                                   N=N+1; X(N)=OBS(ilon,ilat,iyr,imon,iday)
                                endif
                             Enddo
                          Enddo
                       Enddo
                    endif
                 Enddo
              Enddo
              Do ibin=1,Nbin
                 a=Bin(ibin-1); b=Bin(ibin)
                 FrqO(ireg,isea,ibin)=Freq(X,N,a,b)
              Enddo
              Deallocate (X)
           Enddo
        Enddo 
        !Call wr_ncdf_3d (Fname,xreg,xseason,xbin,Nreg,Nseason,Nbin, FrqO, &
        !      "frq_o", "days since 2011-01-01 ","Daily Prec Freq", &
        !       "-", "reg", "season","bin", Miss)
     End Subroutine Calc_FrqO 
     !--------------------------------------------------------------------------------------
     Real Function Freq(X,N,a,b)
        Integer :: N,i,j,ic
        Real :: X(N),a,b,tmp

        ic=0
        Do i=1,N
           if (X(i)>=a.AND.X(i)<b) ic=ic+1
        Enddo
        Freq=real(ic)/real(N)*100.0
     End Function Freq
     !--------------------------------------------------------------------------------------
     Subroutine EXP_for_PCA (iexp)
        Integer :: iexp, ilon,ilat,iyr,imon,iday,it,i,j,k,Ndys,Ndy
        Integer :: Landpoints,Npentad,it1,Nxx
        Real, Allocatable :: X1d(:),Y1d(:),YY(:),xxdat(:,:)
        Real :: tmp

        Ndys=365; Landpoints=0;
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              if (Mask(ilon,ilat)/=Miss) then
                 Landpoints=Landpoints+1;
              endif
           Enddo
        Enddo

        Nxx=Landpoints; Npentad=Ndys/5; Nt=Npentad ! *Nyear
        Nx=Nxx
        Allocate (xxdat(Nt,Nxx))
        Allocate (X1d(Ntime),Y1d(Nt*Nyear))
        Allocate (YY(Npentad))
        i=0
        Do ilat=1,Nlat
           Do ilon=1,Nlon
              if (Mask(ilon,ilat)/=Miss) then
                 X1d(1:Ntime)=EXP3d(ilon,ilat,1:Ntime,iexp)
                 Call Calc_pentad (X1d,Ntime,Y1d)
                 !Call Anomaly (Y1d,Nt,Npentad)
                 !i=i+1;
                 Call Longterm_mean (Y1d,Nt*Nyear,YY,Npentad)
                 tmp=Ave(YY,Nt,Miss)
                 if (tmp/=0.0.AND.tmp/=Miss) then
                    i=i+1
                    xxdat(1:Nt,i)=YY(1:Nt)
                 else
                    Nx=Nx-1; Mask(ilon,ilat)=Miss
                 endif
                 !i=i+1; xdat(1:Nt,i)=Y1d(1:Nt)
              endif
           Enddo
        Enddo
        Allocate (xdat(Nt,Nx))
        Do i=1,Nx; xdat(1:Nt,i)=xxdat(1:Nt,i); Enddo
        print*,"Check Landpoints, Nx, Nt,Ntime : ",Nxx,Nx,Nt,Ntime
        !if (iexp==7.or.iexp==17.or.iexp==28.or.iexp==39) then
        !   print*,xdat; pause
        !endif
        Deallocate (X1d,Y1d,YY)
     End Subroutine EXP_for_PCA
     !--------------------------------------------------------------------------------------
     Subroutine OBS_for_PCA
        Integer :: ilon,ilat,iyr,imon,iday,it,i,j,k,Ndys,Ndy
        Integer :: Landpoints,Npentad,it1
        !Integer, Parameter :: Np=20
        Real, Allocatable :: X3d(:,:,:),X2d(:,:),X1d(:),Y1d(:), Yp(:),YY(:)
        !Real, Parameter :: p(Np)=(/5.,10.,15.,20.,25.,30.,35.,40.,45.,50., &
        !                            55.,60.,65.,70.,75.,80.,85.,90.,95.,99./)
        Real :: tmp, ntmp

        Ndys=365; Landpoints=0;   ! Ntime=total days in Nyear
        Nxx=0
        Allocate (X3d(Nlon,Nlat,Ndys),X2d(Nyear,Ndys),X1d(Nyear))
        Allocate (Y1d(Ndys))
        X3d=Miss; 
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              if (Mask(ilon,ilat)/=Miss) then
                 Landpoints=Landpoints+1; it1=0
                 Nxx=Nxx+1
                 Do iyr=1,Nyear
                    it=0
                    Do imon=1,Nmon
                       Ndy=DayOfMon_NoLeap(imon)
                       Do iday=1,Ndy
                          it=it+1; X2d(iyr,it)=OBS(ilon,ilat,iyr,imon,iday)
                          !it1=it1+1; Ob3d(ilon,ilat,it1)=OBS(ilon,ilat,iyr,imon,iday)
                       Enddo
                    Enddo
                 Enddo
                 Do it=1,Ndys
                    X1d(1:Nyear)=X2d(1:Nyear,it); 
                    X3d(ilon,ilat,it)=Ave(X1d,Nyear,Miss)  ! 5-year daily average
                 Enddo
                 Y1d(1:Ndys)=X3d(ilon,ilat,1:Ndys)
                 if (Ar1_Empty(Y1d,Ndys, Miss)) then
                    Mask(ilon,ilat)=Miss; X3d(ilon,ilat,:)=Miss
                    Landpoints=Landpoints-1
                 endif
              endif
           Enddo
        Enddo

        ! PCA applied to daily mean:
        ! Nx=Landpoints; Nt=Ndys
        !Allocate (xdat(Nt,Nx)) 
        !i=0
        !Do ilat=1,Nlat
        !   Do ilon=1,Nlon
        !      if (Mask(ilon,ilat)/=Miss) then
        !         i=i+1
        !         xdat(1:Nt,i)=X3d(ilon,ilat,1:Nt)
        !      endif
        !   Enddo
        !Enddo

        print*,"Check Landpoints, Nlon*Nlat : ",Landpoints, Nlon*Nlat
        ! PCA applied to pentad:

        Nx=Landpoints; Npentad=Ndys/5; Nt=Npentad ! *Nyear
        Allocate (xdat(Nt,Nx))
        Deallocate (X1d,Y1d); Allocate (X1d(Ntime),Y1d(Nt*Nyear))
        Allocate (YY(Npentad))
        i=0
        Do ilat=1,Nlat
           Do ilon=1,Nlon
              if (Mask(ilon,ilat)/=Miss) then
                 X1d(1:Ntime)=Ob3d(ilon,ilat,1:Ntime)
                 Call Calc_pentad (X1d,Ntime,Y1d)
                 !Call Anomaly (Y1d,Nt,Npentad)
                 i=i+1; 
                 Call Longterm_mean (Y1d,Nt*Nyear,YY,Npentad)
                 xdat(1:Nt,i)=YY(1:Nt)
                 !i=i+1; xdat(1:Nt,i)=Y1d(1:Nt)
              endif
           Enddo
        Enddo
        print*,"Check Landpoints, i, Nt,Ntime : ",Landpoints, i,Nt,Ntime
        Deallocate (X1d,Y1d,YY)
     End Subroutine OBS_for_PCA
     !--------------------------------------------------------------------------------------
     Subroutine Longterm_mean (X,N,Y,Ncyle)
        Integer :: N,Nyear,Ncyle,i,j,k
        Real :: X(N), Y(Ncyle)

        Nyear=N/Ncyle
        Y(:)=0.0
        if (N/=Ncyle*Nyear) then; print*,"Check Longterm_mean: ",N,Nyear,Ncyle; stop; endif
        Do i=1,Ncyle
           Do j=1,N
              k=mod(j,Ncyle); if (k==0) k=Ncyle
              if (k==i) Y(i)=Y(i)+X(j)
           Enddo
           Y(i)=Y(i)/real(Nyear)
        Enddo
     End Subroutine Longterm_mean
     !--------------------------------------------------------------------------------------
     Subroutine Calc_pentad (X,N,Y)
        Integer :: N, Npentad, i,j
        Real :: X(N),Y(N/5), tmp

        Npentad=N/5
        if (N/=Npentad*5) then; print*,"Check pentad N,Npentad= ",N,Npentad; stop; endif
        Do i=1,Npentad
           tmp=0
           Do j=(i-1)*5+1,i*5
              tmp=tmp+X(j)
           Enddo
           Y(i)=tmp
        Enddo 
     End Subroutine Calc_pentad
     !--------------------------------------------------------------------------------------
     Subroutine Read_Exp (Fname,vname,ivar)
        Character (Len=*) :: Fname,vname
        Integer :: ilon,ilat,it,i,j,k,ivar
        Integer :: Nlo, Nla, Nt
        Real, Allocatable :: X(:,:,:),xlon(:),xlat(:),xtime(:)
        Real, Allocatable :: X1d(:)

        Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X,trim(vname))
        !if (ivar==7.or.ivar==17.or.ivar==28.or.ivar==39) then
        !    print*,"Check Read_Exp: ivar, Nlo,Nla,Nt=",ivar,Nlo,Nla,Nt;  pause
        !    print*,X; pause
        !endif
        Where (IsNaN(X)) X=Miss
        Where (X/=Miss) X=X*86400.0   ! Convert mm/s to mm/day
        !Allocate (EXP3d(Nlon,Nlat,Ntime,Nexp)
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              if (Mask(ilon,ilat)/=Miss) then
                 it=365; i=0                      ! Bo nam 2010
                 Do iyr=Yr1,Yr2
                    Do imon=1,Nmon
                       Ndys=DayOfMon(imon,iyr)
                       Do iday=1,Ndys
                          it=it+1; EXPs(ilon,ilat,iyr-Yr1+1,imon,iday,ivar)=X(ilon,ilat,it)
                          !i=i+1; EXP3d(ilon,ilat,i,ivar)=X(ilon,ilat,it)
                       Enddo
                    Enddo
                 Enddo
                ! print*,"Check: Nt,it, i: ",Nt,it-365,i; pause
              endif
           Enddo
        Enddo
        it=0
        Do iyr=1,Nyear
           Do imon=1,Nmon
              Ndy=DayOfMon_NoLeap(imon)
              Do i=1,Ndy
                 it=it+1; EXP3d(:,:,it,ivar)=EXPs(:,:,iyr,imon,i,ivar)
              Enddo
           Enddo
        Enddo
        !if (ivar==7.or.ivar==17.or.ivar==28.or.ivar==39) then
        !    print*,"Check Read_Exp: ivar, Nlo,Nla,Nt=",ivar,Nlo,Nla,Nt;  pause
            !print*,EXP3d; pause
        !endif

        Deallocate (xlon,xlat,xtime,X)
     End Subroutine Read_Exp
     !--------------------------------------------------------------------------------------
     Subroutine Read_TRMM_GPCC
        Real, Parameter :: Missing_TRMM=-9999.9, Missing_GPCC=-9999.0
        Real :: Missing_value
        Integer :: ilon,ilat,it,i,j,k,iyr,imon
        Integer :: Nlo, Nla, Nt, Ndy
        Real, Allocatable :: X(:,:,:),xlon(:),xlat(:),xtime(:)
        Real, Allocatable :: X1d(:)


        !Fname="Data/TRMM_2011_2015.nc";          Missing_value=Missing_TRMM
        !Fname="Data/pr_0.25_GPCC_2011_2015.nc";  Missing_value=Missing_GPCC
        Fname="Data/TRMM.2011.2015.nc";          Missing_value=Missing_TRMM

        Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X,"pr")
        !print*,"Checking in read Prec: Nlo,Nla,Nt=",Nlo,Nla,Nt; pause
        Where (IsNaN(X)) X=Miss
        Where (X==Missing_value) X=Miss

        Allocate (OBS(Nlon,Nlat,Nyear,Nmon,Nday)); OBS=Miss
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              !if (Mask(ilon,ilat)/=Miss) then
                 it=0
                 Do iyr=Yr1,Yr2
                    Do imon=1,Nmon
                       Ndys=DayOfMon(imon,iyr)
                       Do iday=1,Ndys
                          it=it+1; OBS(ilon,ilat,iyr-Yr1+1,imon,iday)=X(ilat,ilon,it) !X(ilon,ilat,it)
                       Enddo
                    Enddo
                 Enddo
                 !print*,"Check: Nt,it: ",Nt,it; pause
              !endif
           Enddo
        Enddo
        !Ntime=Nyear*365
        Allocate (Ob3d(Nlon,Nlat,Ntime)) ! Daily data, no leap
        it=0
        Do iyr=1,Nyear
           Do imon=1,Nmon
              Ndy=DayOfMon_NoLeap(imon)
              Do i=1,Ndy
                 it=it+1; Ob3d(:,:,it)=OBS(:,:,iyr,imon,i)
              Enddo
           Enddo
        Enddo

        Deallocate (X,xlon,xlat,xtime)

     End Subroutine Read_TRMM_GPCC
     !--------------------------------------------------------------------------------------
     Subroutine Read_APH
        ! Return: APH,Nlon,Nlat,Lon,Lat,Mask
        Real, Parameter :: Miss_APH=-99.9
        Integer :: ilon,ilat,it,i,j,k
        Integer :: Nlo, Nla, Nt
        Real, Allocatable :: X(:,:,:),xlon(:),xlat(:),xtime(:)
        Real, Allocatable :: X1d(:)

        Fname="/work/projects/G09/SEACLID/APH_2021/tmp/All.nc"  ! APH
        Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X,"precip")
        Where (IsNaN(X)) X=Miss
        Where (X==Miss_APH) X=Miss
        Nlon=Nlo; Nlat=Nla
        Allocate (Mask(Nlon,Nlat),X1d(Nt),Lon(Nlon),Lat(Nlat)); Mask=Miss
        Lon(:)=xlon(:); Lat(:)=xlat(:)
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              X1d(1:Nt)=X(ilon,ilat,1:Nt)
              if (.not.Ar1_Empty(X1d,Nt, Miss)) Mask(ilon,ilat)=1
           Enddo
        Enddo
        Allocate (OBS(Nlon,Nlat,Nyear,Nmon,Nday)); OBS=Miss
        Do ilon=1,Nlon
           Do ilat=1,Nlat
              if (Mask(ilon,ilat)/=Miss) then
                 it=0
                 Do iyr=Yr1,Yr2
                    Do imon=1,Nmon
                       Ndys=DayOfMon(imon,iyr)
                       Do iday=1,Ndys
                          it=it+1; OBS(ilon,ilat,iyr-Yr1+1,imon,iday)=X(ilon,ilat,it)
                       Enddo
                    Enddo
                 Enddo
                 !print*,"Check: Nt,it: ",Nt,it; pause
              endif
           Enddo
        Enddo
        Deallocate (X,X1d,xlon,xlat,xtime)
     End Subroutine Read_APH
     !--------------------------------------------------------------------------------------
     Subroutine Read_Mask (Fname)
        Character (Len=*) :: Fname

        Call rd_ncdf_2d (Fname,lon,lat,nlon,nlat,Mask,"mask")
        Where (IsNaN(Mask)) Mask=Miss
        Where (Mask==0) Mask=Miss
        Where (Mask>0)  Mask=1
     End Subroutine Read_Mask
      !------------------------------------------------------------------------------
      Subroutine Anomaly (X,Ntime,NCycle)
      ! Ntime: Do dai chuoi
      ! NCycle: Do dai chu ky (Nmon, NdayOfYear, etc)
      Integer :: Ntime, NCycle
      Real :: X(Ntime), Y(Ntime), TB(NCycle)
      Integer :: ic, it, imon, imon1
      Real :: tmp

      if (mod(Ntime,NCycle)/=0) print*,"Something wrong !"
      Do imon=1,NCycle
         tmp=0; ic=0
         Do it=1,Ntime
            imon1=mod(it,NCycle); if (imon1==0) imon1=NCycle
            if (imon1==imon) then
               if (X(it)/=Miss) then
                  ic=ic+1; tmp=tmp+X(it)
               endif
            endif
         Enddo
         if (ic>0) then; TB(imon)=tmp/ic; else; TB(imon)=Miss; endif
      Enddo
      Do it=1,Ntime
         imon=mod(it,NCycle); if (imon==0) imon=NCycle
         Y(it)=X(it)-TB(imon)
      Enddo
      X(:) = Y(:)
      End Subroutine Anomaly
     !--------------------------------------------------------------------------------------
      Logical Function Ar1_Empty(X,N, Miss)
        Integer :: N, i
        Real :: X(N), Miss
        Logical :: Empty
        Empty=.false.
        Do i=1,N
           if (X(i)==Miss) then
              Empty=.true.; exit
           endif
        Enddo
        Ar1_Empty=Empty
      End Function Ar1_Empty
      !-------------------------------------------------------------------------         
      Subroutine Write_EOF_PCA (F_pca, F_eof, ivar)
         Character (Len=*) :: F_eof, F_pca
         Character (Len=255) :: var_ln, t_unit, IDX_name, Vname
         Integer :: i,j,k,iidx,it, ivar
         Real :: tmp
         Real, Allocatable :: PCA(:,:,:), xmode(:), EOF1(:,:), X1(:), Xtmp(:)
         Real, Allocatable :: xtime(:)

         Allocate (PCA(Nlon,Nlat,Nmode), xmode(Nmode), Xtmp(Nmode)) !, EOF1(Nmode,Nt))
         Allocate (EOF1(Nt,Nmode))
         !-------------------------------------------------------------------------         
         print*,"Checking in write PCA: Nlon,Nlat, Nx=",Nlon,Nlat, Nx
         !print*,lon
         !print*,lat; pause
         ! Write to PCs file:
         if (ivar==0) then    ! OBS
            Vname="obs"
         else                 ! Exps
            write (Vname,"(A,I2.2)") "exp",ivar
         endif

         PCA(:,:,:)=Miss; k=0
         Do j=1, Nlat
            Do i=1, Nlon
               if (Mask(i,j)>0) then
                  k=k+1; PCA(i,j,1:Nmode) = PCs(k,1:Nmode)
               endif
            Enddo
         Enddo
         Do i=1, Nmode; xmode(i)=i; Enddo         ! PCA's modes
         var_ln="Spatial distribution of PCAs of Daily Prec"
         write (t_unit,"(A,I4.4,A)") "days since ",Yr1,"-01-01 "
         if (ivar==0) then    ! OBS
            !Vname="obs"
            Call wr_ncdf_3d (F_pca,lon,lat,xmode,nlon,nlat,Nmode, PCA, &
              trim(Vname), trim(t_unit),trim(var_ln), &
               "-", "lon", "lat","mode", Miss)
         else                 ! Exps
            !write (Vname,"(A,I2.2)") "exp",ivar
            Call wr_ncdf_3d_apd (F_pca,PCA,trim(Vname),nlon,nlat, &
                 Nmode, trim(var_ln), "-", "lon", "lat", "mode", Miss)
         endif
         !-------------------------------------------------------------------------
         Xtmp(1:Nmode) = Eigenvalues(1:Nmode)
         Call wr_ncdf_1d_apd (F_pca, Xtmp, "eig_"//trim(Vname), Nmode, &
                "Eigen values of Prec", "-", "mode", Miss)

         tmp=0.0
         Do i=1,Mcount
            tmp=tmp+Eigenvalues(i)
            if (I<=Nmode) then
               ECumulate(i)=0
               Do j=1,I
                  ECumulate(i)=ECumulate(i)+Eigenvalues(j)
               Enddo
            endif
         Enddo
         Do i=1,Nmode
            ECumulate(i)=ECumulate(i)/tmp; Xtmp(i)=Eigenvalues(i)/tmp
         Enddo

         Call wr_ncdf_1d_apd (F_pca, Xtmp, "var_"//trim(Vname),Nmode, &
                "Percent Contribution of Prec", &
                "-", "mode", Miss )

         Xtmp(1:Nmode) = ECumulate(1:Nmode)
         Call wr_ncdf_1d_apd (F_pca, Xtmp, "cum_"//trim(Vname),Nmode, &
                "Cumuluted Eigenvalues of Prec", &
                "-", "mode", Miss )
         !-------------------------------------------------------------------------
         ! Write to EOF file:
         if (Allocated(xtime)) Deallocate (xtime); Allocate (xtime(1:Nt))
         Do it=1, Nt; xtime(it)=it; Enddo     ! Months
         Do i=1, Nmode; EOF1(1:Nt,i)=EOFs(1:Nt,i); Enddo
         var_ln="Timeseries of eofs of Prec"
         if (ivar==0) then
            Call wr_ncdf_2d (F_eof,xtime,xmode, Nt,Nmode, EOF1,trim(Vname), &
                     trim(var_ln), "-", "mon","mode", Miss)
         else
            Call wr_ncdf_2d_apd (F_eof,EOF1, trim(Vname), Nt,Nmode, &
                     trim(var_ln), "-", "mon", "mode", Miss)
         endif
         print*,"Write_EOF_PCA: Done for case: ",iidx
      End Subroutine Write_EOF_PCA
      !--------------------------------------------------------------------
      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
      !--------------------------------------------------------------------
     Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
     Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
     Include "/home/tanpv/Mylibs/DayOfMon.f90"
     Include "/home/tanpv/Mylibs/DayOfMon_NoLeap.f90"
     Include "/home/tanpv/Mylibs/Ave.f90"
     Include "/home/tanpv/Mylibs/Sub_PCA_2019.f90"
     !Include "/home/tanpv/Mylibs/Percentile.f90"
     END

     Include "/home/tanpv/Mylibs/Sub_Trend_MK_Linear.f90"