Use NetCDF
      Real, Parameter :: Miss=-99.
      Integer, Parameter :: Nmon=12, Nreg=11
      Integer, Parameter :: NPeriod=3, NSce=2, Nmodels=7   ! 7 = 6 models + Ens
      !---------------------------------------------------------------------------------
      Real, allocatable :: xlon(:), xlat(:), xtime(:)
      Integer :: Nlon, Nlat, Ntime, Nyear
      Logical :: Status

      Real, allocatable :: Mask(:,:), Regs(:,:)                   ! (Nlon,Nlat)

      ! Input Data (Reading from files): 
      Real, allocatable :: Ob_T2m(:,:,:), Ob_Pr(:,:,:)               ! (Nlon,Nlat,Ntime) - Aug 2020
      Real, allocatable :: BL_T2m(:,:,:,:), BL_Pr(:,:,:,:)           ! (Nlon,Nlat,Ntime,Nmodels)
      Real, allocatable :: MF_T2m(:,:,:,:,:), MF_Pr(:,:,:,:,:)       ! (Nlon,Nlat,Ntime,Nmodels,NSce)
      Real, allocatable :: FF_T2m(:,:,:,:,:), FF_Pr(:,:,:,:,:)       ! (Nlon,Nlat,Ntime,Nmodels,NSce)

      Real, allocatable :: T2m(:,:,:), Pr(:,:,:)
      !---------------------------------------------------------------------------------
      ! Characteristics:
      ! Truong Trung binh thang
      Real, allocatable :: Ob_T2m_Mean(:,:,:), Ob_Pr_Mean(:,:,:)            ! (Nlon,Nlat,Ntime) - Aug 2020
      Real, allocatable :: BL_T2m_Mean(:,:,:,:), BL_Pr_Mean(:,:,:,:)        !(Nlon,Nlat,Nmon,Nmodels)
      Real, allocatable :: MF_T2m_Mean(:,:,:,:,:), MF_Pr_Mean(:,:,:,:,:)    !(Nlon,Nlat,Nmon,Nmodels,NSce)
      Real, allocatable :: FF_T2m_Mean(:,:,:,:,:), FF_Pr_Mean(:,:,:,:,:)    !(Nlon,Nlat,Nmon,Nmodels,NSce)  
      ! Truong Trung binh nam
      Real, allocatable :: Ob_T2m_Ann(:,:), Ob_Pr_Ann(:,:)            !(Nlon,Nlat) - Aug 2020
      Real, allocatable :: BL_T2m_Ann(:,:,:), BL_Pr_Ann(:,:,:)        !(Nlon,Nlat,Nmodels)
      Real, allocatable :: MF_T2m_Ann(:,:,:,:), MF_Pr_Ann(:,:,:,:)    !(Nlon,Nlat,Nmodels,NSce)
      Real, allocatable :: FF_T2m_Ann(:,:,:,:), FF_Pr_Ann(:,:,:,:)    !(Nlon,Nlat,Nmodels,NSce)        
      ! Trung binh thang cho tung vung KH
      Real, allocatable :: Ob_T2m_RMon(:,:), Ob_Pr_RMon(:,:)          !(Nlon,Nlat) - Aug 2020
      Real, allocatable :: BL_T2m_RMon(:,:,:), BL_Pr_RMon(:,:,:)        !(Nreg,Nmon,Nmodels)
      Real, allocatable :: MF_T2m_Rmon(:,:,:,:), MF_Pr_RMon(:,:,:,:)    !(Nreg,Nmon,Nmodels,NSce)
      Real, allocatable :: FF_T2m_Rmon(:,:,:,:), FF_Pr_RMon(:,:,:,:)    !(Nreg,Nmon,Nmodels,NSce)  
      ! Trung binh nam cho tung vung KH
      Real, allocatable :: Ob_T2m_RAnn(:), Ob_Pr_RAnn(:)            !(Nreg) - Aug 2020
      Real, allocatable :: BL_T2m_RAnn(:,:), BL_Pr_RAnn(:,:)        !(Nreg,Nmodels)
      Real, allocatable :: MF_T2m_RAnn(:,:,:), MF_Pr_RAnn(:,:,:)    !(Nreg,Nmodels,NSce)
      Real, allocatable :: FF_T2m_RAnn(:,:,:), FF_Pr_RAnn(:,:,:)    !(Nreg,Nmodels,NSce)        

      ! Changes of Annual mean:
      Real, allocatable :: D_BL_T2m_Ann(:,:,:), D_BL_Pr_Ann(:,:,:)        !(Nlon,Nlat,Nmodels): BL-Obs - Aug 2020
      Real, allocatable :: D_MF_T2m_Ann(:,:,:,:), D_MF_Pr_Ann(:,:,:,:)    !(Nlon,Nlat,Nmodels,NSce)
      Real, allocatable :: D_FF_T2m_Ann(:,:,:,:), D_FF_Pr_Ann(:,:,:,:)    !(Nlon,Nlat,Nmodels,NSce)        
      ! P-value for Changes of Annual mean:
      Real, allocatable :: P_BL_T2m_Ann(:,:,:), P_BL_Pr_Ann(:,:,:)        !(Nlon,Nlat,Nmodels): BL-Obs - Aug 2020
      Real, allocatable :: P_MF_T2m_Ann(:,:,:,:), P_MF_Pr_Ann(:,:,:,:)    !(Nlon,Nlat,Nmodels,NSce)
      Real, allocatable :: P_FF_T2m_Ann(:,:,:,:), P_FF_Pr_Ann(:,:,:,:)    !(Nlon,Nlat,Nmodels,NSce)        
      
      ! Changes of Monthly mean:
      Real, allocatable :: D_BL_T2m_Mon(:,:,:,:), D_BL_Pr_Mon(:,:,:,:)        !(Nlon,Nlat,Nmon,Nmodels): BL-Obs - Aug 2020
      Real, allocatable :: D_MF_T2m_Mon(:,:,:,:,:), D_MF_Pr_Mon(:,:,:,:,:)    !(Nlon,Nlat,Nmon,Nmodels,NSce)
      Real, allocatable :: D_FF_T2m_Mon(:,:,:,:,:), D_FF_Pr_Mon(:,:,:,:,:)    !(Nlon,Nlat,Nmon,Nmodels,NSce)  
      ! P-value for Changes of Monthly mean:
      Real, allocatable :: P_BL_T2m_Mon(:,:,:,:), P_BL_Pr_Mon(:,:,:,:)        !(Nlon,Nlat,Nmon,Nmodels): BL-Obs - Aug 2020
      Real, allocatable :: P_MF_T2m_Mon(:,:,:,:,:), P_MF_Pr_Mon(:,:,:,:,:)    !(Nlon,Nlat,Nmon,Nmodels,NSce)
      Real, allocatable :: P_FF_T2m_Mon(:,:,:,:,:), P_FF_Pr_Mon(:,:,:,:,:)    !(Nlon,Nlat,Nmon,Nmodels,NSce)  


      Real :: DelT(Nmon+1,Nreg,NPeriod-1,NSce,Nmodels)
      Real :: DelP(Nmon+1,Nreg,NPeriod-1,NSce,Nmodels)
      Real :: MeanT(Nmon+1,Nreg,NPeriod,NSce,Nmodels)
      Real :: MeanP(Nmon+1,Nreg,NPeriod,NSce,Nmodels)
      !---------------------------------------------------------------------------------
      Character (Len=100) :: Fname, LandMaskFile, RegionFile
      Character (Len=100) :: Data_Dir, Model_nm(Nmodels), Model_Obs, Out_Dir
      Integer :: YSta(NPeriod), YEnd(NPeriod), itime1(NPeriod), itime2(NPeriod)
      !Integer :: YSta(NPeriod), YEnd(NPeriod)
      Character (Len=5) :: Sce_nm(NSce), Period(NPeriod)
      !---------------------------------------------------------------------------------
      !Namelist /Infor/ Data_Dir, Period, YSta, YEnd, Nyear, Model_Obs, Model_nm, Sce_nm, LandMaskFile, Out_Dir
      Namelist /Infor/ Data_Dir, Period, YSta, YEnd, Nyear, Model_Obs, Model_nm, Sce_nm, LandMaskFile, RegionFile, Out_Dir
      Integer :: ilon,ilat,it,imon,iper,imod,isce
      !---------------------------------------------------------------------------------

      Open (1, file="Namelist_R_T_2020_Aug.txt", Status="Old")
      Read (1, Infor)
      Close (1)

      Call Read_Land_Mask
      Ntime = Nyear*Nmon
      itime1(1) = 1; itime2(1) = Ntime
      itime1(2) = (YSta(2) -YEnd(1)-1)*Nmon+1    ! itime1=(2036-2005-1)*12+1=361
      itime2(2) = itime1(2)+Ntime-1
      itime1(3) = (YSta(3) -YEnd(1)-1)*Nmon+1    ! itime2=(2070-2005-1)*12+1=769
      itime2(3) = itime1(3)+Ntime-1

      Allocate (Ob_T2m(Nlon,Nlat,Ntime), Ob_Pr(Nlon,Nlat,Ntime))
      Allocate (BL_T2m(Nlon,Nlat,Ntime,Nmodels), BL_Pr(Nlon,Nlat,Ntime,Nmodels))
      Allocate (MF_T2m(Nlon,Nlat,Ntime,Nmodels,NSce), MF_Pr(Nlon,Nlat,Ntime,Nmodels,NSce))
      Allocate (FF_T2m(Nlon,Nlat,Ntime,Nmodels,NSce), FF_Pr(Nlon,Nlat,Ntime,Nmodels,NSce))

