Use NetCDF

   Real, Parameter :: Miss=-99.0
   Integer, Parameter :: Nreg=7, Nmon=12, IMiss=-99, NLead=5
   Integer, Parameter :: Yr1=1982, Yr2=2010, Ntime0=(Yr2-Yr1+1)*Nmon, Ntime=Ntime0+NLead
   Integer, Allocatable :: time(:)

   Integer :: i, LenSt, sys
   Character (Len=255) :: Fname, t_unit, DirIn, Str, time_unit, Fname1, Fname2
   Character (Len=255) :: Dir
   Character (Len=10)  :: Model_Name
   Real, Allocatable :: lon(:), lat(:), Regs(:,:)
   !----------------------------------------------------------------------------------------------
   ! Input data:
   Real, Allocatable, dimension(:,:,:,:) :: Pr, T2m, Tx, Tm         ! (Nlon,Nlat,Nlead,Ntime)
   Real, Allocatable, dimension(:,:,:,:) :: h500, h850, olr, q2, rh2, slp, u200, u850, w700

   Real, Allocatable, dimension(:,:,:)   :: Pr_O, T2m_O, Tx_O, Tm_O ! (Nlon,Nlat,Ntime) 
   Real, Allocatable :: xlon(:), xlat(:), xtime(:)
   Real, Allocatable :: X4d(:,:,:,:), X3d(:,:,:) 
   Integer :: Nlon, Nlat, nlo, nla
   !----------------------------------------------------------------------------------------------
   ! Model Climatology:
   Real, Allocatable, dimension(:,:,:,:) :: Pr_M_cli, T2m_M_cli, Tx_M_cli, Tm_M_cli ! (Nlon,Nlat,Nlead,Nmon)
   Real, Allocatable, dimension(:,:,:)   :: Pr_O_cli, T2m_O_cli, Tx_O_cli, Tm_O_cli ! (Nlon,Nlat,Nmon)
   Real, Allocatable, dimension(:,:,:,:) :: Pr_BS, T2m_BS, Tx_BS, Tm_BS ! (Nlon,Nlat,Nlead,Nmon)
   !----------------------------------------------------------------------------------------------
   ! Differences between Model and Obs:
   Real, Allocatable, dimension(:,:,:,:) :: Pr_Df, T2_Df, Tx_Df, Tm_Df ! (Nlon,Nlat,Nlead,Ntime)
   !----------------------------------------------------------------------------------------------
   ! Verification Scores:
   Real, Allocatable, dimension(:,:,:,:) :: Pr_ME, Pr_MAE, Pr_RME, Pr_RMAE, Pr_COR, Pr_pval  ! (Nlon,Nlat,Nlead,Nmon)
   Real, Allocatable, dimension(:,:,:,:) :: T2_ME, T2_MAE, T2_COR, T2_pval           ! (Nlon,Nlat,Nlead,Nmon)
   Real, Allocatable, dimension(:,:,:,:) :: Tx_ME, Tx_MAE, Tx_COR, Tx_pval           ! (Nlon,Nlat,Nlead,Nmon)
   Real, Allocatable, dimension(:,:,:,:) :: Tm_ME, Tm_MAE, Tm_COR, Tm_pval           ! (Nlon,Nlat,Nlead,Nmon)
   !----------------------------------------------------------------------------------------------
   Real, Allocatable :: RMEreg(:,:,:,:), RMAEreg(:,:,:,:), CORreg(:,:,:,:)  ! (Nreg,Mon,LT,Mem)
   Real, Allocatable :: RMEwreg(:,:,:,:), RMAEwreg(:,:,:,:), CORwreg(:,:,:,:)  ! (Nreg,Mon,Nws,Mem)

   Integer :: ilon,ilat,it,imon
   Character (Len=1),Allocatable :: IRegs(:,:) 

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

   !Allocate (IRegs(Nlon,Nlat)); Where (Regs==Miss); IRegs=" "; elsewhere; IRegs="*"; End where
   !do ilat=Nlat,1,-1
   !   print"(200A1)",(IRegs(ilon,ilat),ilon=1,Nlon)
   !enddo
   !pause

   Call Read_OBS   ! Read obs data of Pr, T2m, Tx, Tm

   ! Read model data:
   Dir="/home/tanpv/Data/562_VNU/Rfc_VN/"
