      Use NetCDF
      Real, Parameter :: Miss=-99.
      Integer, Parameter :: Nmon=12, Nreg=8
      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(:,:), Mask_all(:,:)                  ! (Nlon,Nlat)

      ! Input Data (Reading from files): 
      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 :: 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 :: 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 :: 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 :: 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_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)        
      
      ! Changes of Monthly mean:
      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)  


      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
      Character (Len=100) :: Data_Dir, Model_nm(Nmodels), Model_Obs, Out_Dir
      Integer :: YSta(NPeriod), YEnd(NPeriod)
      Character (Len=5) :: Sce_nm(NSce), Period(NPeriod)
      !---------------------------------------------------------------------------------
      Namelist /Infor/ Data_Dir, Period, YSta, YEnd, Nyear, Model_nm, Sce_nm, LandMaskFile, Out_Dir
      Integer :: ilon,ilat,it,imon,iper,imod,isce
      !---------------------------------------------------------------------------------

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

      Call Read_Land_Mask
      Ntime = Nyear*Nmon

      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: 6 models + 1 Ens
            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
            Enddo
         else
            Do isce=1,NSce
               Do imod=1, Nmodels
                  Fname=trim(Data_Dir)//"SCE_"//trim(Period(iper))//"_"//trim(Sce_nm(isce))//"_"//trim(Model_nm(imod))//".nc"
                  Call Read_T_R (Fname)
                  if (iper==2) then   ! MF
                     Do ilon=1,Nlon
                        Do ilat=1,Nlat
                           MF_T2m(ilon,ilat,1:Ntime,imod,isce)=T2m(ilon,ilat,1:Ntime)
                           MF_Pr (ilon,ilat,1:Ntime,imod,isce)=Pr (ilon,ilat,1:Ntime)
                        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,1:Ntime)
                           FF_Pr (ilon,ilat,1:Ntime,imod,isce)=Pr (ilon,ilat,1:Ntime)
                        Enddo
                     Enddo
                     !FF_T2m(:,:,:,imod,isce)=T2m(:,:,:); FF_Pr(:,:,:,imod,isce)=Pr(:,:,:)
                  endif
               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*,FF_Pr(ilon,ilat,1:Ntime,Nmodels,NSce)
      !         pause
      !      endif
      !   Enddo
      !Enddo

      Call Chars
      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
      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 (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 (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))

      BL_T2m_RMon=Miss; BL_Pr_RMon=Miss; MF_T2m_Rmon=Miss; MF_Pr_RMon=Miss; FF_T2m_Rmon=Miss; FF_Pr_RMon=Miss

      ! 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
      ! 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 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"

      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)

      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 (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
   !---------------------------------------------------------------------------------
   Subroutine For_Can_Tho
      Integer :: ilon,ilat,imon,imod,isce, it, iyr, imon1,it1, it2
      Real :: X(Nlon,Nlat), Xmod(Nmodels), Xmon(Nmon), Xyrs(Nyear)
      Real, Parameter :: KDo=105.767, VDo=10.033             ! Can Tho Station

      Real, Allocatable :: BL_T(:,:), BL_P(:,:)              ! (Nmodels,Ntime)
      Real, Allocatable :: MF_T(:,:,:), MF_P(:,:,:)          ! (Nmodels,NSce,Ntime)
      Real, Allocatable :: FF_T(:,:,:), FF_P(:,:,:)          ! (Nmodels,NSce,Ntime)

      Real, Allocatable :: BL_T_Mon(:,:), BL_P_Mon(:,:)      ! (Nmodels,Nmon)
      Real, Allocatable :: MF_T_Mon(:,:,:), MF_P_Mon(:,:,:)  ! (Nmodels,NSce,Nmon)
      Real, Allocatable :: FF_T_Mon(:,:,:), FF_P_Mon(:,:,:)  ! (Nmodels,NSce,Nmon)

      Real, Allocatable :: BL_T_Yrs(:,:), BL_P_Yrs(:,:)      ! (Nmodels,Nyear)
      Real, Allocatable :: MF_T_Yrs(:,:,:), MF_P_Yrs(:,:,:)  ! (Nmodels,NSce,Nyear)
      Real, Allocatable :: FF_T_Yrs(:,:,:), FF_P_Yrs(:,:,:)  ! (Nmodels,NSce,Nyear)

      Real, Allocatable :: X1(:),X2(:),Y1(:),Y2(:), X2Dm(:,:), X2Dy(:,:)
      Character (Len=100) :: Fname, Sce_nm, t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm

      Allocate (BL_T(Nmodels,Ntime),BL_P(Nmodels,Ntime))
      Allocate (MF_T(Nmodels,NSce,Ntime),MF_P(Nmodels,NSce,Ntime))
      Allocate (FF_T(Nmodels,NSce,Ntime),FF_P(Nmodels,NSce,Ntime))

      Allocate (BL_T_Mon(Nmodels,Nmon),BL_P_Mon(Nmodels,Nmon))
      Allocate (MF_T_Mon(Nmodels,NSce,Nmon),MF_P_Mon(Nmodels,NSce,Nmon))
      Allocate (FF_T_Mon(Nmodels,NSce,Nmon),FF_P_Mon(Nmodels,NSce,Nmon))

      Allocate (BL_T_Yrs(Nmodels,Nyear),BL_P_Yrs(Nmodels,Nyear))
      Allocate (MF_T_Yrs(Nmodels,NSce,Nyear),MF_P_Yrs(Nmodels,NSce,Nyear))
      Allocate (FF_T_Yrs(Nmodels,NSce,Nyear),FF_P_Yrs(Nmodels,NSce,Nyear))

      ! Input Data (Interpolating to station location from origin data):
      Do imod=1, Nmodels
         Do it=1, Ntime
            X(:,:)=BL_T2m(:,:,it,imod); BL_T(imod,it)=Get_Point_Value(X, KDo, VDo)
            X(:,:)=BL_Pr (:,:,it,imod); BL_P(imod,it)=Get_Point_Value(X, KDo, VDo)
         Enddo
      Enddo
      Do isce=1, NSce
         Do imod=1, Nmodels
            Do it=1, Ntime
               X(:,:)=MF_T2m(:,:,it,imod,isce); MF_T(imod,isce,it)=Get_Point_Value(X, KDo, VDo)
               X(:,:)=MF_Pr (:,:,it,imod,isce); MF_P(imod,isce,it)=Get_Point_Value(X, KDo, VDo)
               X(:,:)=FF_T2m(:,:,it,imod,isce); FF_T(imod,isce,it)=Get_Point_Value(X, KDo, VDo)
               X(:,:)=FF_Pr (:,:,it,imod,isce); FF_P(imod,isce,it)=Get_Point_Value(X, KDo, VDo)
            Enddo
         Enddo
      Enddo

      ! Monthly mean:
      Allocate (X2(Nyear),Y2(Nyear),X1(Nyear),Y1(Nyear))
      ! 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; X1(iyr)=BL_T(imod,it); X2(iyr)=BL_P(imod,it)
               endif
            Enddo
            BL_T_Mon(imod,imon)=Ave(X1,iyr,Miss); BL_P_Mon(imod,imon)=Ave(X2,iyr,Miss)
         Enddo
      Enddo
      ! Future:
      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; 
                     X1(iyr)=MF_T(imod,isce,it); X2(iyr)=MF_P(imod,isce,it)
                     Y1(iyr)=FF_T(imod,isce,it); Y2(iyr)=FF_P(imod,isce,it)
                  endif
               Enddo
               MF_T_Mon(imod,isce,imon)=Ave(X1,iyr,Miss); MF_P_Mon(imod,isce,imon)=Ave(X2,iyr,Miss)
               FF_T_Mon(imod,isce,imon)=Ave(Y1,iyr,Miss); FF_P_Mon(imod,isce,imon)=Ave(Y2,iyr,Miss)
            Enddo
         Enddo
      Enddo
      Deallocate (X1,X2,Y1,Y2)
      ! Yearly:
      Allocate (X1(Nmon),X2(Nmon),Y1(Nmon),Y2(Nmon))
      Do imod=1, Nmodels
         ! Baseline:
         Do iyr=1, Nyear
            it1=(iyr-1)*Nmon+1; it2=iyr*Nmon
            X1(1:Nmon)=BL_T(imod,it1:it2); X2(1:Nmon)=BL_P(imod,it1:it2)
            BL_T_Yrs(imod,iyr)=Ave(X1,Nmon,Miss); BL_P_Yrs(imod,iyr)=Ave(X2,Nmon,Miss)*Nmon
         Enddo
         ! Future:
         Do isce=1, NSce
            Do iyr=1, Nyear
               it1=(iyr-1)*Nmon+1; it2=iyr*Nmon
               X1(1:Nmon)=MF_T(imod,isce,it1:it2); X2(1:Nmon)=MF_P(imod,isce,it1:it2)
               Y1(1:Nmon)=FF_T(imod,isce,it1:it2); Y2(1:Nmon)=FF_P(imod,isce,it1:it2)
               MF_T_Yrs(imod,isce,iyr)=Ave(X1,Nmon,Miss); MF_P_Yrs(imod,isce,iyr)=Ave(X2,Nmon,Miss)*Nmon
               FF_T_Yrs(imod,isce,iyr)=Ave(Y1,Nmon,Miss); FF_P_Yrs(imod,isce,iyr)=Ave(Y2,Nmon,Miss)*Nmon
            Enddo
         Enddo
      Enddo
      Deallocate (X1,X2,Y1,Y2)
      ! Wite to files:
      Allocate (X2Dm(Nmodels,Nmon), X2Dy(Nmodels,Nyear))
      Do imon=1,Nmon; Xmon(imon)=imon; Enddo
      Do imod=1,Nmodels; Xmod(imod)=imod; Enddo
      Do iyr=1, Nyear; Xyrs(iyr)=iyr; Enddo
      
      Fname=trim(Out_Dir)//"CT_Monthly_Mean.nc"
      var_ln="Monthly mean 2m-temperature for Can Tho"; var_unit="C"; lon_nm="Model"; lat_nm="Mon"
      Call wr_ncdf_2d (Fname,Xmod,Xmon,Nmodels,Nmon,BL_T_Mon,"t_bl", var_ln, var_unit, lon_nm, lat_nm, Miss)
      var_ln="Monthly mean rainfall for Can Tho"; var_unit="mm/mon"
      Call wr_ncdf_2d_apd (Fname,BL_P_Mon,"p_bl",Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, Miss)
      Do isce=1, Nsce
         if (isce==1) then; Sce_nm="45"; else; Sce_nm="85"; endif
         X2Dm(1:Nmodels,1:Nmon)=MF_T_Mon(1:Nmodels,isce,1:Nmon)
         var_ln="Monthly mean 2m-temperature for Can Tho"; var_unit="C"
         Call wr_ncdf_2d_apd (Fname,X2Dm,"t_mf"//trim(Sce_nm),Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, Miss)
         X2Dm(1:Nmodels,1:Nmon)=MF_P_Mon(1:Nmodels,isce,1:Nmon)
         var_ln="Monthly mean rainfall for Can Tho"; var_unit="mm/mon"
         Call wr_ncdf_2d_apd (Fname,X2Dm,"p_mf"//trim(Sce_nm),Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, Miss)
         X2Dm(1:Nmodels,1:Nmon)=FF_T_Mon(1:Nmodels,isce,1:Nmon)
         var_ln="Monthly mean 2m-temperature for Can Tho"; var_unit="C"
         Call wr_ncdf_2d_apd (Fname,X2Dm,"t_ff"//trim(Sce_nm),Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, Miss)
         X2Dm(1:Nmodels,1:Nmon)=FF_P_Mon(1:Nmodels,isce,1:Nmon)
         var_ln="Monthly mean rainfall for Can Tho"; var_unit="mm/mon"
         Call wr_ncdf_2d_apd (Fname,X2Dm,"p_ff"//trim(Sce_nm),Nmodels,Nmon, var_ln, var_unit, lon_nm, lat_nm, Miss)
      Enddo
 
      Fname=trim(Out_Dir)//"CT_Annual_Mean.nc"
      var_ln="Annual mean 2m-temperature for Can Tho"; var_unit="C"; lon_nm="Model"; lat_nm="year"
      Call wr_ncdf_2d (Fname,Xmod,Xyrs,Nmodels,Nyear,BL_T_Yrs,"t_bl", var_ln, var_unit, lon_nm, lat_nm, Miss)
      var_ln="Annual mean rainfall for Can Tho"; var_unit="mm/year"
      Call wr_ncdf_2d_apd (Fname,BL_P_Yrs,"p_bl",Nmodels,Nyear, var_ln, var_unit, lon_nm, lat_nm, Miss)
      Do isce=1, Nsce
         if (isce==1) then; Sce_nm="45"; else; Sce_nm="85"; endif
         X2Dy(1:Nmodels,1:Nyear)=MF_T_Yrs(1:Nmodels,isce,1:Nyear)
         var_ln="Annual mean 2m-temperature for Can Tho"; var_unit="C"
         Call wr_ncdf_2d_apd (Fname,X2Dy,"t_mf"//trim(Sce_nm),Nmodels,Nyear, var_ln, var_unit, lon_nm, lat_nm, Miss)
         X2Dy(1:Nmodels,1:Nyear)=MF_P_Yrs(1:Nmodels,isce,1:Nyear)
         var_ln="Annual mean rainfall for Can Tho"; var_unit="mm/year"
         Call wr_ncdf_2d_apd (Fname,X2Dy,"p_mf"//trim(Sce_nm),Nmodels,Nyear, var_ln, var_unit, lon_nm, lat_nm, Miss)
         X2Dy(1:Nmodels,1:Nyear)=FF_T_Yrs(1:Nmodels,isce,1:Nyear)
         var_ln="Annual mean 2m-temperature for Can Tho"; var_unit="C"
         Call wr_ncdf_2d_apd (Fname,X2Dy,"t_ff"//trim(Sce_nm),Nmodels,Nyear, var_ln, var_unit, lon_nm, lat_nm, Miss)
         X2Dy(1:Nmodels,1:Nyear)=FF_P_Yrs(1:Nmodels,isce,1:Nyear)
         var_ln="Annual mean rainfall for Can Tho"; var_unit="mm/year"
         Call wr_ncdf_2d_apd (Fname,X2Dy,"p_ff"//trim(Sce_nm),Nmodels,Nyear, var_ln, var_unit, lon_nm, lat_nm, Miss)
      Enddo
      Deallocate (X2Dm,X2Dy)

   End Subroutine For_Can_Tho
   !---------------------------------------------------------------------------------
   Subroutine Chars
      Integer :: ilon,ilat,it,imon,imon1,iper,imod,isce
      Integer :: iyr
      Real, Allocatable :: X(:),Y(:),X1(:),Y1(:), X3D(:,:,:), X4D(:,:,:,:), Xmod(:), Xmon(:), XX(:)
      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 (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 (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_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))

      ! Changes of Monthly mean:
      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))

      Allocate (X(Nyear),Y(Nyear),X1(Nyear),Y1(Nyear),Xmon(Nmon), XX(Ntime))

      BL_T2m_Mean=Miss; BL_Pr_Mean=Miss; MF_T2m_Mean=Miss; MF_Pr_Mean=Miss; FF_T2m_Mean=Miss; FF_Pr_Mean=Miss
      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_MF_T2m_Mon=Miss; D_MF_Pr_Mon=Miss; D_FF_T2m_Mon=Miss; D_FF_Pr_Mon=Miss
      D_MF_T2m_Ann=Miss; D_MF_Pr_Ann=Miss; D_FF_T2m_Ann=Miss; D_FF_Pr_Ann=Miss

      Do ilon=1, Nlon
         Do ilat=1, Nlat
            if (Mask(ilon,ilat)>0) then
               !print*,ilon,ilat,Nmodels,Nmon,Ntime
               ! 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)
                           !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)
                  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
               Enddo
               ! 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; 
                              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)

                        ! 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
                  Enddo
               Enddo
            endif
         Enddo
      Enddo
      Deallocate (X,Y,X1,Y1,Xmon)

      ! 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"
      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)
 
      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)
         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)
         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)
         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)
         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) 
      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
      !========================================================
      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)

      !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)
      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
      ! 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)
         
      Enddo

      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, 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
      Call rd_ncdf_2d (trim(LandMaskFile), xlon,xlat, nlon,nlat,Mask,"all") ! all=2..8=7 regs
      Allocate (Mask_all(nlon,nlat))
      Mask_all=Mask
   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"
   !---------------------------------------------------------------------------------

   END  