!      Allocate (Ob_T2m(Nlon,Nlat,Ntime), Ob_Pr(Nlon,Nlat,Ntime))
!      Allocate (BL_T2m(Nlon,Nlat,Ntime,Nmodels), BL_Pr(Nlon,Nlat,Ntime,Nmodels))
!      Allocate (MF_T2m(Nlon,Nlat,Ntime,Nmodels,NSce), MF_Pr(Nlon,Nlat,Ntime,Nmodels,NSce))
!      Allocate (FF_T2m(Nlon,Nlat,Ntime,Nmodels,NSce), FF_Pr(Nlon,Nlat,Ntime,Nmodels,NSce))
!      Allocate (T2m(Nlon,Nlat,Ntime)); Allocate (Pr (Nlon,Nlat,Ntime))

      Do iper=1, NPeriod
         if (iper==1) then  ! BL: Obs +6 models + 1 Ens
            Fname=trim(Data_Dir)//trim(Period(iper))//"_"//trim(Model_Obs)//".nc"
            print*,trim(Fname);
            Call Read_T_R (Fname)
            Do ilon=1,Nlon
               Do ilat=1,Nlat
                  Ob_T2m(ilon,ilat,1:Ntime)=T2m(ilon,ilat,1:Ntime)
                  Ob_Pr (ilon,ilat,1:Ntime)=Pr (ilon,ilat,1:Ntime)
               Enddo
            Enddo
            Deallocate (T2m, Pr)
            Do imod=1, Nmodels
               Fname=trim(Data_Dir)//trim(Period(iper))//"_"//trim(Model_nm(imod))//".nc"
               Call Read_T_R (Fname)
               Do ilon=1,Nlon
                  Do ilat=1,Nlat
                     BL_T2m(ilon,ilat,1:Ntime,imod)=T2m(ilon,ilat,1:Ntime)
                     BL_Pr (ilon,ilat,1:Ntime,imod)=Pr (ilon,ilat,1:Ntime)
                  Enddo
               Enddo
               Deallocate (T2m, Pr)
            Enddo
         else
            Do isce=1,NSce
               Do imod=1, Nmodels
                  Fname=trim(Data_Dir)//"SCE_"//trim(Sce_nm(isce))//"_"//trim(Model_nm(imod))//".nc"
                  Call Read_T_R (Fname)
                  it1=itime1(iper); it2=itime1(iper)
                  if (iper==2) then   ! MF
                     Do ilon=1,Nlon
                        Do ilat=1,Nlat
                           MF_T2m(ilon,ilat,1:Ntime,imod,isce)=T2m(ilon,ilat,it1:it2)
                           MF_Pr (ilon,ilat,1:Ntime,imod,isce)=Pr (ilon,ilat,it1:it2)
                        Enddo
                     Enddo
                     !print*,MF_T2m(ilon,ilat,20:30,imod,isce)
                     !MF_T2m(:,:,:,imod,isce)=T2m(:,:,:); MF_Pr(:,:,:,imod,isce)=Pr(:,:,:)
                  else                ! FF
                     Do ilon=1,Nlon
                        Do ilat=1,Nlat
                           FF_T2m(ilon,ilat,1:Ntime,imod,isce)=T2m(ilon,ilat,it1:it2)
                           FF_Pr (ilon,ilat,1:Ntime,imod,isce)=Pr (ilon,ilat,it1:it2)
                        Enddo
                     Enddo
                     !FF_T2m(:,:,:,imod,isce)=T2m(:,:,:); FF_Pr(:,:,:,imod,isce)=Pr(:,:,:)
                  endif
                  Deallocate (T2m, Pr)
               Enddo
            Enddo
         endif
      Enddo

      print*,"Reading data ... Done !"