!   Model_Name="clWRF"
   Model_Name="RegCM"

   Call Read_Fcs (Dir,Model_Name)
   Call Read_Predictors (Dir,Model_Name)

   Call Mask_Data
   print*,"Writing data ..."
   Dir="/work/users/hmo/tanpv/From_HOME/562/Model_Data/"
   Call Write_Model_Obs (Dir,Model_Name)
   
   CONTAINS
   !----------------------------------------------------------------------------------------------
   Subroutine Write_Model_Obs (Dir,Model_Name)
      Integer :: imon, id
      Character (Len=*) :: Dir,Model_Name
      Character (Len=100) :: var_nm,t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm
      Real, Allocatable :: xtime(:)

      lon_nm="lon"; lat_nm="lat"; t_nm="mon"; t_unit="months since 1982-01-01 "
      Allocate (X3d(Nlon,Nlat,Ntime),xtime(Ntime))
      Do imon=1,Ntime; xtime(imon)=imon; Enddo

      Do ilead=1,NLead
         !write(Fname,"(4A,I1.1,3A)") trim(Dir),trim(Model_Name),"/","lead",ilead,"_",trim(Model_Name),"_VN.nc"
         write(Fname,"(2A,I1.1,3A)") trim(Dir),"lead",ilead,"_",trim(Model_Name),"_VN.nc"
         print*,trim(Fname)
         var_ln="Observed Monthly Prec"; var_unit="mm/month"
         Call wr_ncdf_3d(Fname,lon,lat,xtime,Nlon,Nlat,Ntime,Pr_O,"pr_o",t_unit,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Observed Monthly T2m"; var_unit="C"
         Call wr_ncdf_3d_apd (Fname,T2m_O,"t2_o",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Observed Monthly Tx"; var_unit="C"
         Call wr_ncdf_3d_apd (Fname,Tx_O,"tx_o",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Observed Monthly Tm"; var_unit="C"
         Call wr_ncdf_3d_apd (Fname,Tm_O,"tm_o",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         !
         var_ln="Model Monthly Prec"; var_unit="mm/month"
         X3d(:,:,:)=Pr(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"pr",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Model Monthly T2m"; var_unit="C"
         X3d(:,:,:)=T2m(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"t2",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Model Monthly Tx";
         X3d(:,:,:)=Tx(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"tx",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Model Monthly Tm";
         X3d(:,:,:)=Tm(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"tm",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         !
         var_ln="Model Geopotential height at 500 hPa (m)"; var_unit="m"
         X3d(:,:,:)=h500(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"h500",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Model Geopotential height at 850 hPa (m)"; var_unit="m"
         X3d(:,:,:)=h850(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"h850",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Model TOA OUTGOING LONG WAVE (W m-2)"; var_unit="w/m2"
         X3d(:,:,:)=olr(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"olr",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Model QV at 2 M (kg kg-1)"; var_unit="kg/kg"
         X3d(:,:,:)=q2(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"q2",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Model Relative Humidity at 2m (%)"; var_unit="%" 
         X3d(:,:,:)=rh2(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"rh2",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Model Sea Levelp Pressure (hPa)"; var_unit="hPa"
         X3d(:,:,:)=slp(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"slp",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Model x-wind component at 200 hPa (m s-1)"; var_unit="m/s"
         X3d(:,:,:)=u200(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"u200",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Model x-wind component at 850 hPa (m s-1)"; var_unit="m/s"
         X3d(:,:,:)=u850(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"u850",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
         var_ln="Model z-wind component at 700 hPa (m s-1)"; var_unit="m/s"
         X3d(:,:,:)=w700(:,:,ilead,:)
         Call wr_ncdf_3d_apd (Fname,X3d,"w700",Nlon,Nlat,Ntime,var_ln,var_unit,lon_nm,lat_nm,t_nm,Miss)
      Enddo
         ! h500, h850, olr, q2, rh2, slp, u200, u850, w700
   End Subroutine Write_Model_Obs 
   !----------------------------------------------------------------------------------------------
   Subroutine Mask_Data
      Integer :: ilon,ilat,ilead,it

      Do ilead=1,NLead
         Do ilon=1,Nlon
            Do ilat=1,Nlat
               if (Regs(ilon,ilat)==Miss) then
                  h500(ilon,ilat,ilead,:)=Miss; h850(ilon,ilat,ilead,:)=Miss; olr(ilon,ilat,ilead,:)=Miss
                  q2(ilon,ilat,ilead,:)=Miss;   rh2(ilon,ilat,ilead,:)=Miss;  slp(ilon,ilat,ilead,:)=Miss
                  u200(ilon,ilat,ilead,:)=Miss; u850(ilon,ilat,ilead,:)=Miss; w700(ilon,ilat,ilead,:)=Miss

                  Pr(ilon,ilat,ilead,:)=Miss;   T2m(ilon,ilat,ilead,:)=Miss;  Tx(ilon,ilat,ilead,:)=Miss; 
                  Tm(ilon,ilat,ilead,:)=Miss
               endif
            Enddo
         Enddo
      Enddo
   End Subroutine Mask_Data
   !----------------------------------------------------------------------------------------------
   Subroutine Read_Predictors (Dir,Model_Name)
      Character (Len=*) :: Dir,Model_Name
      Integer :: Nla, Nlo, Nt, Yr0, it1, it2, ilead
      Character (Len=255) Fname
! h500, h850, olr, q2, rh2, slp, u200, u850, w700

      Allocate (h500(Nlon,Nlat,Nlead,Ntime), h850(Nlon,Nlat,Nlead,Ntime), olr(Nlon,Nlat,Nlead,Ntime))
      Allocate (q2(Nlon,Nlat,Nlead,Ntime), rh2(Nlon,Nlat,Nlead,Ntime), slp(Nlon,Nlat,Nlead,Ntime))
      Allocate (u200(Nlon,Nlat,Nlead,Ntime), u850(Nlon,Nlat,Nlead,Ntime), w700(Nlon,Nlat,Nlead,Ntime))
      h500=Miss; h850=Miss; olr=Miss; q2=Miss; rh2=Miss; slp=Miss; u200=Miss; u850=Miss; w700=Miss;

      Do ilead=1,NLead
         write(Fname,"(4A,I1.1,3A)") trim(Dir),trim(Model_Name),"/","lead",ilead,"_",trim(Model_Name),"_baseline_pre_VN.nc"
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"h500")
         it1=ilead+1; it2=it1+Nt-1
         h500(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"h850")
         h850(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"olr")
         olr(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"q2")
         q2(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"rh2")
         rh2(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"slp")
         slp(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"u200")
         u200(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"u850")
         u850(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"w700")
         w700(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
      Enddo
   End Subroutine Read_Predictors
   !----------------------------------------------------------------------------------------------
   Subroutine Read_Fcs (Dir,Model_Name)  ! Read forecast data
      Character (Len=*) :: Dir,Model_Name
      Integer :: Nla, Nlo, Nt, Yr0, it1, it2, ilead
      Character (Len=255) Fname

      !Dir="/home/tanpv/Data/562_VNU/Rfc_VN/"
      Allocate (Pr(Nlon,Nlat,Nlead,Ntime), T2m(Nlon,Nlat,Nlead,Ntime))
      Allocate (Tx(Nlon,Nlat,Nlead,Ntime), Tm(Nlon,Nlat,Nlead,Ntime))         ! (Nlon,Nlat,Nlead,Ntime)
      Pr=Miss; T2m=Miss; Tx=Miss; Tm=Miss
      Do ilead=1,NLead
         write(Fname,"(4A,I1.1,3A)") trim(Dir),trim(Model_Name),"/","lead",ilead,"_",trim(Model_Name),"_baseline_VN.nc"
         !print*,trim(Fname)
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"pr")
         it1=ilead+1; it2=it1+Nt-1 ! Missing data will be filled at "first ilead values" and "last NLead-ilead values"
                                   ! Nt = 348, Ntime = 353
         Pr(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tas")
         T2m(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tasmax")
         Tx(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
         Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tasmin")
         Tm(:,:,ilead,it1:it2)=X3d(:,:,1:Nt)
         Deallocate (xlon,xlat,xtime,X3d)
      Enddo
   End Subroutine Read_Fcs
   !----------------------------------------------------------------------------------------------
   Subroutine Read_OBS  ! Read observed gridded data of Pr, T2m, Tx, Tm
      Integer :: j, it, ilon, ilat, Nla, Nlo, Nt, Yr0, it1, it2, imon
      Character (Len=255) Fname

      Allocate (Pr_O(Nlon,Nlat,Ntime), T2m_O(Nlon,Nlat,Ntime), Tx_O(Nlon,Nlat,Ntime), Tm_O(Nlon,Nlat,Ntime))

      Yr0=1980;              ! Rain data are available from Jan 1980
      it1=(Yr1-Yr0)*Nmon+1   ! Lay tu 1/1982 
      it2=it1+Ntime-1        ! Lay den 5/2011
!print*,it1,it2,it2-it1+1,Yr0+(it1-1)/Nmon; pause

      Fname="/home/tanpv/Data/562_VNU/OBS_VN_Grd/R.1980.2018.mm.nc"
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"rain")
      Pr_O(:,:,1:Ntime)=X3d(:,:,it1:it2)
      Deallocate (xlon,xlat,xtime,X3d)

      Fname="/home/tanpv/Data/562_VNU/OBS_VN_Grd/T2m.1961.2018.mm.nc"
      Yr0=1961;              ! Temp data are available from Jan 1961
      it1=(Yr1-Yr0)*Nmon+1; 
      it2=it1+Ntime-1
!print*,it1,it2,it2-it1+1,Yr0+(it1-1)/Nmon; pause
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"t2m")
      T2m_O(:,:,1:Ntime)=X3d(:,:,it1:it2)
      Deallocate (xlon,xlat,xtime,X3d)

      Fname="/home/tanpv/Data/562_VNU/OBS_VN_Grd/Tx.1961.2018.mm.nc"
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tx")
      Tx_O(:,:,1:Ntime)=X3d(:,:,it1:it2)
      Deallocate (xlon,xlat,xtime,X3d)

      Fname="/home/tanpv/Data/562_VNU/OBS_VN_Grd/Tm.1961.2018.mm.nc"
      Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"tm")
      Tm_O(:,:,1:Ntime)=X3d(:,:,it1:it2)
      Deallocate (xlon,xlat,xtime,X3d)
   End Subroutine Read_OBS
   !----------------------------------------------------------------------------------------------
   Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
   Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
   Include "/home/tanpv/Mylibs/Ave.f90"
   Include "/home/tanpv/Mylibs/DayOfMon.f90"
   END
   !---------------------------------------------------------------------------