Use NetCDF
Real, Allocatable :: XTime(:), T2m(:,:,:), Pr(:,:,:), xLon(:), xLat(:), Mask(:,:)
Real, Allocatable :: lon(:), lat(:), dat(:), AWC(:,:)
Integer :: Nlon, Nlat, Nt, it, ilon, ilat, Nst, ist
Real, Parameter :: Miss=-99.0, Miss_aph=-99.9
Character (Len=255) :: R_Fname
Logical :: status

  !Call Proc_awc      ! Global AWC in nc 

  !Call Create_awc_1  ! Get AWC for SEA from Global AWC 0.25

  !Call Create_awc_VN ! Create AWC form VN from both FAO and VN data 0.25
 
  !Call AWC_VN_HiRes   ! HiRes

  Call Create_awc_VN_From_Global  ! VN from FAO (Global)

CONTAINS
!------------------------------------------------------------------------------
Subroutine AWC_VN_HiRes
   Real, Allocatable :: lon(:), lat(:), dat(:), xlon(:), xlat(:), AWC(:,:), Mask(:,:)
   Integer :: Nlon, Nlat, Nt, it, ilon, ilat, Nst, ist
   integer, parameter :: ncols = 1208, nrows = 1976
   real, parameter :: xcor = 100.6063602, ycor = 7.804953819
   real, parameter :: cellsiz = 0.008333333, iflag = -9999, oflag = -99
   real :: indat(ncols,nrows), x(ncols), y(nrows)
   ! -----------------------------------------------------------<o
   Open (1, file = "/home/tanpv/Data/AWC/AWC_VN.txt", Status="Old")
      read (1,*) ; read (1,*) ; read (1,*) ; read (1,*)
      read (1,*) ; read (1,*) ; read (1,*)
      do irow = 1, nrows
         itmp = nrows - irow + 1
         read (1,*) (indat(icol,itmp),icol = 1, ncols)
      enddo
   close (1)
   Where (indat==-9999) indat=Miss
   do icol = 1, ncols
      x(icol) = xcor + (icol - 1) * cellsiz
   enddo
   do irow = 1, nrows
      y(irow) = ycor + (irow - 1) * cellsiz
   enddo

   Call wr_ncdf_2d ("AWC_VN_HiRes.nc",x,y,ncols,nrows,indat,"awc","Available Water Content from FAO combined with VN","inch", "lon", "lat", Miss)


End Subroutine AWC_VN_HiRes
!------------------------------------------------------------------------------
Subroutine Create_Land_mask

  R_Fname="/home/tanpv/Data/CORDEX_SEACLID/R_T_mon_1970_2100/pr_APH_historical_SRF.nc"
  Call rd_ncdf_3d (trim(R_Fname), xlon, xlat, XTime, Nlon, Nlat, Nt, Pr, "pr")

  Nst=Nlon*Nlat
  ! print*, lon, lat
  ! print*, Pr
  Allocate (Mask(Nlon, Nlat),lon(Nst),lat(Nst),dat(Nst), AWC(Nlon, Nlat))
  ist=0
  Do ilat=1, Nlat
     Do ilon=1, Nlon
        Do it=50, 300
           if (Pr(ilon,ilat,it) /= Miss_aph) then
              Mask(ilon,ilat) = 1.0
           else
              Mask(ilon,ilat) = Miss
           endif
        Enddo
        ist=ist+1
        lon(ist)=xlon(ilon); lat(ist)=xlat(ilat);
     Enddo
  Enddo
  !Call wr_ncdf_2d ("Land_mask_sea.nc",xlon,xlat,nlon,nlat,Mask,"land", "Land mask of SEA", "-", "lon", "lat", Miss)
End Subroutine Create_Land_mask
!------------------------------------------------------------------------------
Subroutine Create_awc_VN
  Real, Allocatable :: lon(:), lat(:), dat(:), xlon(:), xlat(:), AWC(:,:), Mask(:,:)
  Integer :: Nlon, Nlat, Nt, it, ilon, ilat, Nst, ist

  Call rd_ncdf_2d ("Land_Mask_VN.nc",xlon,xlat,nlon,nlat,Mask,"vn")

  Nst=Nlon*Nlat
  Allocate (dat(Nst), lon(Nst), lat(Nst))
  ist=0
  Do ilat=1, Nlat
     Do ilon=1, Nlon
        ist=ist+1; lon(ist)=xlon(ilon); lat(ist)=xlat(ilat)
     Enddo
  Enddo
  Call putAWC_VN(dat,lon,lat,nst)
  Allocate (AWC(Nlon,Nlat))
  ist=0
  Do ilat=1, Nlat
     Do ilon=1, Nlon
        ist=ist+1
        if (Mask(ilon,ilat)>0) then !.and.dat(ist)>=0) then
           AWC(ilon,ilat)=dat(ist)
           if (AWC(ilon,ilat)<0) AWC(ilon,ilat)=6.0
        else
           AWC(ilon,ilat)=Miss
        endif
     Enddo
  Enddo
  Call wr_ncdf_2d ("Avai_Water_Content_VN.nc",xlon,xlat,nlon,nlat,AWC,"awc","Available Water Content from FAO combined with VN","inch", "lon", "lat", Miss)