!      Do ilon=1,Nlon
 !        Do ilat=1, Nlat
 !           if (Mask(ilon,ilat)>0) then
 !              print*,ilon,ilat
 !              print*,Ob_T2m(ilon,ilat,1:Ntime)
 !           endif
 !        Enddo
 !     Enddo
      !pause

      print*,"Calling Chars..."
      Call Chars
      print*,"Calling Chars...Done !"
      print*,"Calling Calc_Regional_Mean..."
      Call Calc_Regional_Mean
      !Call For_Can_Tho

   CONTAINS
   !---------------------------------------------------------------------------------
   Subroutine Calc_Regional_Mean
      Integer :: ilon,ilat,it,imon,imon1,iper,imod,isce, ireg, icount
      Real :: X(Nlon,Nlat)
      Character (Len=100) :: Fname, Sce_nm, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm
      Character (Len=100) :: Fname1, Fname2
      Real, Allocatable :: X2D1(:,:),X2D2(:,:), X3D1(:,:,:), X3D2(:,:,:), Xmod(:), Xmon(:), Xreg(:), XX(:)
      Real, Allocatable :: Y3D1(:,:,:), Y3D2(:,:,:)   ! Dung de sap xep lai cac thang 11,12,1,2,..,10

      Allocate (Ob_T2m_RMon(Nreg,Nmon), Ob_Pr_RMon(Nreg,Nmon))   ! Aug 2020
      Allocate (BL_T2m_RMon(Nreg,Nmon,Nmodels), BL_Pr_RMon(Nreg,Nmon,Nmodels))
      Allocate (MF_T2m_Rmon(Nreg,Nmon,Nmodels,NSce), MF_Pr_RMon(Nreg,Nmon,Nmodels,NSce))
      Allocate (FF_T2m_Rmon(Nreg,Nmon,Nmodels,NSce), FF_Pr_RMon(Nreg,Nmon,Nmodels,NSce))

      Allocate (Ob_T2m_RAnn(Nreg), Ob_Pr_RAnn(Nreg))    ! Aug 2020 
      Allocate (BL_T2m_RAnn(Nreg,Nmodels), BL_Pr_RAnn(Nreg,Nmodels))
      Allocate (MF_T2m_RAnn(Nreg,Nmodels,NSce), MF_Pr_RAnn(Nreg,Nmodels,NSce))
      Allocate (FF_T2m_RAnn(Nreg,Nmodels,NSce), FF_Pr_RAnn(Nreg,Nmodels,NSce))

      Ob_T2m_RMon=Miss; Ob_Pr_RMon=Miss; Ob_T2m_RAnn=Miss; Ob_Pr_RAnn=Miss    ! Aug 2020
      BL_T2m_RMon=Miss; BL_Pr_RMon=Miss; MF_T2m_Rmon=Miss; MF_Pr_RMon=Miss; FF_T2m_Rmon=Miss; FF_Pr_RMon=Miss

      ! Obs Monthly Regional mean:    ! Aug 2020
      Do imon=1, Nmon
         Do ireg=1, Nreg
            X(1:Nlon,1:Nlat)=Ob_T2m_Mean(1:Nlon,1:Nlat,imon); Ob_T2m_RMon(ireg,imon)=Regional_Mean(X,ireg)
            X(1:Nlon,1:Nlat)=Ob_Pr_Mean(1:Nlon,1:Nlat,imon);  Ob_Pr_RMon (ireg,imon)=Regional_Mean(X,ireg)
         Enddo
      Enddo
      ! BL: Monthly Regional mean
      Do imod=1,Nmodels
         Do imon=1, Nmon
            Do ireg=1, Nreg
               X(1:Nlon,1:Nlat)=BL_T2m_Mean(1:Nlon,1:Nlat,imon,imod); BL_T2m_RMon(ireg,imon,imod)=Regional_Mean(X,ireg) !+1)
               X(1:Nlon,1:Nlat)=BL_Pr_Mean (1:Nlon,1:Nlat,imon,imod); BL_Pr_RMon (ireg,imon,imod)=Regional_Mean(X,ireg) !+1) 
            Enddo
         Enddo
      Enddo
      ! SCE: Monthly Regional mean
      Do isce=1, NSce
         Do imod=1,Nmodels
            Do imon=1, Nmon
               Do ireg=1, Nreg
                  X(1:Nlon,1:Nlat)=MF_T2m_Mean(1:Nlon,1:Nlat,imon,imod,isce); MF_T2m_Rmon(ireg,imon,imod,isce)=Regional_Mean(X,ireg) !+1)
                  X(1:Nlon,1:Nlat)=MF_Pr_Mean (1:Nlon,1:Nlat,imon,imod,isce); MF_Pr_RMon (ireg,imon,imod,isce)=Regional_Mean(X,ireg) !+1)
                  X(1:Nlon,1:Nlat)=FF_T2m_Mean(1:Nlon,1:Nlat,imon,imod,isce); FF_T2m_Rmon(ireg,imon,imod,isce)=Regional_Mean(X,ireg) !+1)
                  X(1:Nlon,1:Nlat)=FF_Pr_Mean (1:Nlon,1:Nlat,imon,imod,isce); FF_Pr_RMon (ireg,imon,imod,isce)=Regional_Mean(X,ireg) !+1)
               Enddo
            Enddo
         Enddo
      Enddo
      ! Obs Annual Regional mean    ! Aug 2020
      print*," Obs Annual Regional mean"
      Do ireg=1, Nreg
         X(1:Nlon,1:Nlat)=Ob_T2m_Ann(1:Nlon,1:Nlat); Ob_T2m_RAnn(ireg)=Regional_Mean(X,ireg)
         X(1:Nlon,1:Nlat)=Ob_Pr_Ann(1:Nlon,1:Nlat);  Ob_Pr_RAnn(ireg)=Regional_Mean(X,ireg)
      Enddo
      ! BL: Annual Regional mean
      Do imod=1,Nmodels
         Do ireg=1, Nreg
            X(1:Nlon,1:Nlat)=BL_T2m_Ann(1:Nlon,1:Nlat,imod); BL_T2m_RAnn(ireg,imod)=Regional_Mean(X,ireg) !+1)
            X(1:Nlon,1:Nlat)=BL_Pr_Ann (1:Nlon,1:Nlat,imod); BL_Pr_RAnn (ireg,imod)=Regional_Mean(X,ireg) !+1)
         Enddo
      Enddo
      ! SCE: Annual Regional mean
      Do isce=1, NSce
         Do imod=1,Nmodels
            Do ireg=1, Nreg
               X(1:Nlon,1:Nlat)=MF_T2m_Ann(1:Nlon,1:Nlat,imod,isce); MF_T2m_RAnn(ireg,imod,isce)=Regional_Mean(X,ireg) !+1)
               X(1:Nlon,1:Nlat)=MF_Pr_Ann (1:Nlon,1:Nlat,imod,isce); MF_Pr_RAnn (ireg,imod,isce)=Regional_Mean(X,ireg) !+1)
               X(1:Nlon,1:Nlat)=FF_T2m_Ann(1:Nlon,1:Nlat,imod,isce); FF_T2m_RAnn(ireg,imod,isce)=Regional_Mean(X,ireg) !+1)
               X(1:Nlon,1:Nlat)=FF_Pr_Ann (1:Nlon,1:Nlat,imod,isce); FF_Pr_RAnn (ireg,imod,isce)=Regional_Mean(X,ireg) !+1)
            Enddo
         Enddo
      Enddo
      ! Write Obs regional data to file:
      Allocate (X3D1(Nreg,Nmodels,Nmon),X3D2(Nreg,Nmodels,Nmon), Xmod(Nmodels), Xmon(Nmon), Xreg(Nreg))
      Allocate (Y3D1(Nreg,Nmodels,Nmon),Y3D2(Nreg,Nmodels,Nmon))

      Fname2=trim(Out_Dir)//"Obs_BL_Reg_mean.nc"

      Do ireg=1,Nreg; Xreg(ireg)=ireg; Enddo
      Do imon=1,Nmon; Xmon(imon)=imon; Enddo
      Do imod=1,Nmodels; Xmod(imod)=imod; Enddo

      t_unit="months since 1986-01-01"
      var_ln="Monthly regional mean of Model Temp for period 1986-2005"; lon_nm="reg"; lat_nm="mon"; t_nm="model"
      Call wr_ncdf_3d (Fname2,Xreg,Xmon,Xmod,Nreg,Nmon,Nmodels,BL_T2m_RMon,"t_bl",t_unit, var_ln, "C", lon_nm, lat_nm, t_nm, Miss)
      var_ln="Monthly regional mean of Model Prec for period 1986-2005";
      Call wr_ncdf_3d_apd (Fname2,BL_Pr_RMon,"p_bl",Nreg,Nmon,Nmodels, var_ln, "mm/mon", lon_nm, lat_nm, t_nm, Miss)      

      var_ln="Monthly regional mean of Obs Temp for period 1986-2005"; lon_nm="reg"; lat_nm="mon"
      Call wr_ncdf_2d_apd (Fname2,Ob_T2m_RMon,"tobs",Nreg,Nmon, var_ln,"C", lon_nm, lat_nm, Miss)
      var_ln="Monthly regional mean of Obs Preci for period 1986-2005"
      Call wr_ncdf_2d_apd (Fname2,Ob_Pr_RMon,"pobs",Nreg,Nmon, var_ln, "mm/mon", lon_nm, lat_nm, Miss)
      var_ln="Annual regional mean of Temp for period 1986-2005"
      Call wr_ncdf_1d_apd (Fname2,Ob_T2m_RAnn,"t_ann", Nreg, var_ln, "C", lon_nm, Miss)
      var_ln="Annual regional mean of Preci for period 1986-2005"
      Call wr_ncdf_1d_apd (Fname2,Ob_Pr_RAnn,"p_ann", Nreg, var_ln, "mm/mon", lon_nm, Miss)

      ! Write to file:
      !Allocate (X3D1(Nreg,Nmodels,Nmon),X3D2(Nreg,Nmodels,Nmon), Xmod(Nmodels), Xmon(Nmon), Xreg(Nreg))
      !Allocate (Y3D1(Nreg,Nmodels,Nmon),Y3D2(Nreg,Nmodels,Nmon))

      Do imon=1,Nmon; Xmon(imon)=imon; Enddo
      Do imod=1,Nmodels; Xmod(imod)=imod; Enddo
      Do ireg=1,Nreg; Xreg(ireg)=ireg; Enddo

      t_unit="months since 1986-01-01"; var_unit="C"; lon_nm="reg"; lat_nm="model"; z_nm="scenario"; t_nm="mon"
      var_ln="Monthly regional mean of temperature for the period 1986-2005 projeted by different models"

      Fname=trim(Out_Dir)//"Monthly_Reg_Mean.nc"
      Fname1=trim(Out_Dir)//"Monthly_Reg_Mean_Reshape.nc" ! File co thu tu THANG la 11,12,1,2,...,10

      Do ireg=1,Nreg
         Do imod=1,Nmodels
            X3D1(ireg,imod,1:Nmon)=BL_T2m_RMon(ireg,1:Nmon,imod)
            X3D2(ireg,imod,1:Nmon)=BL_Pr_RMon (ireg,1:Nmon,imod)

            Y3D1(ireg,imod,1:2)   =BL_T2m_RMon(ireg,11:12,imod); Y3D1(ireg,imod,3:12) = BL_T2m_RMon(ireg,1:10,imod)
            Y3D2(ireg,imod,1:2)   =BL_Pr_RMon (ireg,11:12,imod); Y3D2(ireg,imod,3:12) = BL_Pr_RMon (ireg,1:10,imod)
         Enddo
      Enddo
      Call wr_ncdf_3d (Fname,Xreg,Xmod,Xmon,Nreg,Nmodels,Nmon,X3D1,"t_bl",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      Call wr_ncdf_3d (Fname1,Xreg,Xmod,Xmon,Nreg,Nmodels,Nmon,Y3D1,"t_bl",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Fname1

      var_ln="Monthly regional mean of rainfall amount for the period 1986-2005"
      var_unit="mm/mon"
      Call wr_ncdf_3d_apd (Fname,X3D2,"p_bl",Nreg,Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      Call wr_ncdf_3d_apd (Fname1,Y3D2,"p_bl",Nreg,Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)  ! Fname1

      Do isce=1,NSce
         if (isce==1) then; Sce_nm="45"; else; Sce_nm="85"; endif
         Do ireg=1,Nreg
            Do imod=1,Nmodels
               X3D1(ireg,imod,1:Nmon)=MF_T2m_Rmon(ireg,1:Nmon,imod,isce)
               X3D2(ireg,imod,1:Nmon)=MF_Pr_RMon (ireg,1:Nmon,imod,isce)

               Y3D1(ireg,imod,1:2)   =MF_T2m_Rmon(ireg,11:12,imod,isce); Y3D1(ireg,imod,3:12) = MF_T2m_Rmon(ireg,1:10,imod,isce)
               Y3D2(ireg,imod,1:2)   =MF_Pr_RMon (ireg,11:12,imod,isce); Y3D2(ireg,imod,3:12) = MF_Pr_RMon (ireg,1:10,imod,isce)
            Enddo
         Enddo
         var_ln="Monthly regional mean of temperature for the period 2046-2065, RCP"//trim(Sce_nm); var_unit="C"
         Call wr_ncdf_3d_apd (Fname,X3D1,"t_mf"//trim(Sce_nm),Nreg,Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         Call wr_ncdf_3d_apd (Fname1,Y3D1,"t_mf"//trim(Sce_nm),Nreg,Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Fname1

         var_ln="Monthly regional mean of rainfall amount for the period 2046-2065, RCP"//trim(Sce_nm); var_unit="mm/mon"
         Call wr_ncdf_3d_apd (Fname,X3D2,"p_mf"//trim(Sce_nm),Nreg,Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         Call wr_ncdf_3d_apd (Fname1,Y3D2,"p_mf"//trim(Sce_nm),Nreg,Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Fname1

         Do ireg=1,Nreg
            Do imod=1,Nmodels
               X3D1(ireg,imod,1:Nmon)=FF_T2m_Rmon(ireg,1:Nmon,imod,isce)
               X3D2(ireg,imod,1:Nmon)=FF_Pr_RMon (ireg,1:Nmon,imod,isce)

               Y3D1(ireg,imod,1:2)   =FF_T2m_Rmon(ireg,11:12,imod,isce); Y3D1(ireg,imod,3:12) = FF_T2m_Rmon(ireg,1:10,imod,isce)
               Y3D2(ireg,imod,1:2)   =FF_Pr_RMon (ireg,11:12,imod,isce); Y3D2(ireg,imod,3:12) = FF_Pr_RMon (ireg,1:10,imod,isce)
            Enddo
         Enddo
         var_ln="Monthly regional mean of temperature for the period 2080-2099, RCP"//trim(Sce_nm); var_unit="C"
         Call wr_ncdf_3d_apd (Fname,X3D1,"t_ff"//trim(Sce_nm),Nreg,Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         Call wr_ncdf_3d_apd (Fname1,Y3D1,"t_ff"//trim(Sce_nm),Nreg,Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Fname1
         
         var_ln="Monthly regional mean of rainfall amount for the period 2080-2099, RCP"//trim(Sce_nm); var_unit="mm/mon"
         Call wr_ncdf_3d_apd (Fname,X3D2,"p_ff"//trim(Sce_nm),Nreg,Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         Call wr_ncdf_3d_apd (Fname1,Y3D2,"p_ff"//trim(Sce_nm),Nreg,Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) ! Fname1
      Enddo

      Fname=trim(Out_Dir)//"Annual_Reg_Mean.nc"

      ! Baseline:
      var_ln="Annual regional mean of temperature for the period 1986-2005"; var_unit="C"
      Call wr_ncdf_2d (Fname,Xreg,Xmod,Nreg,Nmodels,BL_T2m_RAnn,"t_bl", var_ln, var_unit, lon_nm, lat_nm, Miss)
      var_ln="Annual regional mean of rainfall amount for the period 1986-2005"; var_unit="mm/mon"
      Call wr_ncdf_2d_apd (Fname,BL_Pr_RAnn,"p_bl",Nreg,Nmodels, var_ln, var_unit, lon_nm, lat_nm, Miss)

      ! Obs ann regional mean: Aug 2020
      var_ln="Observed Annual regional mean of temperature for the period 1986-2005"; var_unit="C"
      Call wr_ncdf_1d_apd (Fname,Ob_T2m_RAnn,"tobs", Nreg, var_ln, var_unit, lon_nm, Miss)
      var_ln="Observed Annual regional mean of precipitation for the period 1986-2005"; var_unit="mm/yr"
      Call wr_ncdf_1d_apd (Fname,Ob_Pr_RAnn,"pobs", Nreg, var_ln, var_unit, lon_nm, Miss)

      Allocate (X2D1(Nreg,Nmodels), X2D2(Nreg,Nmodels))
      Do isce=1,NSce
         if (isce==1) then; Sce_nm="45"; else; Sce_nm="85"; endif
         X2D1(:,:)=MF_T2m_RAnn(:,:,isce); X2D2(:,:)=MF_Pr_RAnn(:,:,isce)
         var_ln="Annual regional mean of temperature for the period 2046-2065, RCP"//trim(Sce_nm); var_unit="C"
         Call wr_ncdf_2d_apd (Fname,X2D1,"t_mf"//trim(Sce_nm),Nreg,Nmodels, var_ln, var_unit, lon_nm, lat_nm, Miss)
         var_ln="Monthly regional mean of rainfall amount for the period 2046-2065, RCP"//trim(Sce_nm); var_unit="mm/mon"
         Call wr_ncdf_2d_apd (Fname,X2D2,"p_mf"//trim(Sce_nm),Nreg,Nmodels, var_ln, var_unit, lon_nm, lat_nm, Miss)

         X2D1(:,:)=FF_T2m_RAnn(:,:,isce); X2D2(:,:)=FF_Pr_RAnn(:,:,isce)
         var_ln="Annual regional mean of temperature for the period 2080-2099, RCP"//trim(Sce_nm); var_unit="C"
         Call wr_ncdf_2d_apd (Fname,X2D1,"t_ff"//trim(Sce_nm),Nreg,Nmodels, var_ln, var_unit, lon_nm, lat_nm, Miss)
         var_ln="Monthly regional mean of rainfall amount for the period 2080-2099, RCP"//trim(Sce_nm); var_unit="mm/mon"
         Call wr_ncdf_2d_apd (Fname,X2D2,"p_ff"//trim(Sce_nm),Nreg,Nmodels, var_ln, var_unit, lon_nm, lat_nm, Miss)
      ENddo
   End Subroutine Calc_Regional_Mean
   !---------------------------------------------------------------------------------
   Real Function Regional_Mean (X, ireg)
      Integer :: ilon,ilat,it,imon,imon1,iper,imod,isce, ireg, icount
      Real :: X(Nlon,Nlat), tmp
      Logical :: IsMask

      icount=0; tmp=0
      Do ilon=1, Nlon
         Do ilat=1, Nlat
            !if (ireg>1) then; IsMask=(Mask(ilon,ilat)==ireg); else; IsMask=(Mask(ilon,ilat)>0); endif
            !if (IsMask) then !Mask(ilon,ilat)==ireg) then
            if (Regs(ilon,ilat)==ireg) then
               if (X(ilon,ilat)/=Miss) then
                  icount=icount+1; tmp=tmp+X(ilon,ilat)
               endif
            endif
         Enddo
      Enddo
      if (icount>0) then
         Regional_Mean=tmp/real(icount)
      else
         Regional_Mean=Miss
      endif
   End Function Regional_Mean
   !---------------------------------------------------------------------------------
   Real Function Get_Point_Value(X, KDo, VDo)
      Integer :: ilon, ilat, ilo, ila
      Real :: X(Nlon, Nlat), KDo, VDo
      Logical :: L1, L2

      Do ilon=1, Nlon-1
         Do ilat=1, Nlat-1
            L1=KDo>=xlon(ilon).AND.KDo<=xlon(ilon+1)
            L2=VDo>=xlat(ilat).AND.VDo<=xlat(ilat+1)
            if (L1.AND.L2) then; ilo=ilon; ila=ilat; exit; endif
         Enddo
      Enddo
      Get_Point_Value=X(ilo,ila)
   End Function Get_Point_Value
   !---------------------------------------------------------------------------------
   Real Function P_value(X,Y,N1,N2)   ! Testing whether Xmean and Ymean are different
      ! P_value is P(T>|t|)/2 which is used to compare with alfa=P(T>t_alfa)
      ! P_value < alfa --> |t|>t_alfa --> Ho: Xtb=Ytb is rejected
      Integer :: N1, N2, i, j, n01, n02
      Real :: X(N1), Y(N2)
      Real :: A, t, df, Xtb, Ytb, Sx, Sy, Sx2, Sy2, tmp1, tmp2, tmp3, tmp4, pval

      Xtb=0; n01=0
      Do i=1,N1
         if (X(i) /= Miss) then; Xtb=Xtb+X(i); n01=n01+1; endif
      Enddo
      Ytb=0; n02=0
      Do i=1,N2
         if (Y(i) /= Miss) then; Ytb=Ytb+Y(i); n02=n02+1; endif
      Enddo
      if (n01<5.OR.n02<5) then
         pval=Miss       ! Return with t is missing value
      else
         Sx=Stdev_1(X, N1, Miss); Sy=Stdev_1(Y, N2, Miss)
         if (Sx /= Miss .AND. Sy /= Miss) then
            Sx2=Sx*Sx; Sy2=Sy*Sy     ! S^2
            tmp1=1.0/real(n01)+1.0/real(n02); tmp1=1.0/SQRT(tmp1) ! 1/Sqrt(1/n1 + 1/n2)
            tmp2=SQRT(real(n01+n02-2))
            tmp3=SQRT((n01-1)*Sx2 + (n02-1)*Sy2)
            tmp4=tmp3/tmp2
            A=tmp1/tmp4
            t=ABS(A*(Xtb-Ytb)); df=n01+n02-2
            pval=T_Dist (t, df)
         else
            pval=Miss
         endif
      endif
      P_value=pval
      !print"(A,F8.3,F5.0,F10.5)","t, df, p=",t, df,pval
   End Function P_value
   !---------------------------------------------------------------------------------
   Subroutine Chars
      Integer :: ilon,ilat,it,imon,imon1,iper,imod,isce
      Integer :: iyr
      Real, Allocatable :: X(:),Y(:),X1(:),Y1(:), X3D(:,:,:), X4D(:,:,:,:), Xmod(:), Xmon(:), XX(:), YY(:)
      Real, Allocatable :: XX1(:), YY1(:), X0(:), Y0(:), XX0(:), YY0(:)
      Real :: Tmp
      Character (Len=100) :: Fname, Sce_nm, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm

      !print*,Nlon,Nlat,Nmodels,Nmon,Ntime,NSce, Nyear; pause
      Allocate (Ob_T2m_Mean(Nlon,Nlat,Nmon), Ob_Pr_Mean(Nlon,Nlat,Nmon))   ! Aug 2020
      Allocate (BL_T2m_Mean(Nlon,Nlat,Nmon,Nmodels),BL_Pr_Mean(Nlon,Nlat,Nmon,Nmodels))
      Allocate (MF_T2m_Mean(Nlon,Nlat,Nmon,Nmodels,NSce),MF_Pr_Mean(Nlon,Nlat,Nmon,Nmodels,NSce))
      Allocate (FF_T2m_Mean(Nlon,Nlat,Nmon,Nmodels,NSce),FF_Pr_Mean(Nlon,Nlat,Nmon,Nmodels,NSce))

      Allocate (Ob_T2m_Ann(Nlon,Nlat), Ob_Pr_Ann(Nlon,Nlat))    ! Aug 2020
      Allocate (BL_T2m_Ann(Nlon,Nlat,Nmodels),BL_Pr_Ann(Nlon,Nlat,Nmodels))
      Allocate (MF_T2m_Ann(Nlon,Nlat,Nmodels,NSce),MF_Pr_Ann(Nlon,Nlat,Nmodels,NSce))
      Allocate (FF_T2m_Ann(Nlon,Nlat,Nmodels,NSce),FF_Pr_Ann(Nlon,Nlat,Nmodels,NSce))
      
      ! Changes of Annual mean:
      Allocate (D_BL_T2m_Ann(Nlon,Nlat,Nmodels), D_BL_Pr_Ann(Nlon,Nlat,Nmodels))   ! Aug 2020
      Allocate (D_MF_T2m_Ann(Nlon,Nlat,Nmodels,NSce),D_MF_Pr_Ann(Nlon,Nlat,Nmodels,NSce))
      Allocate (D_FF_T2m_Ann(Nlon,Nlat,Nmodels,NSce),D_FF_Pr_Ann(Nlon,Nlat,Nmodels,NSce))
      ! P-value for Changes of Annual mean:
      Allocate (P_BL_T2m_Ann(Nlon,Nlat,Nmodels), P_BL_Pr_Ann(Nlon,Nlat,Nmodels))   ! Aug 2020
      Allocate (P_MF_T2m_Ann(Nlon,Nlat,Nmodels,NSce),P_MF_Pr_Ann(Nlon,Nlat,Nmodels,NSce))
      Allocate (P_FF_T2m_Ann(Nlon,Nlat,Nmodels,NSce),P_FF_Pr_Ann(Nlon,Nlat,Nmodels,NSce))

      ! Changes of Monthly mean:
      Allocate (D_BL_T2m_Mon(Nlon,Nlat,Nmon,Nmodels), D_BL_Pr_Mon(Nlon,Nlat,Nmon,Nmodels))   ! Aug 2020
      Allocate (D_MF_T2m_Mon(Nlon,Nlat,Nmon,Nmodels,NSce),D_MF_Pr_Mon(Nlon,Nlat,Nmon,Nmodels,NSce))
      Allocate (D_FF_T2m_Mon(Nlon,Nlat,Nmon,Nmodels,NSce),D_FF_Pr_Mon(Nlon,Nlat,Nmon,Nmodels,NSce))
      ! P-value for Changes of Monthly mean:
      Allocate (P_BL_T2m_Mon(Nlon,Nlat,Nmon,Nmodels), P_BL_Pr_Mon(Nlon,Nlat,Nmon,Nmodels))   ! Aug 2020
      Allocate (P_MF_T2m_Mon(Nlon,Nlat,Nmon,Nmodels,NSce),P_MF_Pr_Mon(Nlon,Nlat,Nmon,Nmodels,NSce))
      Allocate (P_FF_T2m_Mon(Nlon,Nlat,Nmon,Nmodels,NSce),P_FF_Pr_Mon(Nlon,Nlat,Nmon,Nmodels,NSce))

      Allocate (X(Nyear),Y(Nyear),X1(Nyear),Y1(Nyear), XX(Ntime), YY(Ntime)) ! Xmon(Nmon), XX(Ntime))
      Allocate (XX1(Ntime), YY1(Ntime), X0(Nyear),Y0(Nyear),XX0(Ntime), YY0(Ntime))

      Ob_T2m_Mean=Miss; Ob_Pr_Mean=Miss   ! Aug 2020
      BL_T2m_Mean=Miss; BL_Pr_Mean=Miss; MF_T2m_Mean=Miss; MF_Pr_Mean=Miss; FF_T2m_Mean=Miss; FF_Pr_Mean=Miss

      Ob_T2m_Ann=Miss; Ob_Pr_Ann=Miss   ! Aug 2020
      BL_T2m_Ann=Miss; BL_Pr_Ann=Miss; MF_T2m_Ann=Miss; MF_Pr_Ann=Miss; FF_T2m_Ann=Miss; FF_Pr_Ann=Miss

      D_BL_T2m_Mon=Miss; D_BL_Pr_Mon=Miss   ! Aug 2020
      D_MF_T2m_Mon=Miss; D_MF_Pr_Mon=Miss; D_FF_T2m_Mon=Miss; D_FF_Pr_Mon=Miss
      P_BL_T2m_Mon=Miss; P_BL_Pr_Mon=Miss   ! Jan 2021: P-value
      P_MF_T2m_Mon=Miss; P_MF_Pr_Mon=Miss; P_FF_T2m_Mon=Miss; P_FF_Pr_Mon=Miss

      D_BL_T2m_Ann=Miss; D_BL_Pr_Ann=Miss   ! Aug 2020
      D_MF_T2m_Ann=Miss; D_MF_Pr_Ann=Miss; D_FF_T2m_Ann=Miss; D_FF_Pr_Ann=Miss
      P_BL_T2m_Ann=Miss; P_BL_Pr_Ann=Miss   ! Jan 2021: P-value
      P_MF_T2m_Ann=Miss; P_MF_Pr_Ann=Miss; P_FF_T2m_Ann=Miss; P_FF_Pr_Ann=Miss

      Do ilon=1, Nlon
         Do ilat=1, Nlat
            if (Mask(ilon,ilat)>0) then
               print*,ilon,ilat,Nmodels,Nmon,Ntime
               ! Obs:   ! Aug 2020
               Do imon=1, Nmon
                  iyr=0
                  Do it=1, Ntime
                     imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon
                     if (imon1==imon) then
                        iyr=iyr+1; X(iyr)=Ob_T2m(ilon,ilat,it); Y(iyr)=Ob_Pr(ilon,ilat,it)
                     endif
                  Enddo
                  Ob_T2m_Mean(ilon,ilat,imon)=Ave(X,iyr,Miss)
                  Ob_Pr_Mean(ilon,ilat,imon) =Ave(Y,iyr,Miss)
               Enddo
               !print "(12F7.1)",Ob_T2m_Mean(ilon,ilat,:); pause
               !Obs Annual mean:
               XX(1:Ntime)=Ob_T2m(ilon,ilat,1:Ntime); Ob_T2m_Ann(ilon,ilat)=Ave(XX,Ntime,Miss)
               XX(1:Ntime)=Ob_Pr (ilon,ilat,1:Ntime); Ob_Pr_Ann (ilon,ilat)=Ave(XX,Ntime,Miss)
               XX1(1:Ntime)=Ob_T2m(ilon,ilat,1:Ntime); YY1(1:Ntime)=Ob_Pr (ilon,ilat,1:Ntime)    ! For annualy P-value
               if (Ob_Pr_Ann (ilon,ilat)/=Miss) Ob_Pr_Ann (ilon,ilat)=Ob_Pr_Ann (ilon,ilat)*Nmon ! Ann total rainfall
               ! Baseline:
               print*,"Baseline"
               Do imod=1,Nmodels
                  Do imon=1, Nmon
                     iyr=0
                     Do it=1, Ntime
                        imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon
                        if (imon1==imon) then
                           iyr=iyr+1; X(iyr)=BL_T2m(ilon,ilat,it,imod); Y(iyr)=BL_Pr(ilon,ilat,it,imod)
                                     X1(iyr)=Ob_T2m(ilon,ilat,it);     Y1(iyr)=Ob_Pr(ilon,ilat,it)      ! For monthly P-value
                           !print*,"iyr=",iyr
                        endif
                     Enddo
                     ! BL Monthly mean:
                     BL_T2m_Mean(ilon,ilat,imon,imod)=Ave(X,iyr,Miss)
                     BL_Pr_Mean(ilon,ilat,imon,imod)=Ave(Y,iyr,Miss)

                     ! P-value for monthly:
                     P_BL_T2m_Mon(ilon,ilat,imon,imod)=P_value(X,X1,iyr,iyr) ! So sanh Trung binh 2 chuoi
                     P_BL_Pr_Mon (ilon,ilat,imon,imod)=P_value(Y,Y1,iyr,iyr)

                     ! BL: Differences BL & Obs monthly:   ! Aug 2020
                     D_BL_T2m_Mon(ilon,ilat,imon,imod)=BL_T2m_Mean(ilon,ilat,imon,imod)-Ob_T2m_Mean(ilon,ilat,imon)
                     D_BL_Pr_Mon (ilon,ilat,imon,imod)=(BL_Pr_Mean(ilon,ilat,imon,imod)-Ob_Pr_Mean(ilon,ilat,imon))/ &
                                                        Ob_Pr_Mean(ilon,ilat,imon)*100.0
                  Enddo
                  ! BL Annual mean:
                  XX(1:Ntime)=BL_T2m(ilon,ilat,1:Ntime,imod); BL_T2m_Ann(ilon,ilat,imod)=Ave(XX,Ntime,Miss)
                  XX(1:Ntime)=BL_Pr (ilon,ilat,1:Ntime,imod); BL_Pr_Ann (ilon,ilat,imod)=Ave(XX,Ntime,Miss)
                  if (BL_Pr_Ann (ilon,ilat,imod)/=Miss) BL_Pr_Ann (ilon,ilat,imod)=BL_Pr_Ann (ilon,ilat,imod)*Nmon ! Ann total rainfall
                  ! BL: Differences BL & Obs Ann:
                  D_BL_T2m_Ann(ilon,ilat,imod)=BL_T2m_Ann(ilon,ilat,imod)-Ob_T2m_Ann(ilon,ilat)
                  D_BL_pr_Ann (ilon,ilat,imod)=(BL_Pr_Ann(ilon,ilat,imod)-Ob_Pr_Ann (ilon,ilat))/ &
                                                Ob_Pr_Ann (ilon,ilat)*100.0
                  ! P-value for Annualy:
                  XX(1:Ntime)=BL_T2m(ilon,ilat,1:Ntime,imod); YY(1:Ntime)=BL_Pr (ilon,ilat,1:Ntime,imod)
                  P_BL_T2m_Ann(ilon,ilat,imod)=P_value(XX,XX1,Ntime,Ntime)
                  P_BL_pr_Ann (ilon,ilat,imod)=P_value(YY,YY1,Ntime,Ntime)
               Enddo
               ! Scenarios:
               print*,"Scenarios"
               Do isce=1,NSce
                  Do imod=1,Nmodels
                     Do imon=1, Nmon
                        iyr=0
                        Do it=1, Ntime
                           imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon
                           if (imon1==imon) then
                              iyr=iyr+1; 
                              X0(iyr)=BL_T2m(ilon,ilat,it,imod);      Y0(iyr)=BL_Pr(ilon,ilat,it,imod)
                              X(iyr) =MF_T2m(ilon,ilat,it,imod,isce); Y(iyr) =MF_Pr(ilon,ilat,it,imod,isce)
                              X1(iyr)=FF_T2m(ilon,ilat,it,imod,isce); Y1(iyr)=FF_Pr(ilon,ilat,it,imod,isce)
                           endif
                        Enddo
                        MF_T2m_Mean(ilon,ilat,imon,imod,isce) = Ave(X, iyr,Miss)
                        MF_Pr_Mean (ilon,ilat,imon,imod,isce) = Ave(Y, iyr,Miss)
                        FF_T2m_Mean(ilon,ilat,imon,imod,isce) = Ave(X1,iyr,Miss)
                        FF_Pr_Mean (ilon,ilat,imon,imod,isce) = Ave(Y1,iyr,Miss)

                        ! P-value for monthly Sce:
                        P_MF_T2m_Mon(ilon,ilat,imon,imod,isce)= P_value(X0,X,iyr,iyr)
                        P_MF_Pr_Mon (ilon,ilat,imon,imod,isce)= P_value(Y0,Y,iyr,iyr)
                        P_FF_T2m_Mon(ilon,ilat,imon,imod,isce)= P_value(X0,X1,iyr,iyr)
                        P_FF_Pr_Mon (ilon,ilat,imon,imod,isce)= P_value(Y0,Y1,iyr,iyr)

                        ! Changes of Monthly mean:
                        D_MF_T2m_Mon(ilon,ilat,imon,imod,isce)= MF_T2m_Mean(ilon,ilat,imon,imod,isce)-BL_T2m_Mean(ilon,ilat,imon,imod)
                        D_MF_Pr_Mon (ilon,ilat,imon,imod,isce)=(MF_Pr_Mean (ilon,ilat,imon,imod,isce)-BL_Pr_Mean (ilon,ilat,imon,imod))/ &
                                                                BL_Pr_Mean (ilon,ilat,imon,imod)*100.0
                        D_FF_T2m_Mon(ilon,ilat,imon,imod,isce)= FF_T2m_Mean(ilon,ilat,imon,imod,isce)-BL_T2m_Mean(ilon,ilat,imon,imod)
                        D_FF_Pr_Mon (ilon,ilat,imon,imod,isce)=(FF_Pr_Mean (ilon,ilat,imon,imod,isce)-BL_Pr_Mean (ilon,ilat,imon,imod))/ &
                                                                BL_Pr_Mean (ilon,ilat,imon,imod)*100.0
                     Enddo
                     XX(1:Ntime)=MF_T2m(ilon,ilat,1:Ntime,imod,isce); MF_T2m_Ann(ilon,ilat,imod,isce)=Ave(XX,Ntime,Miss)
                     XX(1:Ntime)=MF_Pr (ilon,ilat,1:Ntime,imod,isce); MF_Pr_Ann (ilon,ilat,imod,isce)=Ave(XX,Ntime,Miss)
                     if (MF_Pr_Ann (ilon,ilat,imod,isce)/=Miss) MF_Pr_Ann (ilon,ilat,imod,isce)=MF_Pr_Ann (ilon,ilat,imod,isce)*Nmon ! Ann total rainfall
                     XX(1:Ntime)=FF_T2m(ilon,ilat,1:Ntime,imod,isce); FF_T2m_Ann(ilon,ilat,imod,isce)=Ave(XX,Ntime,Miss)
                     XX(1:Ntime)=FF_Pr (ilon,ilat,1:Ntime,imod,isce); FF_Pr_Ann (ilon,ilat,imod,isce)=Ave(XX,Ntime,Miss)
                     if (FF_Pr_Ann (ilon,ilat,imod,isce)/=Miss) FF_Pr_Ann (ilon,ilat,imod,isce)=FF_Pr_Ann (ilon,ilat,imod,isce)*Nmon ! Ann total rainfall 
                     ! Changes of Annual mean:
                     D_MF_T2m_Ann(ilon,ilat,imod,isce)= MF_T2m_Ann(ilon,ilat,imod,isce)-BL_T2m_Ann(ilon,ilat,imod)
                     D_MF_Pr_Ann (ilon,ilat,imod,isce)=(MF_Pr_Ann (ilon,ilat,imod,isce)-BL_Pr_Ann (ilon,ilat,imod))/ &
                                                        BL_Pr_Ann (ilon,ilat,imod)*100.0
                     D_FF_T2m_Ann(ilon,ilat,imod,isce)= FF_T2m_Ann(ilon,ilat,imod,isce)-BL_T2m_Ann(ilon,ilat,imod)
                     D_FF_Pr_Ann (ilon,ilat,imod,isce)=(FF_Pr_Ann (ilon,ilat,imod,isce)-BL_Pr_Ann (ilon,ilat,imod))/ &
                                                        BL_Pr_Ann (ilon,ilat,imod)*100.0

                     ! P-value for Annualy & Sce:
                     XX0(1:Ntime)=BL_T2m(ilon,ilat,1:Ntime,imod);      YY0(1:Ntime)=BL_Pr (ilon,ilat,1:Ntime,imod)
                     XX(1:Ntime) =MF_T2m(ilon,ilat,1:Ntime,imod,isce); YY(1:Ntime) =MF_Pr (ilon,ilat,1:Ntime,imod,isce)
                     XX1(1:Ntime)=FF_T2m(ilon,ilat,1:Ntime,imod,isce); YY1(1:Ntime)=FF_Pr (ilon,ilat,1:Ntime,imod,isce)
                     P_MF_T2m_Ann(ilon,ilat,imod,isce)=P_value(XX0,XX,Ntime,Ntime)
                     P_MF_Pr_Ann (ilon,ilat,imod,isce)=P_value(YY0,YY,Ntime,Ntime)
                     P_FF_T2m_Ann(ilon,ilat,imod,isce)=P_value(XX0,XX1,Ntime,Ntime)
                     P_FF_Pr_Ann (ilon,ilat,imod,isce)=P_value(YY0,YY1,Ntime,Ntime)
                  Enddo
               Enddo
               print*,ilon,ilat
            endif
         Enddo
      Enddo
      print*,"Deallocate..."; 
      Deallocate (X,Y,X1,Y1,XX) ! ,Xmon)

      print*,"Writing in Chars ..."; pause

      ! Write to files: Monthly Mean
      Allocate (X4D(Nlon,Nlat,Nmodels,Nmon), Xmod(Nmodels), Xmon(Nmon))

      Do imon=1,Nmon; Xmon(imon)=imon; Enddo
      Do imod=1,Nmodels; Xmod(imod)=imod; Enddo
      ! BL Mean Temp
      Do ilon=1, Nlon
         Do ilat=1, Nlat
            Do imod=1,Nmodels 
               X4D(ilon,ilat,imod,1:Nmon)=BL_T2m_Mean(ilon,ilat,1:Nmon,imod)
            Enddo
         Enddo
      Enddo
      t_unit="months since 1986-01-01"; var_unit="C"; lon_nm="lon"; lat_nm="lat"; z_nm="model"; t_nm="mon"
      var_ln="Monthly mean temperature for the period 1986-2005 projeted by different models"

      !-------------------------------------------------------------------------
      Fname=trim(Out_Dir)//"Monthly_Mean.nc"
      !-------------------------------------------------------------------------

      Call wr_ncdf_4d (Fname,xlon,xlat,Xmod,Xmon,nlon,nlat,Nmodels,Nmon, X4D,"t_bl", &
                       t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)

      ! BL Mean Prec
      Do ilon=1, Nlon
         Do ilat=1, Nlat
            Do imod=1,Nmodels
               X4D(ilon,ilat,imod,1:Nmon)=BL_Pr_Mean(ilon,ilat,1:Nmon,imod)
            Enddo
         Enddo
      Enddo
      var_ln="Monthly mean rainfall amount for the period 1986-2005 projeted by different models"; var_unit="mm/mon"
      Call wr_ncdf_4d_apd (Fname, X4D,"p_bl", nlon,nlat,Nmodels,Nmon, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)

      ! P-value for differences between BL & Obs of Monthly mean:
      var_ln="P-values for differences between BL & Obs of Monthly Temperature"; var_unit="-"
      Call wr_ncdf_4d_apd (Fname, P_BL_T2m_Mon,"pval_t_bl", nlon,nlat,Nmodels,Nmon, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
      var_ln="P-values for differences between BL & Obs of Monthly Precipitation"; var_unit="-"
      Call wr_ncdf_4d_apd (Fname, P_BL_Pr_Mon,"pval_p_bl", nlon,nlat,Nmodels,Nmon, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
      !---------------------------------------------
      ! Obs Mean Temp - Aug 2020
      var_ln="Observed monthly mean temperature for the period 1986-2005"; var_unit="C"
      Call wr_ncdf_3d_apd (Fname,Ob_T2m_Mean,"tobs",nlon,nlat,Nmon, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      ! Obs monthly Prec - Aug 2020
      var_ln="Observed monthly precipitation for the period 1986-2005"; var_unit="mm/mon"
      Call wr_ncdf_3d_apd (Fname,Ob_Pr_Mean,"pobs",nlon,nlat,Nmon, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      !---------------------------------------------
      print*,"Here ==>"

      Do isce=1,NSce
         if (isce==1) then; Sce_nm="45"; else; Sce_nm="85"; endif
         ! MF mean Temp
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               Do imod=1,Nmodels
                  X4D(ilon,ilat,imod,1:Nmon)=MF_T2m_Mean(ilon,ilat,1:Nmon,imod,isce)
               Enddo
            Enddo
         Enddo
         var_ln="Monthly mean temperature for the period 2046-2065 projeted by different models and scenario of RCP"//trim(Sce_nm); var_unit="C"
         Call wr_ncdf_4d_apd (Fname, X4D, "t_mf"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         ! MF mean Prec
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               Do imod=1,Nmodels
                  X4D(ilon,ilat,imod,1:Nmon)=MF_Pr_Mean(ilon,ilat,1:Nmon,imod,isce)
               Enddo          
            Enddo
         Enddo
         var_ln="Monthly mean rainfall amount for the period 2046-2065 projeted by different models and scenario of RCP"//trim(Sce_nm); var_unit="mm/mon"
         Call wr_ncdf_4d_apd (Fname, X4D, "p_mf"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         ! FF mean Temp
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               Do imod=1,Nmodels
                  X4D(ilon,ilat,imod,1:Nmon)=FF_T2m_Mean(ilon,ilat,1:Nmon,imod,isce)
               Enddo
            Enddo
         Enddo
         var_ln="Monthly mean temperature for the period 2080-2099 projeted by different models and scenario of RCP"//trim(Sce_nm); var_unit="C"
         Call wr_ncdf_4d_apd (Fname, X4D, "t_ff"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         ! FF mean Prec
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               Do imod=1,Nmodels
                  X4D(ilon,ilat,imod,1:Nmon)=FF_Pr_Mean(ilon,ilat,1:Nmon,imod,isce)
               Enddo
            Enddo
         Enddo
         var_ln="Monthly mean rainfall amount for the period 2080-2099 projeted by different models and scenario of RCP"//trim(Sce_nm); var_unit="mm/mon"
         Call wr_ncdf_4d_apd (Fname, X4D, "p_ff"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss) 

         ! P-value for differences between Sce & BL of Monthly mean:
         var_ln="P-values for monthly temp. changes for period 2046-2065 and scenario of RCP"//trim(Sce_nm); var_unit="-"
         Call wr_ncdf_4d_apd (Fname,P_MF_T2m_Mon,"pval_t_mf"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         var_ln="P-values for monthly rainfall changes for period 2046-2065 and scenario of RCP"//trim(Sce_nm); var_unit="-"
         Call wr_ncdf_4d_apd (Fname,P_MF_Pr_Mon,"pval_p_mf"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         var_ln="P-values for monthly temp. changes for period 2080-2099 and scenario of RCP"//trim(Sce_nm); var_unit="-"
         Call wr_ncdf_4d_apd (Fname,P_FF_T2m_Mon,"pval_t_ff"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         var_ln="P-values for monthly rainfall changes for period 2080-2099 and scenario of RCP"//trim(Sce_nm); var_unit="-"
         Call wr_ncdf_4d_apd (Fname,P_FF_Pr_Mon,"pval_p_ff"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
      Enddo 
      ! Write to files: Annual mean
      !-------------------------------------------------------------------------
      Fname=trim(Out_Dir)//"Annual_Mean.nc" 
      !-------------------------------------------------------------------------
      Allocate (X3D(Nlon,Nlat,Nmodels))

      ! Testing only:
      !========================================================
      !X3D(:,:,:)=BL_T2m_Ann(:,:,:)
      !BL_T2m_Ann(:,:,1:3)=X3D(:,:,5:7); BL_T2m_Ann(:,:,4:7)=X3D(:,:,1:4)
      !X3D(:,:,:)=BL_Pr_Ann(:,:,:)
      !BL_Pr_Ann(:,:,1:3)=X3D(:,:,5:7); BL_Pr_Ann(:,:,4:7)=X3D(:,:,1:4)

      !X4D(:,:,:,:)=MF_T2m_Ann(:,:,:,:)
      !MF_T2m_Ann(:,:,1:3,:)=X4D(:,:,5:7,:); MF_T2m_Ann(:,:,4:7,:)=X4D(:,:,1:4,:)
      !X4D(:,:,:,:)=MF_Pr_Ann(:,:,:,:)
      !MF_Pr_Ann(:,:,1:3,:)=X4D(:,:,5:7,:); MF_Pr_Ann(:,:,4:7,:)=X4D(:,:,1:4,:)

      !X4D(:,:,:,:)=FF_T2m_Ann(:,:,:,:)
      !FF_T2m_Ann(:,:,1:3,:)=X4D(:,:,5:7,:); FF_T2m_Ann(:,:,4:7,:)=X4D(:,:,1:4,:)
      !X4D(:,:,:,:)=FF_Pr_Ann(:,:,:,:)
      !FF_Pr_Ann(:,:,1:3,:)=X4D(:,:,5:7,:); FF_Pr_Ann(:,:,4:7,:)=X4D(:,:,1:4,:)
      ! End of testing

      !========================================================
      ! BL Ann:
      var_ln="Annual mean temperature for the period 1986-2005"; t_nm="model"
      X3D(:,:,:)=BL_T2m_Ann(:,:,:)
      Call wr_ncdf_3d (Fname,xlon,xlat,Xmod,nlon,nlat,Nmodels,X3D,"t_bl",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss) 
      !Call wr_ncdf_3d (Fname,xlon,xlat,Xmod,nlon,nlat,Nmodels,BL_T2m_Ann,"t_bl",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Annual total rainfall amount for the period 1986-2005"
      X3D(:,:,:)=BL_Pr_Ann(:,:,:)
      Call wr_ncdf_3d_apd (Fname,X3D,"p_bl",nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

      ! P-value for differences between BL & Obs of Annual mean:
      X3D(:,:,:)=P_BL_T2m_Ann(:,:,:)
      var_ln="P-values for differences between T2m & Obs for baseline period"
      Call wr_ncdf_3d_apd (Fname,X3D,"pval_t_bl",nlon,nlat,Nmodels, var_ln,"-", lon_nm, lat_nm, t_nm, Miss)
      X3D(:,:,:)=P_BL_Pr_Ann(:,:,:)
      var_ln="P-values for differences between Pre & Obs for baseline period"
      Call wr_ncdf_3d_apd (Fname,X3D,"pval_p_bl",nlon,nlat,Nmodels, var_ln,"-", lon_nm, lat_nm, t_nm, Miss)

      ! Obs Ann: Aug 2020
      var_ln="Observed annual mean temperature for the period 1986-2005"
      Call wr_ncdf_2d_apd (Fname, Ob_T2m_Ann,"tannob",nlon,nlat, var_ln, "C", lon_nm, lat_nm, Miss)
      var_ln="Observed annual prec for the period 1986-2005"
      Call wr_ncdf_2d_apd (Fname, Ob_Pr_Ann,"pannob",nlon,nlat, var_ln, "C", lon_nm, lat_nm, Miss)

      !Call wr_ncdf_3d_apd (Fname,BL_Pr_Ann,"p_bl",nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)    

      !Allocate (X3D(Nlon,Nlat,Nmodels))
      Do isce=1,NSce
         if (isce==1) then; Sce_nm="45"; else; Sce_nm="85"; endif
         X3D(:,:,:)=MF_T2m_Ann(:,:,:,isce)
         var_ln="Annual mean temperature for the period 2046-2065, RCP"//trim(Sce_nm); var_unit="C"
         Call wr_ncdf_3d_apd (Fname,X3D,"t_mf"//trim(Sce_nm),nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         X3D(:,:,:)=MF_Pr_Ann(:,:,:,isce)
         var_ln="Annual total rainfall amount for the period 2046-2065, RCP"//trim(Sce_nm); var_unit="mm/mon"
         Call wr_ncdf_3d_apd (Fname,X3D,"p_mf"//trim(Sce_nm),nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         X3D(:,:,:)=FF_T2m_Ann(:,:,:,isce)
         var_ln="Annual mean temperature for the period 2080-2099, RCP"//trim(Sce_nm); var_unit="C"
         Call wr_ncdf_3d_apd (Fname,X3D,"t_ff"//trim(Sce_nm),nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         X3D(:,:,:)=FF_Pr_Ann(:,:,:,isce)
         var_ln="Annual total rainfall amount for the period 2080-2099, RCP"//trim(Sce_nm); var_unit="mm/mon"
         Call wr_ncdf_3d_apd (Fname,X3D,"p_ff"//trim(Sce_nm),nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

         ! P-value for Changes of Annual mean:
         X3D(:,:,:)=P_MF_T2m_Ann(:,:,:,isce)
         var_ln="P-values for changes of Ann. temp. for the period 2046-2065, RCP"//trim(Sce_nm); var_unit="-"
         Call wr_ncdf_3d_apd (Fname,X3D,"pval_t_mf"//trim(Sce_nm),nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         X3D(:,:,:)=P_MF_Pr_Ann(:,:,:,isce)
         var_ln="P-values for changes of Ann. Pre for the period 2046-2065, RCP"//trim(Sce_nm); var_unit="-"
         Call wr_ncdf_3d_apd (Fname,X3D,"pval_p_mf"//trim(Sce_nm),nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         X3D(:,:,:)=P_FF_T2m_Ann(:,:,:,isce)
         var_ln="P-values for changes of Ann. temp. for the period 2080-2099, RCP"//trim(Sce_nm); var_unit="-"
         Call wr_ncdf_3d_apd (Fname,X3D,"pval_t_ff"//trim(Sce_nm),nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         X3D(:,:,:)=P_FF_Pr_Ann(:,:,:,isce)
         var_ln="P-values for changes of Ann. Pre for the period 2080-2099, RCP"//trim(Sce_nm); var_unit="-"
         Call wr_ncdf_3d_apd (Fname,X3D,"pval_p_ff"//trim(Sce_nm),nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      Enddo

      ! Changes in Temp and Prec:
      ! Changes of Annual mean:
      !-------------------------------------------------------------------------
      Fname=trim(Out_Dir)//"Change_Annual_Mean.nc"
      !-------------------------------------------------------------------------
      Do isce=1,NSce
         if (isce==1) then; Sce_nm="45"; else; Sce_nm="85"; endif
         var_ln="Changes in Annual mean temperature for the period 2046-2065, RCP"//trim(Sce_nm); t_nm="model"; var_unit="C"
         X3D(:,:,:)=D_MF_T2m_Ann(:,:,:,isce)
         if (isce==1) then
            Call wr_ncdf_3d (Fname,xlon,xlat,Xmod,nlon,nlat,Nmodels,X3D,"t_mf"//trim(Sce_nm),t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         else
            Call wr_ncdf_3d_apd (Fname,X3D,"t_mf"//trim(Sce_nm),nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         endif
         var_ln="Changes in Annual mean rainfall for the period 2046-2065, RCP"//trim(Sce_nm); t_nm="model"; var_unit="mm/mon"
         X3D(:,:,:)=D_MF_Pr_Ann(:,:,:,isce)
         Call wr_ncdf_3d_apd (Fname,X3D,"p_mf"//trim(Sce_nm),nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

         var_ln="Changes in Annual mean temperature for the period 2080-2099, RCP"//trim(Sce_nm); t_nm="model"; var_unit="C"
         X3D(:,:,:)=D_FF_T2m_Ann(:,:,:,isce)
         Call wr_ncdf_3d_apd (Fname,X3D,"t_ff"//trim(Sce_nm),nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
         var_ln="Changes in Annual mean rainfall for the period 2080-2099, RCP"//trim(Sce_nm); t_nm="model"; var_unit="mm/mon"
         X3D(:,:,:)=D_FF_Pr_Ann(:,:,:,isce)
         Call wr_ncdf_3d_apd (Fname,X3D,"p_ff"//trim(Sce_nm),nlon,nlat,Nmodels, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      Enddo
      ! BL: Obs-BL - Aug 2020: D_BL_T2m_Ann
      var_ln="Differences in Annual mean temperature between Obs and BL for the period 1986-2005"
      Call wr_ncdf_3d_apd (Fname,D_BL_T2m_Ann,"delt",nlon,nlat,Nmodels, var_ln, "C", lon_nm, lat_nm, t_nm, Miss)
      var_ln="Differences in Annual preci between Obs and BL for the period 1986-2005"
      Call wr_ncdf_3d_apd (Fname,D_BL_PR_Ann,"delp",nlon,nlat,Nmodels, var_ln, "mm/yr",lon_nm, lat_nm, t_nm, Miss)

      ! Changes of Monthly mean:
      !-------------------------------------------------------------------------
      Fname=trim(Out_Dir)//"Change_Monthly_Mean.nc"
      !-------------------------------------------------------------------------
      z_nm="model"; t_nm="mon"
      Do isce=1,NSce
         if (isce==1) then; Sce_nm="45"; else; Sce_nm="85"; endif
         ! MF mean Temp
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               Do imod=1,Nmodels
                  X4D(ilon,ilat,imod,1:Nmon)=D_MF_T2m_Mon(ilon,ilat,1:Nmon,imod,isce)
               Enddo
            Enddo
         Enddo
         var_ln="Change in Monthly mean temperature for the period 2046-2065, RCP"//trim(Sce_nm); var_unit="C"
         if (isce==1) then
            Call wr_ncdf_4d (Fname,xlon,xlat,Xmod,Xmon,nlon,nlat,Nmodels,Nmon, X4D,"t_mf"//trim(Sce_nm), &
                       t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
         else
            Call wr_ncdf_4d_apd (Fname, X4D, "t_mf"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         endif
         ! MF mean Prec
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               Do imod=1,Nmodels
                  X4D(ilon,ilat,imod,1:Nmon)=D_MF_Pr_Mon(ilon,ilat,1:Nmon,imod,isce)
               Enddo
            Enddo
         Enddo
         var_ln="Change in Monthly mean rainfall amount for the period 2046-2065,RCP"//trim(Sce_nm); var_unit="mm/mon"
         Call wr_ncdf_4d_apd (Fname, X4D, "p_mf"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         ! FF mean Temp
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               Do imod=1,Nmodels
                  X4D(ilon,ilat,imod,1:Nmon)=D_FF_T2m_Mon(ilon,ilat,1:Nmon,imod,isce)
               Enddo
            Enddo
         Enddo
         var_ln="Change in Monthly mean temperature for the period 2080-2099, RCP"//trim(Sce_nm); var_unit="C"
         Call wr_ncdf_4d_apd (Fname, X4D, "t_ff"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         ! FF mean Prec
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               Do imod=1,Nmodels
                  X4D(ilon,ilat,imod,1:Nmon)=D_FF_Pr_Mon(ilon,ilat,1:Nmon,imod,isce)
               Enddo
            Enddo
         Enddo
         var_ln="Change in Monthly mean rainfall amount for the period 2080-2099, RCP"//trim(Sce_nm); var_unit="mm/mon"
         Call wr_ncdf_4d_apd (Fname, X4D, "p_ff"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)

         ! P-value for differences between Sce & BL of Monthly mean (Re-write to Change file):
         var_ln="P-values for monthly temp. changes for period 2046-2065 and scenario of RCP"//trim(Sce_nm); var_unit="-"
         Call wr_ncdf_4d_apd (Fname,P_MF_T2m_Mon,"pval_t_mf"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         var_ln="P-values for monthly rainfall changes for period 2046-2065 and scenario of RCP"//trim(Sce_nm); var_unit="-"
         Call wr_ncdf_4d_apd (Fname,P_MF_Pr_Mon,"pval_p_mf"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         var_ln="P-values for monthly temp. changes for period 2080-2099 and scenario of RCP"//trim(Sce_nm); var_unit="-"
         Call wr_ncdf_4d_apd (Fname,P_FF_T2m_Mon,"pval_t_ff"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         var_ln="P-values for monthly rainfall changes for period 2080-2099 and scenario of RCP"//trim(Sce_nm); var_unit="-"
         Call wr_ncdf_4d_apd (Fname,P_FF_Pr_Mon,"pval_p_ff"//trim(Sce_nm),nlon,nlat,Nmodels,Nmon, &
                              var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
         
      Enddo
      ! BL: Differences between Obs & BL - Aug 2020:
      print*,"Differences between Obs & BL - Aug 2020"
      var_ln="Differences of Monthly mean temperature between BL and Obs for the period 1986-2005" 
      Call wr_ncdf_4d_apd (Fname, D_BL_T2m_Mon,"delt",nlon,nlat,Nmodels,Nmon, &
                              var_ln,"C", lon_nm, lat_nm, z_nm, t_nm,Miss)
      var_ln="Differences of Monthly preci between BL and Obs for the period 1986-2005"
      Call wr_ncdf_4d_apd (Fname, D_BL_Pr_Mon,"delp",nlon,nlat,Nmodels,Nmon, &
                              var_ln,"C", lon_nm, lat_nm, z_nm, t_nm,Miss)

      ! P-value for differences between BL & Obs of Monthly mean (Re-write to Change file):
      var_ln="P-values for differences between BL & Obs of Monthly Temperature"; var_unit="-"
      Call wr_ncdf_4d_apd (Fname, P_BL_T2m_Mon,"pval_t_bl", nlon,nlat,Nmodels,Nmon, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
      var_ln="P-values for differences between BL & Obs of Monthly Precipitation"; var_unit="-"
      Call wr_ncdf_4d_apd (Fname, P_BL_Pr_Mon,"pval_p_bl", nlon,nlat,Nmodels,Nmon, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)

      Deallocate (X3D, X4D, Xmod, Xmon)
      print*,"End of Differences between Obs & BL - Aug 2020"

      !Deallocate (X3D, X4D, Xmod, Xmon)

   End Subroutine Chars
   !---------------------------------------------------------------------------------
   Subroutine Read_T_R (Fname)
      Character (Len=*) :: Fname
      Real, Allocatable :: xlon(:), xlat(:), xtime(:), T(:,:,:), R(:,:,:)
      Integer :: nlon, nlat, nt

      Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt0, T2m, "t2m")
      Deallocate (xlon, xlat, xtime)
      Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt0, Pr, "pr")
      Deallocate (xlon, xlat, xtime)

      !Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, T, "t2m")
      !Deallocate (xlon, xlat, xtime)
      !Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, R, "pr")
      !Deallocate (xlon, xlat, xtime)

         !if (.not.Allocated(T2m)) Allocate (T2m(Nlon,Nlat,Nt))
         !if (.not.Allocated(Pr )) Allocate (Pr (Nlon,Nlat,Nt))
      !print*,"In file "//trim(Fname),": ",Nlon,Nlat,Nt; pause
      !T2m(:,:,:)=T(:,:,:); Pr(:,:,:)=R(:,:,:)

      !Deallocate (T,R)
      !Ntime = Nt
     !print*,"In file "//trim(Fname),": ",Nlon,Nlat,Nt; pause

   End Subroutine Read_T_R 
   !---------------------------------------------------------------------------------
   Subroutine Read_Land_Mask
      Integer :: nx, ny
      Real, Allocatable :: lonx(:), latx(:)

      Call rd_ncdf_2d (trim(LandMaskFile), xlon,xlat, nlon,nlat,Mask,"land") ! all=2..8=7 regs
      Call rd_ncdf_2d (trim(RegionFile), lonx, latx, nx, ny, Regs, "reg")
   End Subroutine Read_Land_Mask

      !---------------------------------------------------------------------------------
      Logical Function Empty_Ar1(X,N)
         Integer :: N,i
         Real :: X(N)
         Logical :: Empty
         Empty=.False.
         Do i=1,N
            if (X(i)==Miss.OR.IsNaN(X(i))) then
               Empty=.True.; Exit
            endif
         Enddo
         Empty_Ar1=Empty
      End Function Empty_Ar1
      !---------------------------------------------------------------------------------
      Real Function Ave(X, N, Miss)
        Integer :: N, i, icount
        Real :: X(N), Mean, Miss
        Mean=0.; icount=0
        Do i=1, N
           if (X(i)/=Miss) then
              icount=icount+1; Mean=Mean+X(i)
           endif
        Enddo
        if (icount>0) then
           Mean=Mean/real(icount)
        else
           Mean=Miss
        endif
        Ave=Mean
      End Function Ave
   !---------------------------------------------------------------------------------
   Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
   Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
   Include "/home/tanpv/Mylibs/Stdev_1.f90"
   !---------------------------------------------------------------------------------

   END  
   !---------------------------------------------------------------------------------
   Include "/home/tanpv/Mylibs/Student_Prob.f90"
   !---------------------------------------------------------------------------------