Use NetCDF
      Integer, parameter :: Nreg=7
      Real, parameter :: Miss=-99
      Real, Allocatable, Dimension(:,:) :: VN,B1,B2,B3,B4,N1,N2,N3
      Real, Allocatable, Dimension(:)   :: lon, lat

      Integer :: Nlon,Nlat,ilon,ilat,ireg

      Character (Len=255) :: Fname

      Fname="/work/projects/G09/SEACLID/Mask_VN_0.1/All_region.nc"

      Call rd_ncdf_2d (Fname,lon,lat,nlon,nlat,B1,"b1")
      Deallocate (lon,lat)
      Call rd_ncdf_2d (Fname,lon,lat,nlon,nlat,B2,"b2")
      Deallocate (lon,lat)
      Call rd_ncdf_2d (Fname,lon,lat,nlon,nlat,B3,"b3")
      Deallocate (lon,lat)
      Call rd_ncdf_2d (Fname,lon,lat,nlon,nlat,B4,"b4")
      Deallocate (lon,lat)
      Call rd_ncdf_2d (Fname,lon,lat,nlon,nlat,N1,"n1")
      Deallocate (lon,lat)
      Call rd_ncdf_2d (Fname,lon,lat,nlon,nlat,N2,"n2")
      Deallocate (lon,lat)
      Call rd_ncdf_2d (Fname,lon,lat,nlon,nlat,N3,"n3")

      Allocate (VN(Nlon,Nlat))
      VN=Miss
      Where (B1==1) VN=1
      Where (B2==1) VN=2
      Where (B4==1) VN=4
      Where (B3==1) VN=3
      Where (N1==1) VN=5
      Where (N2==1) VN=6
      Where (N3==1) VN=7

      Fname="Mask_VN_0.1.nc"
      Call wr_ncdf_2d (Fname,lon,lat,nlon,nlat,VN,"reg", &
         "VN climate subregions","-","lon","lat", Miss)
     
      CONTAINS

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