Use NetCDF
Implicit None
Real, Parameter :: Miss=-99.0
Real, allocatable :: Lon(:), Lat(:), Reg(:,:)

Real, allocatable, dimension (:)   :: loninfo, latinfo
Real, allocatable, dimension (:,:) :: vietnam ! thong tin cac vung khi hau: B1=2, B2=3,..., N3=8
Integer :: Nlon, Nlat, NlonAPH, NlatAPH, NtAPH
Real, allocatable, dimension (:,:,:) :: PStnAPH, TStnAPH
Real, allocatable, dimension (:,:)   :: PStn, TStn
Real, allocatable, dimension (:)     :: LonAPH, LatAPH, TimeAPH
Integer :: ilon,ilat, ilon1, ilat1

Logical :: Status
Character (Len=255) :: Fname
!-----------------------------------------------------------------------

Fname="/home/tanpv/SEACLID/MASK/VN.grid_0.25.nc"
Call rd_ncdf_2d (trim(Fname),loninfo,latinfo,nlon,nlat,vietnam,"all")
Allocate (Reg(Nlon,Nlat), Lon(Nlon), Lat(Nlat))
Reg(1:Nlon,1:Nlat) = vietnam(1:Nlon,1:Nlat); Lon(1:Nlon) = loninfo(1:Nlon); Lat(1:Nlat) = latinfo(1:Nlat)
Deallocate(loninfo, latinfo, vietnam)
Where (Reg>1) Reg=Reg-1
Allocate (PStn(Nlon,Nlat),TStn(Nlon,Nlat))  ! To be stored VN area's stations

! Precip
Fname="APH_Psta2007.nc"
Call rd_ncdf_3d (trim(Fname), LonAPH, LatAPH, TimeAPH, NlonAPH, NlatAPH, NtAPH, PStnAPH, "rstn")
Where (PStnAPH==0.OR.PStnAPH==2.0000000E+20) PStnAPH=Miss
!print*,PStnAPH;pause
Do ilon=1, Nlon
   Do ilat=1, Nlat
      Do ilon1=1, NlonAPH
         Do ilat1=1, NlatAPH
            if (LonAPH(ilon1)==Lon(ilon).AND.LatAPH(ilat1)==Lat(ilat)) then
               if (PStnAPH(ilon1,ilat1,1)>0.) then
                  PStn(ilon,ilat)=PStnAPH(ilon1,ilat1,1)/4
               else
                  PStn(ilon,ilat)=Miss
               endif
            endif
         Enddo
      Enddo 
   Enddo
Enddo

! Temp
Fname="APH_Tsta2007.nc"
Deallocate (LonAPH,LatAPH,TimeAPH)
Call rd_ncdf_3d (trim(Fname), LonAPH, LatAPH, TimeAPH, NlonAPH, NlatAPH, NtAPH, TStnAPH, "rstn")
Where (TStnAPH==0.OR.TStnAPH==2.0000000E+20) TStnAPH=Miss
!print*,TStnAPH;pause
Do ilon=1, Nlon
   Do ilat=1, Nlat
      Do ilon1=1, NlonAPH
         Do ilat1=1, NlatAPH
            if (LonAPH(ilon1)==Lon(ilon).AND.LatAPH(ilat1)==Lat(ilat)) then
               if (TStnAPH(ilon1,ilat1,1)>0.) then
                  TStn(ilon,ilat)=TStnAPH(ilon1,ilat1,1)/4 
               else
                  TStn(ilon,ilat)=Miss
               endif
            endif
         Enddo 
      Enddo 
   Enddo
Enddo
! Write to file:
Fname="APH_Num_of_Stn.nc"
Call wr_ncdf_2d (trim(Fname), Lon,Lat,Nlon,Nlat, PStn, "PStn", "Num of Prec stations", "-", "DegE", "DegN", Miss)
Call wr_ncdf_2d_apd (trim(Fname), TStn, "TStn", Nlon,Nlat, "Num of Temp stations", "-", "DegE", "DegN", Miss)

CONTAINS
Include "sub_rd_ncdf_Tan.F90"
Include "sub_wr_ncdf_Tan.F90"

END