End Subroutine Create_awc_VN
!------------------------------------------------------------------------------
Subroutine Create_awc_VN_From_Global
  Real, Allocatable :: lon(:), lat(:), dat(:), xlon(:), xlat(:), AWC(:,:), Mask(:,:)
  Integer :: Nlon, Nlat, Nt, it, ilon, ilat, Nst, ist

  Call rd_ncdf_2d ("Land_Mask_VN.nc",xlon,xlat,nlon,nlat,Mask,"vn")

  Nst=Nlon*Nlat
  Allocate (dat(Nst), lon(Nst), lat(Nst))
  ist=0
  Do ilat=1, Nlat
     Do ilon=1, Nlon
        ist=ist+1; lon(ist)=xlon(ilon); lat(ist)=xlat(ilat)
     Enddo
  Enddo
  Call putAWC(dat,lon,lat,nst)
  Allocate (AWC(Nlon,Nlat))
  ist=0
  Do ilat=1, Nlat
     Do ilon=1, Nlon
        ist=ist+1
        if (Mask(ilon,ilat)>0) then !.and.dat(ist)>=0) then
           AWC(ilon,ilat)=dat(ist)
           if (AWC(ilon,ilat)<0) AWC(ilon,ilat)=6.0
        else
           AWC(ilon,ilat)=Miss
        endif
     Enddo
  Enddo
  Call wr_ncdf_2d ("AWC_VN_FAO.nc",xlon,xlat,nlon,nlat,AWC,"awc","Available Water Content from FAO","inch", "lon", "lat", Miss)
End Subroutine Create_awc_VN_From_Global
!------------------------------------------------------------------------------
SUBROUTINE putAWC_VN(dat,lon,lat,nst)

  !ncols   1208    ! Number of columns
  !nrows   1976    ! Number of rows
  !xllcorner       100.6063602     ! Longitude of Lower left corner
  !yllcorner       7.804953819     ! Latitude of Lower left corner
  !cellsize        0.008333333     ! Do phan giai

   integer :: nst, ist
   real :: dat(nst), lon(nst), lat(nst)
   !integer, parameter :: ncols = 1099, nrows = 549
   integer, parameter :: ncols = 1208, nrows = 1976
   integer :: icol, irow, itmp
   !character (len = *), parameter :: fileAWC = "/work/users/thanhnx/Func_Sub/AWC_Global.asc"
   !real, parameter :: xcor = -180.0, ycor = -89.9056530001
   !real, parameter :: cellsiz = 0.327697, iflag = -9999, oflag = -99
   real, parameter :: xcor = 100.6063602, ycor = 7.804953819
   real, parameter :: cellsiz = 0.008333333, iflag = -9999, oflag = -99
   real :: indat(ncols,nrows), x(ncols), y(nrows)
   ! <-~~~=|==-----------------------------------------------------------<o
   ! <-~~~=|== Read asc AWC file
   !open (1,file = fileAWC,status = "old")
   Open (1, file = "/home/tanpv/Data/AWC/AWC_VN.txt", Status="Old")
      ! <-~~~=| Skip reading 7 lines
      read (1,*) ; read (1,*) ; read (1,*) ; read (1,*)
      read (1,*) ; read (1,*) ; read (1,*)
      ! <-~~~=|
      do irow = 1, nrows
         itmp = nrows - irow + 1
         read (1,*) (indat(icol,itmp),icol = 1, ncols)
      enddo
   close (1)
   close (1)
   do icol = 1, ncols
      x(icol) = xcor + (icol - 1) * cellsiz
   enddo
   do irow = 1, nrows
      y(irow) = ycor + (irow - 1) * cellsiz
   enddo
   ! <-~~~=| Compute for new grid
   do ist = 1, nst
      icol = nrtpos(x,lon(ist),ncols,oflag)
      irow = nrtpos(y,lat(ist),nrows,oflag)
      dat(ist) = indat(icol,irow)*0.0393701
      if (dat(ist) == iflag) then
         dat(ist) = oflag   ! Doi -99.9 --> 99.0 
      endif
   !   print*,ist, dat(ist)
   enddo
return
END SUBROUTINE putAWC_VN
!------------------------------------------------------------------------------
Subroutine Create_awc_1
  Real, Allocatable :: lon(:), lat(:), dat(:), xlon(:), xlat(:), AWC(:,:), Mask(:,:)
  Integer :: Nlon, Nlat, Nt, it, ilon, ilat, Nst, ist

  Call rd_ncdf_2d ("Land_mask_sea.nc",xlon,xlat,nlon,nlat,Mask,"land")

  Nst=Nlon*Nlat
  Allocate (dat(Nst), lon(Nst), lat(Nst))
  ist=0
  Do ilat=1, Nlat
     Do ilon=1, Nlon
        ist=ist+1; lon(ist)=xlon(ilon); lat(ist)=xlat(ilat)
     Enddo
  Enddo
  Call putAWC(dat,lon,lat,nst)
  Allocate (AWC(Nlon,Nlat))
  ist=0
  Do ilat=1, Nlat
     Do ilon=1, Nlon
        ist=ist+1
        if (Mask(ilon,ilat)>0) then !.and.dat(ist)>=0) then
           AWC(ilon,ilat)=dat(ist)
           if (AWC(ilon,ilat)<0) AWC(ilon,ilat)=6.0
        else
           AWC(ilon,ilat)=Miss
        endif
     Enddo
  Enddo
  Call wr_ncdf_2d ("Avai_Water_Content_1.nc",xlon,xlat,nlon,nlat,AWC,"awc","Available Water Content from FAO","inch", "lon", "lat", Miss)
End Subroutine Create_awc_1

!------------------------------------------------------------------------------
Subroutine Proc_awc
   Integer, Parameter :: ncols=1099, nrows=549
   Real*8, Parameter :: xllcorner=-180.0, yllcorner=-89.9056530001, cellsize=0.327697, Miss_dat=-9999.0
   Real*8 :: dat(ncols,nrows), dxy
   Character (Len=100) :: Fname
   Integer :: ilon,ilat, Nlon,Nlat
   Real, Allocatable :: awc(:,:), lon(:), lat(:)

   Fname="/home/tanpv/Data/AWC/AWC_Global.asc"
   Open (1, file = Fname, Status="Old")
   Do ilat=1,7; read(1,*); Enddo   ! Skip 7 lines
   Do ilat=1, nrows
      Read (1,*) (dat(ilon,ilat),ilon=1, ncols) 
   Enddo
   Nlat=nrows; Nlon=ncols-1; dxy=360.0/(Nlon)
   Allocate (lon(Nlon),lat(Nlat),awc(Nlon,Nlat))
   Do ilon=1, Nlon
      lon(ilon)=-180.0+(ilon-1)*dxy
      !print*,lon(ilon)
   Enddo
   !pause
   Do ilat=1, Nlat
      !lat(ilat)=-90.0-dxy/2+(ilat-1)*dxy
      lat(ilat)=yllcorner-dxy/2+(ilat-1)*dxy
      !print*,lat(ilat)
   Enddo
   !pause
   Do ilat=1, Nlat
      Do ilon=1, Nlon
         awc(ilon,ilat)=dat(ilon,Nlat-ilat+1)
      Enddo
   Enddo
   Where (awc==Miss_dat) awc=Miss
   Fname="awc_global.nc"
   Call wr_ncdf_2d (Fname,lon,lat,nlon,nlat,awc,"awc","Available Water Content from FAO","mm", "lon", "lat", Miss)

   print*,"Proc_awc terminated."
End Subroutine Proc_awc
!------------------------------------------------------------------------------
SUBROUTINE putAWC(dat,lon,lat,nst)
   integer :: nst, ist
   real :: dat(nst), lon(nst), lat(nst)
   integer, parameter :: ncols = 1099, nrows = 549
   integer :: icol, irow, itmp
   !character (len = *), parameter :: fileAWC = "/work/users/thanhnx/Func_Sub/AWC_Global.asc"
   real, parameter :: xcor = -180.0, ycor = -89.9056530001
   real, parameter :: cellsiz = 0.327697, iflag = -9999, oflag = -99
   real :: indat(ncols,nrows), x(ncols), y(nrows)
   ! <-~~~=|==-----------------------------------------------------------<o
   ! <-~~~=|== Read asc AWC file
   !open (1,file = fileAWC,status = "old")
   Open (1, file = "/home/tanpv/Data/AWC/AWC_Global.asc", Status="Old")
      ! <-~~~=| Skip reading 7 lines
      read (1,*) ; read (1,*) ; read (1,*) ; read (1,*)
      read (1,*) ; read (1,*) ; read (1,*)
      ! <-~~~=|
      do irow = 1, nrows
         itmp = nrows - irow + 1
         read (1,*) (indat(icol,itmp),icol = 1, ncols)
      enddo
   close (1)
   close (1)
   do icol = 1, ncols
      x(icol) = xcor + (icol - 1) * cellsiz
   enddo
   do irow = 1, nrows
      y(irow) = ycor + (irow - 1) * cellsiz
   enddo
   ! <-~~~=| Compute for new grid
   do ist = 1, nst
      icol = nrtpos(x,lon(ist),ncols,oflag)
      irow = nrtpos(y,lat(ist),nrows,oflag)
      dat(ist) = indat(icol,irow)*0.0393701
      if (dat(ist) == iflag) then
         dat(ist) = oflag   ! Doi -99.9 --> 99.0 
      endif
   !   print*,ist, dat(ist)
   enddo
return
END SUBROUTINE putAWC
!--------------------------------------------------------------------
integer FUNCTION nrtpos(x,y,n,flag)
   integer :: n, i, itmp
   real :: flag, x(n), y, minx

   minx = 999999
   itmp = 0
   do i = 1, n
      if (x(i) /= flag .and. y /= flag) then
         if (abs(x(i) - y) < minx) then
            minx = abs(x(i) - y)
            itmp = i
         endif
      endif
   enddo

   nrtpos = itmp
   return
END FUNCTION nrtpos
!--------------------------------------------------------------------

 Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
 Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
 END